Version notes: v2.1.1 [2021.07.15]: Generic FITS + Hinode SP L2 input, fixed bug in spherical geometry v2.1 [2019.11.05]: Energy calculations performed in double precision v2.0 [2018.01.12]: Rewritten in Fortran 90 v1.2 [2016.12.24]: Update to dfti_example_status_print.f90 to account for update to MKL library This code package for automated ambiguity resolution of vector magnetic field data is loosely based on the Minimum Energy ambiguity resolution method described by Metcalf (1994, Solar Phys., 155, 235). Some description of the present version is described in Leka, Barnes & Crouch 2009, Proceedings of the 2nd Hinode Science Meeting. Please email graham@nwra.com with questions about the code package. In this implementation, simulated annealing to find the global minimum to the function E = |div B| + |J_z| where the vertical derivative in the divergence term is approximated by the vertical derivative of the potential field with the observed line of sight component on the boundary. After the ambiguity has been resolved based on the minimum energy state found, pixels below a given threshold in the transverse field strength are revisited using an iterative acute angle to nearest neighbors method, similar to that described in Canfield et al. (1993, Astrophys. J., 411, 362). This is done because it has been found by Leka et al. (2009, Solar Phys.) that the minimum energy state is not the correct ambiguity resolution in the presence of noise. Thus, in areas where the uncertainty in the azimuth is larger, the method looks for a smooth solution. Details of the package, how to compile and run, and an example, are given below. ***** The code package consists of the following files: ambig.f90 anneal.f90 bobs.f90 bounds.f90 boxit.f90 buffer.f90 CalcDE.f90 CalcDE_reconfig_OCBP_PF_dzh_4p.f90 CalcDE_reconfig_spherical_PF_4p.f90 CalcE_OCBP_PF_dzh_4p.f90 CalcE_spherical_PF_4p.f90 constant.f90 cuspl2d.f90 cuspleval2d.f90 cuspleval.f90 cuspl.f90 deriv_coefficients.f90 dfti_example_status_print.f90 disk_center.f90 energy_arrays.f90 ephemeris.f90 fits.f90 flags.f90 get_ran_pix.f90 global.f90 grow1.f90 maskvec.f90 mgram_data.f90 minimise_energy.f90 mollweide.f90 nacute5.f90 nacute6.f90 ntiles.f90 pacute.f90 pad.f90 pix_size.f90 point.f90 potential.f90 pot_field.f90 printerror.f90 ran3.f90 ran_pix.f90 ranseed.f90 rd_field.f90 rd_fits.f90 rd_param.f90 rd_sotsp.f90 recon.f90 reconfig.f90 revmollweide.f90 setup_mask.f90 setup_OCBP_PF_dzh_4p.f90 setup_spherical_PF_4p.f90 sizes.f90 smooth.f90 sortrx.f90 spherical_deriv_coefficients.f90 spherical_position.f90 tile.f90 transform.f90 trnsfrm.f90 update_fits.f90 upmap.f90 verbose.f90 WeightingFactor.f90 write_field.f90 par makefile ***** By default, it requires the Intel Math Kernel library (MKL) for performing FFTs, and the FITSIO library for I/O. The paths containing these files must be correctly set in the make file. The par file contains the details of the case being running in the following form: filename irflag ibflag iaflag igflag ipflag npad nap ntx nty athresh bthresh nerode ngrow iaunit ipunit incflag iseed iverb neq lambda tfac0 tfactr filename is the name of the file containing the input data. irflag determines whether the input data are in generic FITS format (irflag=1), in the FITS format with extensions found in the Hinode SOT/SP pipeline (irflag=2), or in a formatted text file (irflag=0). In the Hinode (irflag=2) case, output will be to a (L2.1) FITS file, with no azimuth file written. ibflag determines whether the field is specified as line of sight and transverse components plus the azimuthal angle (ibflag=0) or as the magnitude of the field, plus the inclination and azimuthal angles (ibflag=1). iaflag determines the direction of zero azimuthal angle: CCD+x (iaflag=0), CCD+y (iaflag=1), CCD-x (iaflag=2), CCD-y (iaflag=3). igflag determines whether to use planar (igflag=1) or spherical geometry (igflag=2). For a given number of pixels, the planar version is approximately a factor of two faster, but if the area in the field of view is a significant fraction of the solar disk, spherical geometry should be used. ipflag determines whether to perform a potential field acute angle ambiguity resolution (ipflag=1) or not (ipflag=0). npad determines the number of pixels to add on each side of the field of view, to mitigate the effects of the period boundary conditions in using FFTs to calculate the potential field. For full disk images, it is recommended that npad be set to approximately 10% of the number of pixels across the disk. nap determines the number of pixels over which to smoothly transition the field to zero outside the field of view (planar geometry), or number of pixels outside a tile to include in the potential field calculation (spherical geometry). It must be less than npad else will be reset to equal npad. ntx sets the number of tiles in the x-direction (longitude) used to construct the potential field derivatives when the code is run using spherical geometry. nty sets the number of tiles in the y-direction (latitude) used to construct the potential field derivatives when the code is run using spherical geometry. athresh is the transverse field strength above which the annealing will be used. Typically a value of 200G < athresh < 400G produces reasonable results, but athresh should be set lower than bthresh. It may be useful to start with a value about the 6-sigma detection level in the transverse field. However, if speed is not a concern, then setting athresh=0 may produce the best results. bthresh is the transverse field strength below which the nearest neighbor acute angle method will be uses. Typically a value of 200G < bthresh < 400G produces reasonable results, but bthresh should be set higher than athresh. It may be useful to start with a value about the 6-sigma detection level in the transverse field. nerode is the number of pixels by which to erode the mask, to prevent the inclusion of isolated above-threshold pixels in the annealing. nerode=1 is generally sufficient. ngrow is the number of pixels by which to subsequently grow the mask, so allow a buffer around above-threshold pixels in the annealing. ngrow=1 may be all that is necessary. iaunit specifies whether azimuth is in radians (iaunit=0) or degrees (iaunit=1). ipunit specifies whether pointing is in radians (ipunit=0) or degrees (ipunit=1). incflag specifies whether inclination of 0 is vertical (incflag=0) or horizontal (incflag=1). iseed is used for initializing the random number generator. Must be a positive integer. iverb controls additional output: iverb=0 supresses all print statements; iverb=1 will print a statement indicating when each stage of the code finishes; iverb=2 will print the energy at each temperature in the annealing, plus number of pixels flipped during each iteration in smoothing. The latter option produces a large amount of output, but can be useful to see how well the solution is converging. neq determines the number of times pixels are visited at a given temperature. Must be a positive integer. lambda is the weighting factor between the divergence term and the vertical current density term. Must be positive (or zero). tfac0 scales the initial temperature. Must be positive (and typically should not be smaller than of order 1). tfactr determines the cooling rate. Must have 0 < tfactr < 1, with slower cooling for values closer to 1. Slower cooling is more likely to find the global minimum, but takes longer to run. If the input data are in generic FITS format, the following information must also be present in the par file: iblong ibtrans ibazim longitude keyword latitude keyword x-center keyword y-center keyword p-angle keyword b-angle keyword radius keyword pixel size in x-direction keyword pixel size in y-direction keyword iblong, ibtrans, ibazim specify the array element in the FITS file containing the line of sight component, transverse component and azimuth or magnitude of field, inclination and azimuth. longitude keyword specifies the keyword name (max 8 character string) in the FITS file containing the longitude. If keyword is not found, the code checks for xcen and ycen keywords instead. latitude keyword specifies the keyword name (max 8 character string) in the FITS file containing the latitude. If keyword is not found, the code checks for xcen and ycen keywords instead. x-center keyword specifies the keyword name (max 8 character string) in the FITS file containing the x-position of the field of view relative to disk center. If the keyword is not found (and the longitude/latitude were not found), then the code assumes the image is at disk center. y-center keyword specifies the keyword name (max 8 character string) in the FITS file containing the y-position of the field of view relative to disk center. If the keyword is not found (and the longitude/latitude were not found), then the code assumes the image is at disk center. p-angle keyword specifies the keyword name (max 8 character string) in the FITS file containing the p-angle. If keyword is not found, a warning is printed and a value of p-angle=0. is used. b-angle keyword specifies the keyword name (max 8 character string) in the FITS file containing the b-angle. If keyword is not found, a warning is printed and a value of b-angle=0. is used. radius keyword specifies the keyword name (max 8 character string) in the FITS file containing the solar radius. Must be in same units as pixel size. If keyword is not found, a warning is printed and a value of radius=959. is used. pixel size in x-direction keyword specifies the keyword name (max 8 character string) in the FITS file containing the pixel size in the x-direction. Can be in any units, provided x- and y-pixel sizes are in the same units. pixel size in y-direction keyword specifies the keyword name (max 8 character string) in the FITS file containing the pixel size in the y-direction. Can be in any units, provided x- and y-pixel sizes are in the same units. If the input data are from the Hinode SP pipeline, such that all the images are contained in extensions, the following information must be present in the par file (in place of the information above for a generic fits file): jfs jfi jfa jsf jci jfs is the HDU number containing the field strength (typically jfs=2, for the first extension). jfi is the HDU number containing the field inclination (typically jfi=3, for the second extension). jfa is the HDU number containing the field azimuth (typically jfa=4, for the third extension). jfs is the HDU number containing the straylight fraction (typically jsf=13, for the twelfth extension). jci is the HDU number containing the continuum intensity (typically jci=33, for the thirty-second extension). If the input data are not in FITS format, the input data must be in the following form, in the file specified by "filename" in the par file: nx ny xpix ypix b p radius theta phi bl(nx,ny) bt(nx,ny) ba(nx,ny) nx, ny are the dimension of the input arrays. xpix, ypix are the pixel sizes in the x- and y-directions. b is the solar b-angle. p is the solar p-angle. radius is the solar radius. theta is the central meridian angle of the center of the field of view. phi is the latitude of the center of the field of view. bl is an nx by ny array of the line of sight component of the field bt is an nx by ny array of the transverse component of the field ba is an nx by ny array of the azimuthal angle, containing the ambiguity. If the input data are in generic FITS format, then the azimuthal angle in the FITS file will be updated to reflect the ambiguity resolution, and a comment will be added to the file, indicating that the ambiguity was resolved using this method. If the input data are not in FITS format, a new file, named "azimuth.dat" will be created, containing a formatted array of the ambiguity-resolved azimuth angles, of the same size as the original. In either case, the angles will be given in the same units (degrees/radians) as the input angles. ***** To run the code: 1) Familiarize yourself with the options listed above regarding file type, components, and angle definitions for your data. 2) compile the code: -edit the file "makefile" for locations of the MKL libraries, and other specifics. -on the command line: %> make This should produce an executable called "ambig". 3) In the directory where you will run the code, copy the file "par" and edit to reflect your data. 4*) on the command line: %> ambig (it looks for the file "par"; see the previous section for additional command-line usage options and examples). The result will be stored as described above depending on your data type. *See the section at the very end for alternative calling options using command line arguments, and Hinode-centric examples. ***** EXAMPLE CALL AND OUTPUT: leka@kolea%> ls ambig.tgz leka@kolea%> tar -xvf ambig.tgz ambig.f90 anneal.f90 bobs.f90 bounds.f90 boxit.f90 buffer.f90 CalcDE.f90 CalcDE_reconfig_OCBP_PF_dzh_4p.f90 CalcDE_reconfig_spherical_PF_4p.f90 CalcE_OCBP_PF_dzh_4p.f90 CalcE_spherical_PF_4p.f90 constant.f90 cuspl2d.f90 cuspleval2d.f90 cuspleval.f90 cuspl.f90 deriv_coefficients.f90 dfti_example_status_print.f90 disk_center.f90 energy_arrays.f90 ephemeris.f90 fits.f90 flags.f90 get_ran_pix.f90 global.f90 grow1.f90 maskvec.f90 mgram_data.f90 minimise_energy.f90 mollweide.f90 nacute5.f90 nacute6.f90 ntiles.f90 pacute.f90 pad.f90 pix_size.f90 point.f90 potential.f90 pot_field.f90 printerror.f90 ran3.f90 ran_pix.f90 ranseed.f90 rd_field.f90 rd_fits.f90 rd_param.f90 rd_sotsp.f90 recon.f90 reconfig.f90 revmollweide.f90 setup_mask.f90 setup_OCBP_PF_dzh_4p.f90 setup_spherical_PF_4p.f90 sizes.f90 smooth.f90 sortrx.f90 spherical_deriv_coefficients.f90 spherical_position.f90 tile.f90 transform.f90 trnsfrm.f90 update_fits.f90 upmap.f90 verbose.f90 WeightingFactor.f90 wr_sotsp.f90 write_field.f90 par makefile leka@kolea%> make ifort -g -c /opt/local/mkl/include/mkl_dfti.f90 ifort -g -c dfti_example_status_print.f90 ifort -extend_source -w ....... The output for running the code will depend on what options are set, particularly whether the code is using spherical or planar geometry. leka@kolea%> ./ambig parameters read field read transform calculated buffering done apodizing done potential field done optimization done smoothing done results written ***END*** If the file is a generic FITS (irflag=1), the FITS file will be updated. If a Hinode FITS or formatted array file, results written to a formatted output file named "azimuth.dat', of the same size as the original arrays. Call the size nx, ny. To read (for example in IDL): IDL> azim=fltarr(nx,ny) IDL> openr,1,'azimuth.dat' IDL> readf,1,azim IDL> close,1 ***** Alternative calling sequence using command line arguments These options were developed for use with Hinode L2.1 fits files, but are useful for general scripting, where par file, input/output filenames, and random number seeds may be specified on the command line. By default, the code is run without any command line arguments, and reads a parameter file named 'par', where the input filename is specified on the first line as described. Command line arguments supercede defaults and/or values specified within the par file. These are strictly positional, so they must be given in the following order, with a placeholder for azimuth output file filled in as necessary. 1. parameters filename [default is 'par', if not given] 2. azimuth output file [default is 'azimuth.dat' if not given] 3. Input file [overrides what is specified in the par file] 4. Hinode output file [output FITS filename for the case of irflag=2] 5. Random number seed [overrides what is what is given in the par file] Two examples follow, using the included 'par_hinode' parameters file and a 20151122_185049.fits FITS file. Note that 'junk' is specified for the azimuth file as a placeholder, as no azimuth file is written in this case. 1. Run with input from 20151122_185049.fits, and (default) output to be written to 20151122_185049_L2.1.fits: ./ambig par_hinode junk 20151122_185049.fits 2. Input from 20151122_185049.fits, and run with a random number seed of 3, with output to be written to 20151122_185049_L2.1_s3.fits (specified): ./ambig par_hinode junk 20 20151122_185049.fits 20151122_185049_L2.1_s3.fits 3