Summary of Image Reconstruction Programs

Beta version (Actually works.) Useful but temporary. To be replaced. Rudimentary verion being revised.

function hsi_1d_pattern

PURPOSE:
Returns the modulation pattern cross section for a set of harmonics (e.g. 1,2,...,hmax) and grid params
(Done once for each collimator)
This is meant to replace mk_1d_pattern.pro, and is (hopefully) consistent
with Gordon Hurford's 7-element structure representing the subcollimator calibration.

CATEGORY:
Imaging

EXPLANATION:
Returns Fourier-series representation of a HESSI
subcollimator #Ncoll using coefficients in structure grm

SYNTAX:
Result=hsi_1d_pattern(Ncoll, grm1 [,harmonics] [,/positive] [/verb])

EXAMPLE:
tri= hsi_1d_pattern(8,grm_8,1+indgen(3))
The array tri is a triangle waveform (possibly distorted)

INPUTS:
Ncoll = collimator number (scalar 0-8)
grm1 = structure containing at least 7 coefficients:

grm.(0) = gridtran (scalar)
grm.(1) = modampl hmax-element vector of amplitudes
grm.(2) = modphz hmax-element vector of phases

harmonics = vector of positive integers

for example [1,3,5] or [2,4,6], Any 0 element will be removed
if it exists in the harmonics array.


 
            harmonics   grm.(1)     what you'll get
              [0,1]     whatever       [1]
              [1,2,3]   whatever       [1,2,3]
              [1,3,5]   [1,2,3]        [1,3]
              [2,4,6]   [1,2,3,4,5]    [2,4]
              [2,3,4]   [1,2,3]        [2,3]
              [2,3,4]   [1,2,3]        [2,3]
              vector1   vector2        Intersection(vector1,vector2)

   

OPT. INPUTS:
harmonics need not be specified. Default is to use all that appear in the stucture grm

KEYWORDS:
Use "/POSITIVE" to enforce positivity of the output profile.
. (The default is to set no lower bound)
Use /VERB to get explicit info on harmonics actually used
TBL_SIZE=tablesize can be used to increase the table size

COMMON :
NONE
WARNING:
Result can have negative values if slit < pitch/2 and harmonics
are few. This can have undesirable consequences in some
mapping programs. Use /positive to enforce positivity.

HISTORY:
Version 2, Based on mk_1d_patrn.pro
1998/12/02 Revised to use input GH's structure grm
Arbitrary set of harmonics (e.g. odd only)

NOTES:
Be sure pitch is in units of angle (asec?)
x is in asec (1 sec per pixel)
angpitch=? hardwired or from hessi_params? Use hessi_params for now

CONTACT:
Ed Schmahl code 682

SOURCE: hsi_1d_pattern.pro

function mk_1d_pattern

PURPOSE: Makes 1-D modulation lookup table (Done once for each collimator)

Use as input for mk_modul_patrn for speedy creation of modulation patterns

INPUTS:

Ncoll = collimator number (0-8)
s2p = slit/pitch ratio

OPTIONAL INPUTS:

harmonics (vector) -- default is [0,1,2,3]
-- may be any subset of the non-negative integers.
POSITIVE keyword: set if positivity is to be enforced
TBL_SIZE keyword: set to any value > 3000 if desired

OUPUTS:

Vector containing the cross section of a modul_patrn

CALLING SEQUENCE:

F_1d=mk_1d_pattern(8,0.55)
F_1d=mk_1d_pattern(0,0.45,tbl_size=4000)
F_1d=mk_1d_pattern(4,0.50,[0,1,3],/positive)

NOTES:

Calls hessi_params to get pitch
Phase is zero at the center of the table
Default pixel size is 1 arcsec
Default array size is 3000 (phase=0 between pixels 1499 and 1500)

function hessi_params

PURPOSE: Returns a structure with important HESSI parameters
CALLED BY: mk_1d_pattern.pro, mk_modul_patrn for speedy creation of modul_patrn arrays
INPUTS: None explicit.
OUTPUTS: A structure containing:

pitch,orient,coll_phase,slit2slat,spin_period,thickness,delta

NOTES:

Won't be needed after higher level tools exist

function mk_modul_patrn

PURPOSE: Makes modulation pattern array of size Nangles x NX x NX
INPUTS:

F1d array = cross section of pattern (from mk_1d_pattern)
theta = vector or scalar of rotation angle(s) (radians)
aspp = arc sec per pixel
Offset = [Xoffset,Yoffset} = Heliographic position of map center (arcsec)
pixels = vector of pixel x,y coordinates around the offset point

OPTIONAL INPUTS:

pixels = scalar causes a modul pattern of size pixels * pixels to be returned
Aspect =dxdy solution vectors = FLOAT(2,Nangles)
Default is to use all zeroes: dxdy = fltarr(2,Nangles)

OUTPUTS:

Example 1:

If theta is a scalar and pixels is a vector, the result is a
modulation pattern = FLOAT( n_elements(pixels)/2)

Example 2:

If both theta and pixels are vectors, the result is an array of
mod patterns = FLOAT(n_elements(theta), n_elements(pixels)/2).

Example 3:

If pixels is a 2-element vector representing a single
pixel, the result is a vector = FLOAT(n_elements(theta)).

CALLING SEQUENCE:

For a 64 x 64 map one may use:
xyarray=[(findgen(64^2) mod 64)-64/2+.5, fix(findgen(64^2)/64)-64/2+.5]
mod_pat=mk_modul_patrn(F1D,theta,aspp,offset,xyarray)

The following returns the same result:

mod_pat=mk_modul_patrn(F1D,theta,aspp,offset,64)

The modulation pattern at a single point [3,-5] is returned by:

mod_pat=mk_modul_patrn(F1D,theta,aspp,offset,[3,-5])

CALLED BY: model_to_score, score_to_bproj, max_entropy, ..
NOTES:

The patterns produced here are arrays containing a "flattened" (reformed) version of the map. This is useful for applications which perform matrix multiplications with the modulation patterns (back-projections, count-rate models, etc.)
No interpolation of the simulated aspect solution is done in the present version.

function model_to_score_params

PURPOSE: Returns a structure with the following parameters needed by model_to_score.pro :
INPUTS: Selected by user
OUTPUTS: parameter structure:

coll=list of collimators (vector), e.g. [0 8]
aspp = arc sec per pixel in model map, e.g., 2.0
offset = vector offset of map center from spin axis (arcsec)
stdout = -1 for screen output, N>0 for output to disk,
dev ='x', 'ps', 'null'
model_name = full path name of model file
time_range = Start and end time, e.g. [0,4]
factor = Multiplicative factor for map, e.g. 1000.
chan_range ; Channel range (2 integers from 0 to 8191)

CALLING SEQUENCE: params=model_to_score_params()
NOTES:

Will be replaced by RS&KT upper-level tasks

function bproj_params

PURPOSE: Returns a structure with the parameters needed by score_bproj.pro:
INPUTS: Selected by user
OUTPUTS: parameter structure:

score_name = full path name of score file
coll=list of collimators (vector), e.g. [0 8]
aspp = arc sec per pixel in model map, e.g., 2.0
offset = vector offset of map center from spin axis (arcsec)
stdout = -1 for screen output, N>0 for output to disk,
dev ='x', 'ps', 'null'
time_range = Start and end time, e.g. [0,4]

CALLING SEQUENCE: params=bproj_params(0)
NOTES:

Awaiting glorious, full-featured RS&KT version.

pro model_to_score

PURPOSE: Makes a score from a given model flare map set
INPUTS:

input_struct: The structure returned by model_to_score_params containing the following:

coll=list of collimators (vector), e.g. [0,8]
aspp = arc sec per pixel in model map, e.g., 2.0
offset = vector offset of map center from spin axis (asec)
stdout = -1 for screen output, N>0 for output to disk,
dev ='x', 'ps', 'null'
model_name = full path name of model file
time_range = Start and end time, e.g. [0,4]
factor = Multiplicative factor for map, e.g. 1000.
chan_range ; Channel range (2 integers from 0 to 8191)


OPTIONAL INPUTS:
SASZERO: If set, the simulated aspect solution will be zeroed out.

OUTPUTS:

A FITS file containing:

Time tags, collimators, pulse heights
ID of source (in FITS header)
ID of aspect (in FITS header)
Name contains a unique Magic Number

An IDL SAVE file containing:

Simulated aspect solution with t, dx,dy array,
Aspect save file has same unique Magic Number as score FITS file

CALLING SEQUENCE:
Example 1: model_to_score,input_structure
Example 2: model_to_score,input_structure /SASZERO

CALLED BY: User
MODIFICATION HISTORY

Fixed a bug in an internal function which modified the pulse-height distribution;
Enabled use of min and max pulse-height channels. 9/23/98

pro score_bproj

PURPOSE: Select score file and make a back-projection image.
EXPLANATION:
The score is a FITS file containing photon time tags, detector #s,
and amplitudes as a 3 X N long integer array.
It was produced by the program model_to_score.pro from a flare model
Time tags (score(0,*)) are long integers in units of micro seconds
Detector numbers (score(1,*)) run from 1 to 9
Amplitudes (score(2,*)) are (currently) random integers in the range 0-16383,
In later versions of score_bproj.pro, the amplitudes will be synthetic pulse heights.

INPUTS:

bproj_params_structure
(See bproj_params.pro)

OPTIONAL INPUTS:

If ASPECT is set, [dx,dy] is taken from the savefile created by model_to_score
The Default is to use [dx,dy]=0 vector

METHOD:
The variable "score" is retrieved from the fits file, and
the nine binned modulation profiles are displayed on device dev.
A full-sun backprojection is then performed using the specified
detectors by adding on-the-fly modul_patrns together.

CALLED BY: User

pro score_binner

PURPOSE: To bin the time tags of a score
INPUTS:

Select score file
Choose time range

OUTPUTS: A structure containing:

Histograms of timetags=count-rate profiles
ID of original score
Parameters of the binning

NOTES:

Default bins: 3645, 2104, 1215, 701, 405, 233, 135, 77, 44
Default is 9 collimators

function likelihood

PURPOSE:
To determine goodness of fit of f_model to f_obs (Similar to chi-square statistic)

INPUTS:
f_model = current model count rate vector
f_obs = observed count rate vector (all integers or near-integers)
modul_patrn = array of modulation patterns

OPTIONAL INPUTS:
KEYWORD epsilon can be set to replace zeroes of f_obs with a small value.
This is to help with reconstruction programs that fail ehen f_obs has zeroes.

OUTPUTS:
Output is vector of log-likelihood terms

CALLED BY:
Max Entropy program
...

CALLING SEQUENCE:
Result = likelihood(f_fmodel,f_obs [epsilon=tiny_value])
chi_squared=total(result) [approximately]

NOTES:

Total is asymptotic to chi-squared statistic
See W. Cash (Ap. J. 1979) for details
Has been revised to provide use of zero count rates
Warnings are given if f_obs has non-integer values.
Alog(factorial(x)) has been replaced by lngamma(x+1)

function grad_likelihood

PURPOSE: Computes gradient of the total likelihood
INPUTS:

f_model
f_obs
modul_patrn

OUTPUTS:

Map of log-likelihood gradient
Actually looks better than back-projected map

CALLED BY:

Max Entropy to make the iterated map

NOTES:

Needs to be revised to provide use of zero count rates
Used by :
Fixed-point (Sakao/Willingale) MEM
conjugate gradient technique (in prep)

Can be used to make a fast dirty map.

function aspect_sim

PURPOSE:

Quick and dirty aspect simulator
Computes arrays of simulated displacements of the telescope axis from sun center.
Creates high-freq jitter of small amplitude, circle of large radius

INPUTS:

times = A vector of times (in sec)
z is optional: if z=0, dx,dy will be all zeroes

OUTPUTS:

[dx,dy} as a function (arcsec) of vector times (sec)

NOTES:

Presently outputs a deterministic function for debugging reasons
Next version will produce pseudo-random jitter.

pro Sakao-MEM

PURPOSE: Maximizes entropy using Sakao/Willingale fixed-point method
INPUTS:

Observed count rate
Modulation patterns


NOTES:

On hold while testing problems caused by zeroes in count rate.

pro conj_grad_mem

PURPOSE: Maximizes entropy using conjugate gradient method
INPUTS:

Same as for Sakao-MEM

NOTES:

Presently vaporware
Working test version not converted to HESSI use (yet)

phi2score:

PURPOSE:

Converts a 2-D array of angles phi to a sorted time-tagged lists with first column time and 2nd column detector number.
Called repeatedly to add data computed from model pixels

INPUTS:

phi = float(Nphi,Ncoll)
-- 1st argument = angle, 2nd argument = collimator number
pulseht = float(Nphi,Ncoll) -- pulse height associated with phi
spin_period (seconds) -- assumes time and angle are exactly proportional

RESTRICTIONS:

assumes there is at least 1 non-negative angle for each collimator
assumes there are at least 2 collimators
collimator numbers start with 0, detector numbers start with 1
pulse ht array is repeated for each pixel

OUTPUT:

array=float(3,N)
time tag(usec) detector pulse height (random)

function time_order,score

PURPOSE: Puts a score into time order
INPUTS: Array score(0,*)
OUTPUTS:time-sorted score array

Last updated September 23, 1998.