Calling Sequence:

Optimization_fff, bx_in, by_in, bz_in, $
                  filename = filename, $
                  x_in = x_in, y_in = y_in, z_in = z_in, $
                  Restart_flag = restart_flag, $ 
                  Bc_flag = bc_flag, $
                  potlin = potlin, iterations = iterations, $
                  Dt = dt, nsave = nsave, double = double, $
                  convergence_level = convergence_level, $
                  pr_iter = pr_iter, $
                  alpha_in = alpha_in, $
                  slow_way = slow_way, $
                  spherical = spherical, $
                  dont_kill_input = dont_kill_input
                  use_fortran = use_fortran, $
                  fortran_exe = fortran_exe, $
                  use_sdf = use_sdf,$
                  wfile=wfile, weight_in=weight_in, wtype=wtype

Input:

bx_in, by_in, bz_in = the input x, y, and z (or r, theta and phi) components of the magnetic field. Each array can be two or three dimensions; two dimensional input arrays are considered to be the fields at z=0, the surface of the sun. Three dimensional input fields define the full extrapolation volume.

Keyword Input:

filename = The name of a file that contains the input magnetic field. For this option, then input B arrays must be 3D. The file should be an ascii file, with format: xsize ysize zsize bx by bz where xsize, ysize, and zsize are the grid sizes of the x, y, and z directions, bx, by, and bz are the 3 components of the field. This option is useful for restarting calculations.

x_in, y_in, z_in = input spatial grids. If the input B's are 3D, then the number of elements in each input grid must match the appropriate size in the B array:e.g., If you input 128x64x32 B fields, then x_in must be 128 elements, y_in must be 64 elements, and z_in must be 32 elements. If you pass in 2D arrays for B, then only x_in and y_in must match. The default for grid creation is to use uniform size grids in all 3 directions, with the largest transverse grid (x or y) evenly divided between -1.0 and 1.0, and equal grid spacing for all three directions: E.g., in the 128x64x32 example, the x grid will stretch from -1 to 1 in 128 steps, the y grid will be from -0.5 to 0.5 in 64 steps and the zgrid will be from 0 to 0.5 in 32 steps. The default is to use uniformly spaced grids, but the user may pass in any grids. Note that, while it is not necessary to pass in all three grids if you are passing in custom grids, it is a good idea to do so, especially if the grids do not range from 1 to -1 (or 0 to 1 for z). Updated on 3-Nov-2008.

If the input is 2D, then the z grid is chosen to be the same size as the minimum transverse grid, e.g., for 128x64 input, the 3D fields will be 128x64x64, with the xgrid from -1.0 to 1.0 in 128 steps, the y grid from -0.5 to 0.5 in 64 steps, and the z grid from 0 to 1 in 64 steps. Note that the code has been changed from the original version in order for this to be true, Updated on 28-Nov-2005.

restart_flag = if set, then the program will read in the B field from a file (You must use the filename option), and use this field as the starting point.

bc_flag = if set, then the B field on all boundaries is set to the boundary value for the input field. (This means that the input field must be 3D.) The boundary field remains unchanged for the entire calculation. If not set (or set to zero), then only the bottom boundary is used for the calculation, and the field on the upper and side boundaries is given by the initial potential or linear force-free approximation. Note that the boundary field remains unchanged for the entire calculation for this case too. The default is 0, lower BC's only.

potlin = If 0 or unset, then the initial field is a potential extrapolation (from Tom Metcalfs FFF.pro). If set to 1 then the initial field is set to a linear FFF extrapolation. If set to 2 to the the field is linear FFF near the origin, but varies gradually to be potential on the outer boundaries. The Default is 0

iterations = The maximum number of iterations allowed. The default is 10000.

dt = The initial step size in t. The default is 1.0e-5. The value of dt will increase by 1% for each successful iteration, for unsuccessful iterations, dt is cut in half.

nsave = all of the output is save every nsave iterations, the default is 500.

double = if set, use double precision

convergence_level = If the fractional difference objective function L is less than this then the program stops.

pr_iter = Every pr_iter iterations, an entry to the 'lhist.dat' file is made. The columns in this file are (1) the iteration, (2) the value of dt, (3) the current value of L, the objective function.

alpha_in = An initial value of alpha for use in the initial linear FFF approximation, only used for potlin = 1. The default is to use the average value of alpha calculated on the lower boundary.

slow_way = If set, then use the program calc_lin_fff to do the initial field (potential or linear FFF) instead of fff.pro. Calc_lin_fff uses a Green's function calculation (Chiu & Hilton, 1977) to get the initial field, while fff.pro uses a Fourier method. This is very slow, especially for large arrays.

dont_kill_input = If set, don't overwrite the input arrays for the case where the input is 3D. (2D input is not overwritten.

spherical = If set, then use spherical coordinates, with x, y, z becoming r, theta, phi. For this case, the input field must be 3D, and the initial conditions must be present, (bc_flag = 1 and restart_flag = 1, must be set). There is currently no method for initializing the field in spherical coordinates. We're working on it...

Added Feb-2008

use_fortran = if set, this will do the iterations using the new fortran program. Documentation for the fortran program does not yet exist.

fortran_exe = the name of the fortran routine used, the default is 'optimization_fff.exe'. Full path, please.

use_sdf = if set, all data i/o is dealt with using SDF ('simple data format') files. This includes the stuff that is passed to the FORTRAN program. These files will have .sdf as a suffix. Note that the output "lhist" files are ASCII files. To get the Software needed for use of the SDF format, go to: http://solarmuri.ssl.berkeley.edu/~fisher/public/software/SDF/ The SDF idl directory must be in the IDL path.

For weighting function options:

weight_in = an input function, for weighting the objective function, must be (rsize, tsize, vsize), i.e., the same dimensions of the 3d B field.

wfile = an input file, from which you can read in a weight function. An ASCII File needs a format

 printf, unit, rsize, tsize, vsize
 printf, unit, wf

The use_sdf functionality of this has not been tested, but those variable names are needed, rsize is the number of x (or r) grid points, tsize is the number of y (or theta) grid points, vsize is the number of z (or phi) grid points, and wf is the function name.

wtype = a number for a given functional type that will be created:

         1: Function varies as  cos((!pi/2)*(z/max(z)) Or r for spherical
         2: Function varies as cos((!pi/2)*(dist/max(dist))) Where dist
         is sqrt(x^2+y^2+z^2), not used for spherical -- if spherical
         is set then you get the same as wtype=1
         3: Function varies as the inverse of dz for variable grids
            (or dr for spherical)
         4: Function varies as the inverse of dV, from the origin 

Output:

There are no output parameters in the argument list. Instead, the magnetic field is output into a file named "field.dat". This file can be read using the program read_bfield_fff.pro, which is included in the software.

 
Read_bfield_fff, 'field.dat', xsize, ysize, zsize, bx, by, bz, x, y, z 

The file, "lhist.dat", contains the iteration history of the extrapolation, The columns in this file are (1) the iteration, (2) the value of dt, (3) the value of L, the objective function.

Other files are also created. The file "pfield.dat" contains the initial potential or linear force-free field. This does not include the lower boundary or other boundaries (for the bc_flag = 1 case). It can be read in the same way as "field.dat". Another file, "field0.dat" is created that contains the initial potential or linear FF field, but has the lower and other boundaries (for bc_flag=1) set to those that will be used in the calculation. This file is used to initialize the FORTRAN program.

If the /use_fortran keyword is set, then a file "fff_input_pars.dat" is created, for input to the FORTRAN program. This file contains:

dt                 ;the initial value of dt
convergence_level  ;the convergence_level
iterations         ;the max number of iterations
qsphere            ;1 for spherical coordinates, 0 for cartesian
slow               ; set to 0 this is not used!!!
alpha              ;the value of alpha for the init field (usually 0.0)
sdfu               ;1 if using sdf, 0 if not

If a weighting function is used, and the /use_fortran keyword is set, then a file weight0.dat is created to pass to the FORTRAN routine.

If the /use_sdf option has been set, then the field and weight files are sdf's. The lhist file, and fff_input.dat are always ascii files.

Examples:

Using optimization_fff is pretty straightforward. Given valuse for bx, by, and bz on the surface, with no non-default options, just type in:

optimization_fff, bx, by, bz

Some sample data has been included with the software in the test_data directory. The file "bfield_nlfff1.sav" contains an actual non-linear force free field created using the method of Low and Lou (1990, Apj 352, 343.). Here is how you might use the field in this file to test to see how well the extrapolation works on a known field:

restore, '/ssw/packages/nlfff/test_data/bfield_nlfff1.sav'
help, bx0, by0, bz0                       ;these will be 64x64x64 arrays
optimization_fff, bx0, by0, bz0, $        ;The input is 3D
                  bc_flag = 1, $          ;This means that all BC's are present
                  /double, $              ;Double precision works better with test data, takes longer...
                  /slow, $                ;Uses Green's function for initialization
                  convergence = 1.0e-8, $ ;this is extreme, the default is 1.0e-6
                  iter = 50000l, $        ;A high limit
                  nsave = 5000            ;save results every 5000 iterations

(This command took 816 minutes on a 2 GHz processor pc, with 2 Gbytes of RAM. But probably the calculation could have been cut off earlier, we don't have a good handle on how long to allow a test calculation to go before the answer is "good enough". Real magnetogram data, which is noisy, actually converges much more quickly.)

Since the output file name is hardwired, you might want to move the output files. (Note -- The spawn commands are for a linux/unix system. I have no clue how to do this in windows or Mac. The extrapolation code should work for any operating system, though.)

spawn, '/bin/cp lhist.dat lhist__nlfff1.dat'      ;history file
spawn, '/bin/cp field.dat field__nlfff1.dat'      ;output field
spawn, '/bin/cp pfield.dat pfield__nlfff1.dat'    ;initial (potential) field

Or you may want to save the output in a save file:

read_bfield_fff, 'field.dat', rsize, tsize, vsize, bx, by, bz, x, y, z

save, rsize, tsize, vsize, bx, by, bz, x, y, z, file = 'field__nlfff1.sav'

restore, 'bfield_nlfff1.sav'
help, bx0, by0, bz0                       ;these will be 64x64x64 arrays
optimization_fff, bx0, by0, bz0, $        ;The input is 3D, but only z=0 field is used
                  bc_flag = 0, $          ;Lower BC's are present, (default)
etc....

 fieldlines_fff, bx_in, by_in, bz_in, ix_in, iy_in, x_in, y_in, z_in, $
                 angx = angx, angz = angz, $
                 thick = thick, $
                 oplot = oplot, linecolor = linecolor, $
                 title = title, pick_lines = pick_lines, $
                 n_lines = n_lines, window_number = window_number

So say that you want to view the sample B field from the file "field__nlfff1.sav".

restore, 'field__nlfff1.sav'
fieldlines_fff, bx, by, bz, ix, iy, x, y, z, /pick_lines

You'll see the message:

'Pick      20  Fieldline starting points:'

Then click on the image 20 times, the fieldlines are drawn on the image, and the x and y starting points are saved in the arrays, ix and iy. You'll get an image:

You can change the view of the image by changing the angles:

fieldlines_fff, bx, by, bz, ix, iy, x, y, z, angx=45, angz=45

Say you want to see what the potential field for this case looked like. You didn't save it before, so you need to read in the data.

read_bfield_fff, 'pfield__nlfff1.dat', rsize, tsize, vsize, pbx, pby, pbz, x, y, z

fieldlines_fff, pbx, pby, pbz, ix, iy, x, y, z, angx=45, angz=45

Yes, the fields do look different.