MPOLE Programs

Functional List

Alphabetcal List: Level 1

Alphabetcal List: Level 2

all_nulls

Purpose:

try to find as many nulls as possible using a combination of
reasonable initial guess locations. remove duplicate nulls.
label the nulls Ai or Bj for negative/positive nulls
(using relabel_nulls). return an array of null structures.



Calling sequence:

nls = all_nulls( pls )


Keywords:

• noasymp: unless set a zero will be sought based on
the asymptotic field

Outputs:

• nls: the array of null structures.

Modification History:

DWL 8/01: found the assymptotic null
DWL 1/02: relabeled the nulls to be sequential
Graham Barnes 01/13/03: created most_nulls.pro
DWL 11/03: merged all_nulls and most_nulls
DWL 3/05: added call to index_check

all_seprs

Purpose:

given poles & nulls find all of the separators and
return them as an array of spr structures.
fill in structure with sepr_info



Calling sequence:

sprs = all_seprs( pls, nulls )


Keywords:

• ref_max: the number of refinements in the null-scan
• no_info: if set, exclude the calculation of the sepr_info

Modification History:

8/01: DWL converted to structure format of spr
9/02: DWL added ref_max as keyword
12/03: DWL added automatic call to sepr_info to
fill in fields in spr structure.
3/05: DWL skipped upright nulls in pot_sepr_list

cm2csv

Purpose:

given a set of poles and an array containing the connectivity
matrix, output in connectivity matrix
as a .csv (comma separated values) file.



Calling sequence:

cm2csv, pls, cm


Inputs:

• pls: the poles structure ( pls.phi & pls.lab are used )
• cm: connectvity matrix -- int matrix which is 0 for no connection
and >0 for connection.



Keywords:

• minflux=minflux: the minimum flux to include (poles
whose connections total less than this are not included)
• file=file: the name of an output file [ 'cm.csv' ]

cm_list

Purpose:

given a set of poles and an array containing the connection
matrix (returned by connectivity) print out a summary
of the connections



Calling sequence:

cm_list, cm, [ pls ]


Keywords:

• ntop=ntop: print, the ntop largest connections. [ 25 ]
• err=err: if set, this is the error matrix from connectivity
• format=format: the format for printing the fluxes

connectivity

Purpose:

calculate domain fluxes for a set of sources. uses a
Monter-Carlo method to determine which poles
are connected to which other poles. returns an array
cm(i,j)
which the approximate flux connecting i to j
i=0 refers to infinity.
the array returned is symmetric: cm(i,j) = cm(j,i)
uses a bayesian estimate of the flux connecting i to j
when m_i of N_i and m_j of N_j line connect the connecting
value with maximum liklihood is
cm(i,j) = ( m_i + m_j )/( N_i/Phi_i + N_j/Phi_j )



Calling sequence:

cm = connectivity( pls )


Keywords:

• nlines=nlines: total number of lines to initiate from each
source. [ 300 ]
• phicut=phicut: if set nlines is set to 1.5*phi/phicut
so 95! of all domains w/ psi_{ij} > phicut will be found. this overides nlines.
• nopos: if set field lines from positive poles will not
be used
• noneg: if set field lines from negative poles will not
be used
• doonly = name_list: initialize field lines from
only these poles.
• err = err: returns the statistical error in each value
• max_steps = max_steps: the maximum number of steps to take
[10000]

Modification History:

DWL - 7/03 made the indexing consistent w/ that in
domain_matrix: i=0 is infinity (used to be i=np was
infinity)
DWL - 5/04 re-wrote w/ Bayesian Max. Likelihood estimate
DWL - 4/05 added phicut keyword.
DWL - 11/05 added doonly keyword

cribbon

Purpose:

Calculate the ribbon along the separator. This is
a structure containing several (3,n) arrays.
cr.x - is the separator itself a (3,n) array
cr.w - an array of the nomralized widths
s.t. Delta = w * sqrt( |I|/I_0 )
cr.pos - array of vectors chosen to represent the
edge of a current sheet w/ current +I_*
cr.neg - array of vectors chosen to represent the
edge of a current sheet w/ current -I_*
cr.i0 - the fiducial current I_* for the separator



Calling sequence:

cr = cribbon( pls, fl ) or cr = cribbon( pls, nls, spr )


Inputs:

• pls: structure of poles
• nls: an array of separator structures
• spr: a single separator structure (refering to pls & nls )
• fl: the field line array to use instead of the spr



Modification History:

2004-04-28 DWL fixed dependency on sepr which prevented
use in mode 2.

current_solve

Purpose:

solve the equation governing the separator currents:
i * L alog( e i0 / |i| ) + M.i + psi = psi0
where M is the mutual induction matrix, psi are the separator
vacuum fluxes and psi0 is the value they are required to have.
this cannot be solved if any element of the vector
( psi0 - psi )/( i0*L ) has magntiude greater than one.
in this case psi0 is altered to make the equation solvable
and psi0 WILL BE CHANGED UPON RETURN.



Calling sequence:

i = current_solve( psi0, si [, m ] )


Inputs:

• psi0: a vector of desired fluxes. if this is
unobtainable it is changed to one which can be obtained.
• si: the structure of separator info. returned from sepr_info



Optional Inputs:

• m: the mutual induction matrix. returned by mut_induct_mat



Keywords:

• sweep: if set a solution is constructed by sweeping the
matrix m to its final value from 0.0. this will provide a unique solution related to the one w/o induction matrix.
• init=i: an initial guess for the current. this will be
disabled by the sweep keyword.
• erg=erg: will return the energy of the current system.

Modification History:

Dana Longcope 2005-07-08

disk2tan_plane

Purpose:

translate a set of poles from disk to tangent plane coordinates.
poles pls have x, y coordinates in terms of location on the
disk (radius rad) & a z coordinate giving distance above
surface (same units as x & y). translate this to coordinates
(in Mm) on a tangent plane centered at heliographic coordinates
( lat, cmd ) and disk location ( x0, y0 )
return the inverse transformation as a structure
view containing the matrix, the displacement & the
heliographic coodinates (in degrees) of the transformation
the disk location xd is related to tangent_plane location
xtp by
xd = view.mat # xtp + view.disp
xd now contains a z coordinates indicating vertical location



Calling sequence:

view = disk2tan_plane( pls )


Inputs:

• pls: - the poles. will be changed upon return



Keywords:

• rad = rad: radius of disk [ default = 925 arcsec ]
• b = b: the b-angle (deg) for heliographic coods. [ default = 0 ]
• p = p: position angle (deg) for heliographic coods. [ default = 0 ]
• center_hel = ( lat, cmd ): force center to be at there
heliographic coordinates
• center_disk = ( x0, y0 ): force center to be at this
disk location ( default = center of charge )
• date: if set this string contains the date of the
transformation
• vert: if set the fluxes in disk corrdinates are assumed to be
vertical field Bz integrated over pixels rather than line-of-sight field.
• mg = mg: a data structure which will be used to replace
keywords rad, b, p and date in one call. the structure
• must contain fields: mg.b0, mg.p, mg.rad, mg.date

Modification History:

8/01 - DWL added vert keyword
2004-04-23: DWL fixed bug finding lat & cmd from x0

dm2tex

Purpose:

given a set of poles and an array containing the domain
matrix, produce a LaTeX table (on the screen) with the
tabel.



Calling sequence:

dm2tex, pls, dm


Inputs:

• pls: the poles structure ( pls.phi & pls.lab are used )
• dm: domain matrix -- int matrix which is 0 for no connection
and >0 for connection. its value will be used to determined the pinted character.



Keywords:

• vstring = vstring: a string array w/ N elements
where N = max(dm)+1. vstring(dm(i,j)) will be printed at each location. default = [ '0', '1', '2', ... ]

dm_print

Purpose:

given a set of poles and an array containing the domain
matrix, print a table



Calling sequence:

dm_print, pls, dm


Inputs:

• pls: the poles structure ( pls.phi & pls.lab are used )
• dm: domain matrix -- int matrix which is 0 for no connection
and >0 for connection. its value will be used to determined the pinted character.



Keywords:

• vstring = vstring: a string array w/ N elements
where N = max(dm)+1. vstring(dm(i,j)) will be printed at each location. default = [ '0', '1', '2', ... ]
• col_space=cs: spacing between column. default is 2
setting to 0 will make the heading labels run together

dm_schematic

Purpose:

graph the domain matrix schematically.
the graph will be two columns (positive left of negatives)
with the lines connecting between them.



Calling sequence:

dm_schematic, pls, dm


Keywords:

• title=title: the title to print at the top
• nosum: if set the summary legent is suppressed
• linestyles=linetsyles: an array of length max(dm) st.
connection i,j is plotted with line style linestyle( dm(i,jI)-1 ) [ default is all 0s ]
• min_deg=min_deg: plot only poles with degree above min_deg
min_deg=1 omits unconnected poles
• nophi: if set the total fluxes will be suppressed
• format=format: the format for printing the fluxes

Modification History:

DWL: 2004-05-21
DWL: 2004-06-07 made character size scale w/ number of rows.

domain_char

Purpose:

calculate flux and volume of
the flux domain connecting poles ia to ib. monte carlo
integration is used to perform the integrals. field lines
are initialized at random on pole ia and traced to their
terminus.



Calling sequence:

res = domain_char( pls, ia, ib, add=add )


Keywords:

• n = n: the number of lines in the monte carlo integral [100]
• graph: if set each of the lines connecting ia->ib are rendered
• threed: same as graph, but for 3d rendering
• rad: radius to begin field lines [ 2*pls.drmax ]

domain_matrix

Purpose:

list the flux domains for a set of sources. Uses the
nulls to decide which separatrices exist. This is used
to fill in a matrix
dm(i,j) = dm(j,i)
with 1 when i connects to j and 0 otherwise.
row/column 0 stands for infinity. row/column i stands
for pole i-1



Calling sequence:

dm = domain_matrix( pls, nls, [ spr ] )


Optional Inputs:

• spr: the separator structure. if present, use
spr.anull and spr.bnull to deduce connections



Keywords:

• planar: if set use only the photospheric separatrices
• leaf: if this is set AND spr is present the leaf domains
are added using one fan trace from each null *not* connected to a separator.

Modification History:

DWL: 2004-05-21 added spr as optional input
DWL: 2005-03-22 added the leaf keyword

fill_dmn_fp

Purpose:

overplot a specific domain footprint with a filled region.
the footprint must be a simple 4-sided region with opposite
poles on diagonal vertices and prone nulls (of any type) at
the other vertices. the boundaries of the region consist
of spines and separatrix traces.
*** if the region does not fit these criteria nothing will
be plotted ***



Calling sequence:

fill_dmn_fp, pls, i1, i2, null1, null2


Inputs:

• pls: the complete poles structure
• i1, i2: indices from pls which form vertices of the region.
OR then names of the poles given as strings. these must match an entry in pls.lab exactly.
• null1, null2: the poles (complete structures) which form
the other vertices



Keywords:

• shrink=fctr: shrink the region by a factor fctr (0<fctr!1)
for aesthetic purposes.
• shade=shade: the color with which to shade the region.
• view=view: a viewing transformation to apply.

Example:

IDL> show_fp, pls, nls
IDL> fill_dmn_fp, pls, 'P03', 'N06', nls(2), nls(7)



Modification History:

5-2006: DWL

index_check

Purpose:

check poles and nulls against two Poincare
indices involving the numbers of various sources & null points:
P2d = prone - upright - sources + 2
P3d = B-type - A-type - positive + negative
if all nulls are present both indices will be zero.
if P2d is positive there are at least P2d prone nulls
missing. if it is negative, there are at least -P2d upright
nulls missing. etc.
[ see Longcope & Klapper, ApJ v.579, p.468 (2002) ]
the program prints out the numbers of each type & reports
on the conclusions derived from the indices.
if the separator is passed in it also reports the number
of domains topologically required.
domains = separators + sources - coronal-nulls - 1
[ Beveridge & Longocpe Sol Phys v.227, p.193 (2005) ]



Calling sequence:

index_check, pls, nls, [ spr ]


Inputs:

• pls: the poles structure
• nls: the array of nulls structures



Optional Inputs:

• spr: the array of separator structures



Modification History:

Dana Longcope 2005-03

inv_lam_func

Purpose:

return the inverse of function
\Lamda(u) = u ln( e / |u| ) = 12.745 f(u/12.745)



Calling sequence:

u = inv_lam_func( lambda )


Procedure:

\Lamda(x) = x ln( e / |x| ), \Lambda'(x) = - ln( |x| )
|x| = 1 are extrema, and solutions exist
only for -1 < x ! 1.
initial guess from approximation:
L(x) = y = sqrt[ x * ( 2 - x ) ]
==> x = 1 - sqrt( 1 - y^2 )

latex_poles

Purpose:

output a table in LaTeX format summarizing the poles



Calling sequence:

latex_poles, pls


Inputs:

• pls: the structure of poles



Keywords:

• view = view: is used to translate poles into heliographic
coords
• tan_plane: if set coordinates are recoreded in Mm
• pwr = pwr: write fluxes in units of 10^{pwr} Mx
• tot: if set place a total at the bottom of each side

Modification History:

DWL 12/8/00

line_info

Purpose:

return a 3 X N array (float) where each point contains
[ len, B, b^ . z^ ]



Calling sequence:

info = line_info( pls, fl )

mut_induct_mat

Purpose:

return a an array containing the mutual induction matrix
for the separators given.



Calling sequence:

m = mut_induct_mat( pls, nls, sepr )


Inputs:

• pls: poles structure
• nls: array of nulls structures
• spr: array of separator structures



Keywords:

• verbose: if set will display a graph showing the progress

Outputs:

• m: the an NxN mutual inductance matrix

Modification History:

DWL 11/07/00
DWL 11/26/03 - updated to use sepr structures

new_tan_plane

Purpose:

translate a set of poles from one tangent plane to
a new tangent plane. change pole the pole structure
and the viewing structure. the transformation is
specified by its cantact-point



Calling sequence:

new_tan_plane, pls, view, new_center [, options ] new_tan_plane, pls, view, match=view2


Inputs:

• pls: the poles. will be changed upon return.
• view: the viewing transformation - changed upon return.
• new_center: the new center float(2) specified in heliographic or
disk coordinates



Keywords:

• match = view2: advance tangent plane of this view to
the date/time of view, and use this center to do the tangent plane. if this option is set the new_center argument is not used.
• hel: if set interpret center as ( lat, cmd )
in heliographic coordinates [defult]
• disk: if set interpret center as ( x0, y0 ) in arcseconds
from center of disk
• tp: if set interpret center = ( x0, y0 ) in the tangent plane

null_graph

Purpose:

use the null & separator structures to plot a schematic graph
with pos. nulls in left column, negative nulls in right column
connected by a line for each separator



Calling sequence:

null_graph, nls, spr


Inputs:

• nls: an array of null structures
• spr: an array of separators. spr.anull and spr.bnull are
unsed to index the nls connected by each separator



Keywords:

• title=title: the title to print at the top
• nosum: if set the summary legent is suppressed
• noorphans: if set the graph will omit nulls w/o separators.

Modification History:

DWL 2004-05-21
DWL 2004-06-16 : fixed bug for npos=1 or nneg=1
DWL 2005-03-22 : changed algorithm & added multi-separator
capabilities. added noorphan keyword

pixel2tan_plane

Purpose:

translate a set of poles from to tangent plane coordinates from
coordinates given in terms of pixels relative to some arbitrary
point. the location on the disk will be inferred by specifying
(lat,cmd) for some reference pixel, ref_pixel, in the image
(assumed to be center-of-charge) and the pixel scale (assumed
to be 1 arcsecond). poles pls have x, y coordinates in terms
of location on the disk (radius rad) &
once again the z coordinate gives distance above
surface (same units as x & y). translate this to coordinates
(in Mm) on a tangent plane centered at heliographic coordinates
(lat,cmd) and pixle-location ref_pixel.
return the inverse transformation as a structure
view containing the matrix, the displacement & the
heliographic coodinates (in degrees) of the transformation
the pixel-location xp is related to tangent_plane location
xtp by
xp = view.mat # xtp + view.disp
xp now contains a z coordinates indicating vertical location



Calling sequence:

view = pixel2tan_plane( pls, lat, cmd )


Inputs:

• pls: the poles. will be changed upon return
• lat: heliographic latitude of ref_pixel
• cmd: heliographic longitude of ref_pixel



Keywords:

• rad = rad: radius of disk [ default = 925 arcsec ]
• b = b: the b-angle (deg) for heliographic coods. [ default = 0 ]
• p = p: position angle (deg) for heliographic coods. [ default = 0 ]
• ref_pixel: 2 element array containin the pixel to which
lat and cmd refer
• pix_scale: size of pixels in arcseconds

pls_phi_hist

Purpose:

plot a histogram summarizing the fluxes in all poles
a particular structure.



Calling sequence:

pls_phi_plot, pls


Inputs:

• pls: the poles stucture to be summarized. pls.phi is
used (and possibly scaled) and pls.lab to label the bins



Keywords:

• scl=scl: multiply fluxes by scl. [ 1.0 ]
• mx: if set then scl=1.0e16 meaning the fluxes are assumed
be in G Mm^2 and the plot will be in Mx.
• phimax=phimax: the range of the plot.
• nlab=nlab: label the largest nlab poles along the bottom. [5]
• nmax=nmax: include only nmax largest positive and nmax larges
negative sources. default is to plot all of them.
• title=title: title for the top of the plot.

Modification History:

2005-06-22 Dana Longcope

polka_map

Purpose:

plot locations of each pole in a poles structure as
disks whose area is proportional to their flux.



Calling sequence:

polka_map, pls


Keywords:

• over: if not set a bounding box will be constructed.
and a grey-scale will be rendered
• nolab: omit labels
• view = view: will transform according to transformation in
structure view. title=title
• maxrad = maxrad: if set this will define the maximum
radius in data units [ default = 0.1*xrange ]
• flux_scl=flux_scl: the value of the flux to assign the maximum
radius [ default is maximum flux ]
• axis=ax: if set to 1 no axis will be draw
if set to 2 a box will be drawn

rd_pns

Purpose:

read in a .pns file containg poles, nulls & separators.
return all of these as the appropriate structures/arrays.
a file may contain only poles, or poles & nulls. and the
remaining structures will not be set upon return. old-style
.mp files, containing only poles may be read in.
a line @--- will separate multiple entries in a single file.
if keyword multi is set then read through multi-1 of these
lines.



Calling sequence:

rd_pns, file_name, pls [ , nulls [ , sepr ] ]


Keywords:

• view = view: will return the viewing transformation (if
present)
• dm = dm: will return the domain matrix (if present)
• multi: is set read through multi-1 different entries
if there are not this many entries return multi set to -1.
• norefine: if set the null positions will not be refined
• info: if set the separator information will be stored
in the sturcture. this takes additional time.
• norelabel: if set the nulls will not be relabeled
• noends: if set the end-points of spines will not be found.
• rmax = rmax: (fix rmax)

Modification History:

1/02 - DWL call to relabel_nulls to allow sequential labeling
11/14/03 - DWL fixed bug which crashed program if nulls was
not passed and separators exisited.
5/24/04: added dm option
9/09/04: added call to sepr_info to fill in spr fields

sepr_info

Purpose:

return an array (or structure)
containing magnetic information on each of
the magnetic separators in the array sepr.
the information: [ length, flux, I_0, s, hel ]
I_0: I_*/c in units of flux/length
For Mx & cm I_0 has units of Mx/cm = 10 Amps
(i.e. 10*I_0 is the current in Amps)
For 1.0e16 Mx & Mm I_0 has units of 1.0e9 Amps
hel: relative helicty per current. The relative helicity
of the entire field is the sum over
all separators as H_R = \sum_i hel(i) cur(i)



Calling sequence:

info = sepr_info( pls, nulls, sepr )


Keywords:

• structure: if set the information is handed back as
a structure, rather than an array.
• lineclose: if set the flux is calculated using
a return path making a straight line between the null points

Procedure:

The self flux is expressed in the form
\Psi = (I/c) L \ln[ e I_0 / |I| ]
= \Psi_0 ( I / I_0 ) \ln[ e I_0 / |I| ]
= \Psi_0 \Lamda( I/I_0 )
\Psi_0 = (I_0/c) L and \Lambda(x) = x \ln( e/|x| )
( Lambda(x) is inverted by inv_lam_func.pro )
where L is the length of the separator and \Psi refers to the flux
through the separator. (these are each 1/2 the value for
a circular loop for which s = 1) in general
\ln(s) = {1\over L}\int d\ell \int d\ell'\left[
{{\bf\hat{n}}\cdot{\bf\hat{n}}' \over |{\bf x}-{\bf x}'|}
-{C(\ell-\ell') \over L}\right]
where C(\elll-\elll') is the normalized term for a circular loop
C(x) = \pi \csc(\pi|x|/L) - 2\pi\sin(\pi|x|/L)
C(x) = {\pi \cos(2\pi|x|/L)\over\sin(\pi|x|/L)}
and the fiducial current is defined
I_0/c = s { 64 e^{-4} \over pi^2 } \langle B'\rangle L^2
where \langle B' \rangle is geometric average
\langle B'\rangle=\exp\{ L^{-1} \int \ln[ B'(\ell) ] d\ell \}
and B'(\ell) is returned by bp_eval.pro



Modification History:

10-10-02 - DWL rescaled the definition of i0.
03-10-03 - DWL added the helicity piece which calls eval_z
12-15-05 - DWL mutliplied helicity by 4pi, so it
is correct in cgs-emu (as is everything else)

show_3dnulls

Purpose:

graph all the nulls in an array of null-structures



Calling sequence:

show_3dnulls, nulls


Keywords:

• label: if set, display the labels of each null

show_3dpoles

Purpose:

plot locations of each pole in a poles structure --- in 3d



Calling sequence:

show_3dpoles, pls


Keywords:

• xrange = [ x0, x1 ]: force the xrange of plot window
• yrange = [ y0, y1 ]: force the yrange of plot window
• zrange = [ z0, z1 ]: force the zrange of plot window
• over: if not set a bounding box will be constructed.
• nolab: omit labels
• zaxis = za: the zaxis will go from 0 to za
default is the full zrange
• ax = ax: the x-rotation for the 3d projection
• az = az: the x-rotation for the 3d projection

show_connectivity

Purpose:

graph the connectivity on an exisiting graph of poles



Calling sequence:

show_connectivity, pls, c_mat


Keywords:

• min_flux = min_flux: the minimum flux to indicate [ 0.5 ]
• count = count: return the number of connection
• inf: if set, connections to infinity are depicted as arrows
linestyle=linestyle

Modification History:

DWL: 7/03 made the indexing consistent w/ that in
domain_matrix: i=0 is infinity (used to be i=np
was infinity)
DWL: 12/03 added \inf flag (modeled on show_dm)

show_cr_skel

Purpose:

draw the skeleton of a current ribbon over an exisiting plot



Calling sequence:

show_cr_skel, cr


Inputs:

• cr: structure containing the current ribbon as returned by cribbon



Keywords:

• i = i : the current (unormalized)
• norm_i = ni: the normalized current i/i_*
• view = view: will transform according to transformation in
structure view.
• threed: show field line in 3d
• ribs=n: if set the "ribs" of the skeleton are shown.
show every n ribs

show_dm

Purpose:

graph the domain matrix on an exisiting graph of poles



Calling sequence:

show_dm, pls, dm


Keywords:

• linestyle = linestyle: linestyle used to l]plot lines
• mindegree = md: will only show edges connecting to
vertices of degree md or greater
• inf: if set arrows are used to show linkages to infinity

Modification History:

DWL: 12/03 - fixed so inf-inf is not rendered

show_domain

Purpose:

show field lines from the domain ia/ib. nlines field lines
are started at ia & traced until they hit another pole.
only those which terminate at ib are shown.



Calling sequence:

show_domain, pls, ia, ib


Keywords:

• ntries = ntries: the number of lines in the monte carlo
integral [100]
• nlines = nlines: the number of lines to draw. if
unspecified then every one fromthe MC integral is drawn
• threed: same as graph, but for 3d rendering
• rad = rad: radius to begin field lines [ 2*pls.drmax ]
• view = view: project with view transfrom
• max_segs = max_segs: set to the maximum number of segments in a line
[ 1000 ]

show_flux_tube

Purpose:

graph a set of field lines which start in a ring about the
point x0. the ring has radius r and is oriented normal to the
magnetic field at point x0. there are nlines field lines in
the ring



Calling sequence:

show_flux_tube, pls, x0


Inputs:

• pls: the poles structure
• x0: 3-element array (float) of the initial point.



Keywords:

• nlines=nlines: number of lines to use [ 30 ]
• radius=radius: radius of circle [ 10*pls.drmax ]
• threed: if set, render in 3d
• view=view: use view structure
• circle: if set draw the circle around x0
• clip: if set clip the field lines
• dir: if set to +1 ot -1 only that direction will be traced

show_fp

Purpose:

graph the locations of sources, nulls and footprints of
the field. footprints are shown as dashed and solid lines.
solid lines denot photospheric spines, dashed lines are
the footprints of separatrices. dotted lines
are the spines of coronal null points.



Calling sequence:

show_fp, pls, nls


Keywords:

• view = view: will transform according to transformation in
structure view.
• nosepx: do not show the separatrix field lines
• nogamma: do not show the gamma-lines (spines)
• view = view: perform the viewing transformation
• threed: show the field lines in 3d
• clip: if set, each of the lines in clipped
• max_segs: set to the maximum number of segments in a line
[ 1000 ]
• label: label the nulls
• axis=ax: if set to 1 no axis will be draw
if set to 2 a box will be drawn
• over: plot is made over the exisiting plot
• color: if set the pos/neg poles are shown in white/black
this is useful when plotting over a grey-scale magnetogram
• mg: setting this is equivalent to /over, /color
which is useful when overplotting on a magnetogram

Modification History:

5-2004: Written by Dana Longcope
1-2004: DWL fixed bug in viewing transform

show_gamma_lines

Purpose:

graph all the gamma lines for a set of nulls



Calling sequence:

show_gamma_lines, pls, nulls


Keywords:

view = view - will transform according to transformation in structure view. nosepx - do not show the separatrix field lines nogamma - do not show the gamma field lines view - perform the viewing transformation threed - show the field lines in 3d clip - if set, each of the lines in clipped max_segs - set to the maximum number of segments in a line [ 1000 ]

show_loop_skeleton

Purpose:

given a field line and a set of radii, draw a sekeloon
consisting of rings arounf each point in the field line



Calling sequence:

show_loop_skeleton, fl, rad


Keywords:

• view = view: will transform according to transformation in
structure view.
• threed: show field line in 3d
• skip=skip: if set then only every skip ribs are rendered. [1]

show_nulls

Purpose:

graph all the nulls in an array of null-structures



Calling sequence:

show_nulls, nulls


Keywords:

• label: if set, display the labels of each null
• new: draw new set of axes.
• view = view: will transform according to transformation in
structure view.

Modification History:

DWL 2004-06-17: adjusted label locations for prone nulls

show_poles

Purpose:

plot locations of each pole in a poles structure



Calling sequence:

show_poles, pls


Keywords:

• xrange = [ x0, x1 ]: force the xrange of plot window
• yrange = [ y0, y1 ]: force the yrange of plot window
• over: if not set a bounding box will be constructed.
• color: if set the pos/neg poles are shown in white/black
this is useful when plotting over a grey-scale magnetogram
• nolab: omit labels on the poles
• axis=ax: if set to 1 no axis will be draw
if set to 2 a box will be drawn
• view = view: - will transform according to transformation in
structure view.
• mg: setting this is equivalent to /over, /color
which is useful when overplotting on a magnetogram

Modification History:

2004/04/09: DWL made plots isotropic

show_random_3dlines

Purpose:

graph nlines different field lines



Calling sequence:

show_random_3dlines, pls


Keywords:

• nlines=nlines: the number of lines to draw. [ 30 ]

show_sepr

Purpose:

draw all the separators in array sepr



Calling sequence:

show_sepr, pls, nls, spr


Keywords:

• view = view: will transform according to transformation in
structure view.
• threed: show field line in 3d
• shadow: if set (& in 3d mode) a projection onto z=0 is also
shown
• clip: clip the field line

show_sepx

Purpose:

graph field lines from the separatrix of null.



Calling sequence:

show_sepx, pls, null


Inputs:

• pls: the poles structure
• null: a single null structure



Keywords:

• view = view: set to the viewing transformation
• thrnge = [ th0, th1]:
a 2-element array defining the range of theta to be used. default is [ null.theta0, null.theta0 + !pi ]
• thvals = thv: an array of theta values
• nlines = n: the number of lines to draw. default is 30
• threed: if set the field lines are rendered in 3d.
• clip: if set this will clip the field lines

solar_rotate_view

Purpose:

given a structure for viewing at a tangent plane. advance the
view by a fixed time interval (in seconds) according to
the solar rotation. the center of the
of the tangent plane is ( lat, cmd ), b angle and p angle
will be given in degrees. the view structure contains the
following elements
mat: the 3 x 3 transformation matrix
disp: the displacement vector. transformation from
tan-plane coordinates xtp to disk coordinates xd:
xd = view.mat # xtp + view.disp
b: the b angle (degrees)
p: the p angle (degrees)
rad: apparent solar radius (arcseconds)



Calling sequence:

new_view = solar_rotate_view( view, time )


Inputs:

• view = view: the old view
• time: time (in seconds)



Keywords:

• to = date: if set advance the view to a date specified
by the string date from view.date
• rad=rad: this new radius will be used rather thant the one
in view. this permits a change in perspective as between SOHO and an Earth-orbiting satellite
• b0=b0: this new b0 will be used rather thant the one
in view. this permits a change in perspective as between SOHO and an Earth-orbiting satellite
• shift=[dx,dy]: shift the view by a [dx,dy] arcseconds.
this permits co-alignment

Modification History:

2004-07-12 DWL: added rad and b0 keyords to permit
changes between SOHO and Earth-orbiting S/C

sum_graph

Purpose:

draw the domain graph and null graph for a set of poles, nulls
and separators



Calling sequence:

sum_graph, pls, nls, spr, [ dm ]


Inputs:

• psl: the poles structure
• nls: an array of null structures
• spr: an array of separators. spr.anull and spr.bnull are
unsed to index the nls connected by each separator. this is used to generate the domain matrix if it is not passed in separately



Optional Inputs:

• dm: the domain matrix



Keywords:

• title=title: a title for the top of the plot
• noorphans: if set the null graph will omit nulls w/o separators.
• leaf: if keyword is set, and spr is present, use the slower
domain_matrix algorithm which finds leaf domains. from nulls w/ unbroken fans.

Modification History:

DWL: 2004-05-21

tan_plane2disk

Purpose:

translate a set of poles from the tangent plane to disk coordinates
(arcseconds from disk center). the transoforamtion is
specified in the structure view. returns a new viewing
transformation which does the composite transformation.



Calling sequence:

new_view = tan_plane2disk( pls, view )


Inputs:

• pls: the poles. will be changed upon return
• view: the structure specifying the transformation

view_xform

Purpose:

returns a structure for viewing at a tangent plane. the origin
of the tangent plane is ( lat, cmd ), b angle and p angle
will be given in degrees. the view structure contains the
following elements
mat: the 3 x 3 transformation matrix
disp: the displacement vector. transformation from
tan-plane coordinates xtp to disk coordinates xd:
xd = view.mat # xtp + view.disp
b: the b angle (degrees)
p: the p angle (degrees)
rad: apparent solar radius (arcseconds)
scale: the linear scale factor.
date: a string containing the date & time of the
observation



Calling sequence:

view = view_xform( lat, cmd )


Inputs:

lat - the latitude cmd - the longitude



Keywords:

• rad = rad: radius of disk [ default = 925 arcsec ]
• b = b: the b-angle (deg) for heliographic coods. [ default = 0 ]
• p = p: position angle (deg) for heliographic coods. [ default=0 ]
degrees - if set lat & cmd are taken to be in degrees date = date_str - the string containg date & time

Modification History:

2000-11-15 - DWL fixed p-angle
2004-04-26 - DWL made viewing transform a named structure
of type "view"

voronoi_pls

Purpose:

a graph of the voronoi tesselation formed
by a set of poles



Calling sequence:

voronoi_pls, pls


Keywords:

• fill: if set the positive regions will be filled using
the default line-fill pattern.

Modification History:

DWL 2006-02-26: adjusted label locations for prone nulls

write_pns

Purpose:

output a file containing poles, nulls & separators



Calling sequence:

write_pls, file_name, pls [ , nls [ , sepr ] ]


Inputs:

• file_name: string the name of the file to be created.
file_name = '-' means write to std out
• pls: the structure containing all poles.
• nls: array of structures containing the nulls
• sepr: float array ( 3, n ) contaning the separators



Keywords:

• view = view: the viewing transformation.
• append: if set the information is added to the end
of the file, and can be read using the multi option of rd_pns
• dm=dm: the domain matrix

Modification History:

2004-06-11 DWL: made a check so that a sepr=0.0 will not crash
the program

xadd_nulls

Purpose:

put up a window displaying poles and nulls (and possibly
gamma lines), and add new potential nulls.



Calling sequence:

xadd_nulls, pls, nls


Inputs:

• pls: the list of poles
• nls: the null array. will be changed on output.



Keywords:

• z=z: put the initial guess at his level [ default: z=0 ]
• view=view: a viewing transformation

Modification History:

1/15/2003: DWL - added viewing xform capability

xselect_fl

Purpose:

on a window showing poles etc. use the mouse to select
field lines at various heights.



Calling sequence:

xselect_fl, pls


Inputs:

• pls: (optional) the list of poles



Keywords:

• z0 = [ zlist ]: the values of z to use for the possible
field lines. default is [10,20,30,40]. one field line will be integrated for each height.
• view=view: nulls are displayed using viewing trasnform view.
• quiet: if set the lines are plotted w/o listing the initial points.

Modification History:

11/04 - DWL

xselect_null

Purpose:

put up a window displaying poles and nulls (and possibly
gamma lines), and select a subset of the nulls with the
cursor. return only that subset.



Calling sequence:

new_nulls = xselect_nulls( nls [, pls ] )


Inputs:

• nls: the complete list of nulls (from which to select
• pls: (optional) the list of poles



Keywords:

• gamma: draw the gamma lines (spines)
• over: nulls are already displayed
• omit: rather than selecting nulls to include, select
nulls to omit
• view=view: nulls are displayed using viewing trasnform view.

Modification History:

3/02 - DWL - added viewing xform capability

advance_dif_rot

Purpose:

advance cmd according to differential rotation law
\[ \Omega = a + b \sin^2(lat) \]
a = 14.552, b = -2.84 degrees per day (sidereal)



Calling sequence:

advance_diff_rot, lc, time


Inputs:

• lc: [ lat, cmd ] in radians unless degree flag is set
• time: time (in seconds)
• KEYWORD PARAMETRS:
• degree: if set the lat and cmd are interpretted in degrees
(default is radians)

asym_fl_map

Purpose:

given a point at or beyond the asymptotic field limit
|r1| > pls.rmax, return the point at the same radius
which the asymptotic field maps to. use the multi-pole
expansion in pls. if the field line does not return
then closes = 0 upon return



Calling sequence:

r2 = asym_fl_map( pls, r1, closes=closes )


Keywords:

closes - set to 0 if the field does not close

Procedure:

The external field can be written B = -grad(chi)
where chi(x) has only the monopole & dipole terms
Q p.x
chi(x) = --- + ---
|x| |x|^3
the field is axisymmetric about an axis || to p. using
coordinates where p = |p| z^2 and the axi-symmetric flux
function f(r,theta), s.t. B = grad(f) X grad(phi)
f(r,th) = - Q cos(th) + |p| sin^2(th)/r
two points ( R, th1 ) & ( R, th2 ) are on the same field line
if f(R,th1) = f(R,th2)
0 = -Q [ cos(th1) - cos(th2) ] + |p|/R [ sin^2(th1) - sin^2(th2) ]
= (QR/|p|) [ cos(th2) - cos(th2) ] + [ cos^2(th2) - cos^2(th1) ]
= [ cos(th2) - cos(th1) ] * [ (QR/|p|) + cos(th2) + cos(th1) ]
the solutions are cos(th2) = cos(th1) ==> th2 = th1 (trivial)
Q R
cos(th2) = - --- - cos(th1) [1]
|p|
if Q = 0, the asymptotic field is purely dipolar: th2 = pi - th1
if R > |2p/Q|, the distance to the null in the assymptotic
field, then no solution th2. in this case all field lines
crossing |x| = R will map to infinity.
otherwise, there are a range of angles th1 for which a solution
th2 exists. for Q > 0 this will be a null point at r = 2|p|/Q
along the th = pi direction. to be solvable the rhs >= -1
cos(th1) <= 1 - QR/|p| ==> th1 >= acos( 1 - QR/|p| )
and 0 <= th1 != pi since it is a polar angle.
note that if R = 2|p|/Q th1 = pi is the only angle permitted,
and it maps to th2 = pi.

avg_bp

Purpose:

return <B'>



Calling sequence:

ab = avg_bp( pls, fl )

calc_sepr_info

Purpose:

calculate the properties of each separator and add it
to the spr structure. if there are multiple separators
do all of them.



Calling sequence:

calc_sepr_info, pls, nulls, sepr

clip_fl

Purpose:

clip a field line to fit in 3d window



Calling sequence:

fl = clip_fl( fl )

closest_pole

Purpose:

return the index of the pole closest to the point x



Calling sequence:

i = closest_pole( pls, x )

create_null

Purpose:

locate a null-point, given an initial guess. return a
null structure null:
null.x - float array(3) the x,y,z coordinares of the null
null.b - float |B(x)|
null.lam - float array(3) the eigenvalues of bprime (sorted)
null.e0 - float array(3) the spine direction
null.e1 - float array(3) one of the span direcitons
null.e2 - float array(3) other spine direction
null.type - character, either 'A' or 'B'
null.char - character, 'P' (prone), 'U' (upright), 'C' (coronal)
of 'M' (mirror)
null.ends - int array(2) the indices of charges at ends of
the spines. -1 for infinity. sorted in
ascending order
null.label - string giving the name of the null
null.theta0 - float the angle begining the circle on the
plane. after that circle increases by !pi



Calling sequence:

null = create_null( pls, x ) KEYWORD: noend - if set the spines will not be traced


Modification History:

8/5/99 - DWL made coordinate unambiguous. both set_null_ends
and set_theta0 can chenge directions of e0, e1 & e2
8/01 - DWL accounted for mapping of asymptotic field
1/02 - DWL changed labeling to use relabel_nulls
5/03 - DWL added char field

create_pls

Purpose:

create a blank poles structure, with np poles.
pls.x - float array (np) of x-locations
pls.y - float array (np) of y-locations
pls.z - float array (np) of z-locations
pls.q - float array (np) of charges of the poles
pls.phi - float array (np) of fluxes = 2*!pi*pls.q
pls.lab - string array (np) of labels
pls.rmax - float --- characterisitic maximum distance
pls.bmax - float --- characterisitic maximum field
strength
pls.drmax - float --- characterisitic small distance
pls.coc - float array (3) is the
center-of-(unsigned)charge
pls.mpole - float array with the monopole moment
of set of poles
pls.dpole - float array (3) with the dipole moment
of set of poles relative to the COC



Calling sequence:

pls = create_pls( np )

create_spr

Purpose:

create a blank separator structure
spr structure:
spr.anull - (int) index of the A-type null
spr.bnull - (int) index of the B-type null
spr.atheta - (float) angle from the A-null
spr.btheta - (float) angle fromt he B-null
spr.poles - (int array 4) list of poles ++--
spr.label - (string) the label for the separator
fields filled in by spr_info
spr.psi - the flux enclosed between separator and z=0
spr.len - the length of the separator
spr.zmax - maximum height of separator
spr.i0 - the characteristic current



Calling sequence:

spr = create_spr( [ array ] )


Inputs:

• array (optional): the 4-element float array
used in v1.0. this will be used to fill in anull, bnull, atheta and btheta fields of structure



Modification History:

8/01 - DWL

eval_a

Purpose:

at a point, x(0:2) evaluate the vector potential due to
the set of poles given by pls. for a point charge at the
origin the vector potential is
1 - cos(th) ( 1 - z/r )
A = ---------- phi^ = ----------- z^ X r
r sin(th) x^2 + y^2



Calling sequence:

a = eval_a( pls, x )

eval_b

Purpose:

at a point, x(0:2) evaluate the magnetic field due to
the set of poles given by pls



Calling sequence:

b = eval_b( pls, x )

eval_b2p

Purpose:

at a point, x(0:2) evaluate the 2nd derivative tensor of
the magnetic field due to the set of poles given by pls
d^2 B_i d^3
d(i,j,k) = --------- = - -------------- chi(x)
dx_j dx_k dx_i dx_j dx_k



Calling sequence:

b = eval_b2p( pls, x )

eval_bprime

Purpose:

at a point, x(0:2) evaluate the derivative matrix of
the magnetic field due to the set of poles given by pls



Calling sequence:

b = eval_b( pls, x )

eval_z

Purpose:

at a point, x(0:2) evaluate the Auxilliary field Z(x) s.t.
curl(Z) = A the vector potential.
for a point charge at the origin the vector potential is
1 - cos(th) ( 1 - z/r )
A = ---------- phi^ = ----------- z^ X r
r sin(th) x^2 + y^2
1 - cos(th) r z - r^2 z^
Z = ---------- theta^ = --------------
sin(th) |r|( |r| + z )



Calling sequence:

z = eval_z( pls, x )

find_null

Purpose:

Newton-Raphson to find a null, given a guess. Return the
converged null in place of the guess



Calling sequence:

find_null, pls, x


Keywords:

twod - restricts search to x,y plane stat - return the status of the search. stat=0 if failed

Modification History:

DWL - 7/01 check changes in location for convergence.
not absolute magnetic field

fl_from_point

Purpose:

return a single field line, traced from from point x0.
default is to plot in both directions.



Calling sequence:

fl = fl_from_point( pls, x0 )


Keywords:

• max_segs: set to the maximum number of segments
• seg_size: number of integration steps per segment [ 10 ]
• dir: if set to +1 ot -1 only that direction will be traced

fl_intgrt

Purpose:

integrate the field line from initial point x0
until termination by one of several criteria



Calling sequence:

done = fl_intgrt( pls, x0 )


Inputs:

pls - the structure of poles x0 - the initial point for the field line will be change on return
• RESULT:
done has the following value 0 = failure 1 = len_max is exceeded 2 = cang_min exceeded 3 = bmax exceeded 4 = plane crossed 5 = dist_max exeeded 6 = rad_max exeeded



Keywords:

dir - trace in direction dir (default is +1.0 ) max_step - the maximum integration step to take len_max - maximum integrated distance of line dist_max - maximum distance from x0 step_max - maximum size of a step (default = pls.drmax) max_steps - maximum number of steps (default = 1000) cang_min - the minimum cosine of the normal angle cross = cpv = float array(3,2) = [ [ xctr ], [ dir ] ] if set end the tracing when the field line crosses the plane defined by point xctr and direction vector dir

fl_view_xform

Purpose:

given a field line perfom the specified viewing transformation.
the field line is returned transformed.



Calling sequence:

fl_view_xform, fl, view

gamma_line

Purpose:

return one of the gamma lines from null, including its two ends



Calling sequence:

fl = gamma_line( pls, null )

init_sepx_line

Purpose:

begin a field line at position theta along the given null.
locally the field has the structure
B = lam1 (r.e1) e1 + lam2 (r.e2) e2
where |lam1| < |lam2| are the eigenvalues of bprime
the field lines can then be defined in terms of
a = lam1/lam2, so that 0 <= a != 1
r = ( cos(th) )^a rho^a e1 + sin(th) rho e2
where rho increases outward from a field line & is the largest
physical distance at an
d r / d rho = rho^{-1} [ a (r.e1) e1 + (r.e1) e2 ]
= ( rho lam2 )^{-1} B



Calling sequence:

x = init_spex_line, pls, null, theta


Modification History:

3/9/02 - DWL fixed pathelogical situation when theta = 0.5*!pi
3/26/02 - DWL fixed pathelogical situation when theta = -!pi
5/25/03 - DWL tinkered w/ a for cases of small ratios

is_pole

Purpose:

return the index of the pole identified by label lab
if no pole matches return -1.
print, pls.lab( is_pole( pls, 'N01' ) )
N01



Calling sequence:

i = is_pole( pls, lab )

line_piece

Purpose:

return a section of a field line, with a fixed number of
points, uniformly distributed. input and output are float
arrays fl(3,*)



Calling sequence:

fl = line_piece( ofl, zmin=zmin, n=n ) ARGUUMENTS: ofl - the old field line


Keywords:

zmin - the lowest altitude for the field line n - the number of points in the field line ref = ref. fi non-zero the points will be exponentially concentrated at the ends with maximum ratio ref (> 1.0)

null_incidence_matrix

Purpose:

return an array listing the number of separaors connecting
a given pair of null points.
mat(i,j): the number of separators linking null i to null j



Calling sequence:

mat = null_incidence_matrix( nls, spr )


Inputs:

• nls: the array of null structures
• spr: the array of separator structures



Keywords:

• count=count: an array where count(i) gives the number
of nulls incident on null i

Outputs:

• mat(i,j): the number of separators linking null i to null j

Modification History:

Dana Longcope 2005-03-22

null_scan

Purpose:

given a single null compile a list of all charge-pairs (by index)
contacted by the null. if the keyword th is set, it will be
returned with a list of all angles theta bracketing the
separator.;+



Calling sequence:

inds = null_scan( pls, null )


Inputs:

pls - the structure of poles null - a single null structure



Keywords:

twopi - scan over all 2pi radians, rather than simply pi nth - the number of segments to check ref_max - the number of times to refine a boundary [ 5 ]

Outputs:

• inds - intarr(2,n): lists the pairs of poles bracketed
by separators.

Modification History:

9/26/02 - DWL added keywords: nth & ref_max.
converted angles to double-precision
5/24/03 - DWL if null.char eq 'C' do a 2pi scan

pole_of_theta

Purpose:

return the pole index found from integrating away from a given
null at angle theta. If the trace extends past the afl radius
continue using asym_fl_map.



Calling sequence:

i = pole_of_theta( pls, null, th_val )


Inputs:

• pls: the poles structure
• null: a single null structure
• th_val: the angle from which to sent a fan field line
from the null point.

random_field_line

Purpose:

return a single field line, selected at random.
probability is based on flux.



Calling sequence:

fl = random_field_line( pls )

relabel_nulls

Purpose:

define the null labels according to a given scheme.
sequential:
nls(5).lab = 'A06' or 'B06'
polar:
nls(5).lab = 'P03-P15' etc.



Calling sequence:

relabel_nulls, nls, [ pls ], [ keywords ]


Keywords:

• sequential: use sequenctial labeling [ default ]
• polar: use polar labeling
• sort: put the nulls into an order based on spine sources

Modification History:

2004-06-22 DWL: fixed bug in /polar setting

relabel_poles

Purpose:

assign pole labels in sequential order
pls.lab(5) = 'P06' or 'N06'



Calling sequence:

relabel_poles, pls,


Keywords:

• sort: if set the poles are sorted in decreasing order
of magnetitude order = order - if initially set to an array the array will be returned with values in the new order.

sepr_line

Purpose:

return a separator field line from nullA - nullB. defined by
the angle thA and thB



Calling sequence:

fl = sepr_line( pls, nulls, sepr )

sepr_plane

Purpose:

return a point, and a triad of vectors defining a plane normal
to the separator. this is returned as a 3 X 4 float array.
result(*,0) is the origin of the new coordinate system.
result(*,1) and result(*,2) are vectors in the normal plane
and result(*,2) is the normal vector (parallel to B).



Calling sequence:

result = sepr_plane( pls, nulls, sepr )


Modification History:

2004-04-24: DWL increased the maximum number of steps

sepx_line

Purpose:

return a separatix field line from null along direction theta



Calling sequence:

fl = sepx_line( pls, null, theta )


Keywords:

max_segs=max_segs - the maximum number of segments to use

tan_plane_matrix

Purpose:

a 3 x 3 matrix translating neighbors of the point ( lat, cmd )
on the disk, to a tangent plane.



Calling sequence:

tpmat = tan_plane_matrix( lat, cmd ). matrix takes a tangent vector in disk coordinates [ vx, vy, vz ] = tmat # [ vw, vn, vr ]


Inputs:

lat - latitude of point cmd - the longitude



Keywords:

• b = b: the b-angle (deg) for heliographic coods. [ default=0 ]
• p = p: position angle (deg) for heliographic coods. [default=0 ]
degrees; if set the lat & cmd are in degrees (default is radians)

update_nulls

Purpose:

given a set of poles and a set of nulls from an related set of
poles update the nulls. use their location as initial guesses.



Calling sequence:

new_nulls = update_nulls( pls, nls )


Modification History:

3/02 - DWL relabel the nulls w/ sequential labels

view2time

Purpose:

from a single viewing trasnformation or an array return
a float (or array of floats) denoting the time from
some reference point. unless otherwise specified the time
will be in hours and the reference time will be the begining
of the day of the first element.



Calling sequence:

time = view2time( view )


Inputs:

• view: a view-xform structure, or array of structures.



Keywords:

• days: if set the time is converted to decimal days
• from=datestr: a string giving the reference tim

Outputs:

• time: a float or float array

Modification History:

2005-04-26