;+ ; Returns an image of the SXT point-spread-function. ; FUNCTION sxt_psf, index, xc, yc, res=res, dim=dim, $ spower=spower, fourier=fourier, $ moffat_only=moffat_only, delta=delta, log10t=log10t, $ no_shadow=no_shadow, scat_fraction=scat_fraction, $ noverbose=noverbose ; ; REQUIRED INPUT PARAMETERS: ; index = SXT image index, used to get image center, filter, ; and resolution. ; OR ; filter = number of the X-ray analysis filter, GT-FILTB convention. ; In this case you must specify res. ; xc, yc = center of FOV in CCD Full Resolution pixels. ; res = The resolution [0,1,2] of the image that this PSF will ; be compared to. ; ; OUTPUT: ; psf = image of the psf. Default is: ; odd dimension (65x65), array centered, ; normalized to unit integral. ; ; OPTIONAL INPUT KEYWORDS: ; dim = size of output PSF image in [res]-resolution pixels. ; Default = 65. ; ; fourier = Return psf ready for use by IDL procedure FFT. ; Psf image is dimensioned power of 2 (64 x 64), ; corner centered, normalized to unit mean. ; ; delta = returns scattering wings plus a delta at the center ; only, Moffat not included. (Default is Moffat core ; plus scattering wings) ; ; moffat_only = returns Moffat psf only. Scattering wings not ; included. ; ; no_shadow = The default mode is to use a scattering wing model ; that includes the "ray" shadows of the filter supports. ; This is currently our best fit to the scattering wing, ; with a slope of -2 in the radial dependence. ; Selecting NO_SHADOW will defeat this, and assume a ; cylindrically symmetric, smooth scattering wing with ; no "rays". ; ; spower = Slope of the scattering wing, with or without "rays", ; is set by SPOWER. Default SPOWER is -2. ; ; scat_fraction = fraction of total signal to be attributed to the ; scattering wings. The default values are based on ; McKenzie's analysis of May 2002; but with this optional ; parameter you can force X percent of the signal into ; the wings; e.g., scat_fraction=0.073 for 7.3% scattering ; ; log10t = Log10 of temperature. The Moffat parameters ; are wavelength, and therefore, temperature ; dependent. Default uses values derived ; assuming equal emission measures at all ; temperatures 5.5 <= log10t <= 8.0. ; Also applies wavelength dependency to scattering ; wings, via (ref. wavelength/log10T wavelength)^2 ; ; CALLING EXAMPLES: ; psf = SXT_PSF(index) ; psf = SXT_PSF(filter, xc, yc) ; psf = SXT_PSF(filter, xc, yc, /delta, scat_fraction=0.152) ; ; METHOD: ; Uses the Moffat function fits described by Piet Martens. ; Adds scattering wings described by Acton. ; Wavelength weighting done by SXT_WEIGHT procedure. (Not so) ; Creates arrays with equivalent field of view in full-res, ; and calculates the PSF in full-res. Then resizes the PSF ; array to the desired output size and resolution. ; ; RESTRICTION: ; Scattering has been empirically estimated for all filters, ; using full-res unsaturated flare images. Some uncertainty ; still exists, but the defaults are the best current estimates. ; (May 2002, DMcK) ; ; PROGRAMS CALLED: ; SIZE, GT_FILTB, GT_CENTER, GT_RES, FINDGEN, REBIN, TRANSPOSE, ; ATAN, SXT_MWAVE, DELVARX, SXT_SCATWINGS, TOTAL, STDEV, SHIFT ; ; HISTORY: ; Written April 4, 1994 Barry LaBonte ; Added /DIM keyword May 23, 1994 BJL ; Added /FOURIER, /MOFFAT keywords, scattering wings. June 27, 1994 BJL ; Modified to use Moffat parameters weighted ; by SXT passband - Mewe spectrum product, ; added /LOG10T keyword. March 24, 1995 BJL ;V2.00 08-Feb-2001, DMcK - Incorporated many options that ; were added to hacked versions of the last two years, or ; which grew out of deconvolution workshop: resolution ; flexibility; optional SLEVEL, SPOWER; log10T wavelength ; dependency in the scattering wings; delta function replacement ; for Moffat core; SCAT_FRACTION option. Updated version to ; 2.00. ;V2.10 09-Apr-2001, DMcK - Incorporated Acton's nifty program ; SXT_SCATWINGS to get a scattering wing model that includes ; the "ray" shadows of the filter supports. Also removed ; one call to REBIN that was disrupting the symmetry of the ; "ray-less" scattering wing. Made SXT_SCATWINGS the ; default scattering model, with introduction of /no_shadow ; optional switch. Updated version to 2.10. ;V3.00 23-May-2002, DMcK - Major restructuring. ; Changed /MOFFAT to /MOFFAT_ONLY for clarity; added ; documentation; removed redundant SLEVEL input; ; Moffat core PSF is default; set scattering level ; hierarchy to SCAT_FRACTION > LOG10T > defaults. ; Default SCAT_FRACTIONs based on analysis of 35 over-the-limb ; flares throughout the mission, and analysis of spikey arcade ; from 18-dec-98. Promoted to Version 3.00. ; ;- progverno = 3.000*1000 if (keyword_set(spower) eq 0) then spower = -2.00 ; April 2001, Acton finds that a slope of -2. is the best fit. if not(keyword_set(noverbose)) then begin print,' scattering index =',spower endif ; Handle inputs sz = SIZE(index) type = sz(N_ELEMENTS(sz)-2) IF( type eq 8 ) THEN BEGIN filter = GT_FILTB(index(0)) xy = GT_CENTER(index(0)) xc = xy(0) yc = xy(1) ; Find the image resolution if n_elements(res) LE 0 then begin res = gt_res(index) endif else res=res if not(keyword_set(noverbose)) then print,' Resolution = ',gt_res(res,/st) ENDIF ELSE BEGIN filter = index ENDELSE ; How big to make output? orsize = 65 IF( KEYWORD_SET(fourier) NE 0 ) THEN orsize = 64 IF( KEYWORD_SET(dim) NE 0 ) THEN orsize = dim orhalf = orsize/2 ; Coordinate images FRE_size = orsize * 2.^(res) ; Equivalent FOV in full-res FRE_half = FRE_size / 2 x = findgen(FRE_size) - float(FRE_half) x = rebin(x, FRE_size, FRE_size) y = transpose(x) if not(keyword_set(delta)) then begin ; MOFFAT Least-squares Fit parameters for equal emission measure ; Filter order is given by GT_FILTB(): ; 1 Open Op ; 2 Al.1 A1 ; 3 AlMg AM ; 4 Be119 Be ; 5 Al12 A2 ; 6 Mg3 Mg aeqem = [0.77, 0.76, 0.76, 0.96, 0.71, 0.75] beqem = [1.53, 1.52, 1.52, 1.66, 1.48, 1.51] xceqem = [529.0, 528.8, 528.7, 533.0, 527.8, 528.7] yceqem = [585.7, 584.9, 583.7, 586.3, 579.6, 583.8] eps0 = [-5.54e-4, -3.18e-3, -5.05e-3, 8.00e-2, -2.47e-2, -6.82e-3] eps1 = [1.40, 1.41, 1.42, 1.12, 1.49, 1.42] eps2 = [-3.66, -3.69, -3.71, -3.32, -3.86, -3.72] eps3 = [2.69, 2.70, 2.73, 2.72, 2.81, 2.72] ; Temperature correction coefficients atc = [ [ 0.223, 0.230, 0.123, -0.208, 0.174, 0.072], $ [-0.447, -0.475, -0.235, 0.118, -0.465, -0.162], $ [ 0.180, 0.194, 0.096, 0.005, 0.208, 0.074] ] btc = [ [ 0.218, 0.224, 0.124, -0.155, 0.170, 0.074], $ [-0.426, -0.453, -0.225, 0.084, -0.443, -0.155], $ [ 0.168, 0.182, 0.088, 0.007, 0.196, 0.066] ] xctc = [ [ 5.39, 5.53, 3.01, -4.41, 4.19, 1.79], $ [-10.65, -11.32, -5.61, 2.45, -11.08, -3.86], $ [ 4.24, 4.59, 2.25, 0.15, 4.93, 1.71] ] yctc = [[ 31.93, 32.78, 19.59, -5.25, 24.88, 12.43], $ [-59.11, -62.77, -31.45, 0.60, -60.86, -21.45], $ [ 22.06, 23.93, 10.81, 1.50, 26.11, 7.68]] e0tc = [[0.108, 0.111, 0.060, -0.088, 0.084, 0.036], $ [-0.213, -0.226, -0.112, 0.049, -0.222, -0.077], $ [0.085, 0.092, 0.045, 0.003, 0.099, 0.034]] e1tc = [[-0.371, -0.381, -0.207, 0.309, -0.289, -0.123], $ [0.734, 0.781, 0.387, -0.172, 0.765, 0.266], $ [-0.293, -0.317, -0.156, -0.010, -0.341, -0.118]] e2tc = [[0.934, 0.959, 0.548, -0.452, 0.727, 0.338], $ [-1.786, -1.898, -0.946, 0.216, -1.849, -0.648], $ [0.689, 0.746, 0.352, 0.035, 0.808, 0.259]] e3tc = [[-0.667, -0.685, -0.413, 0.066, -0.520, -0.264], $ [1.226, 1.302, 0.653, 0.017, 1.261, 0.445], $ [-0.454, -0.493, -0.221, -0.033, -0.539, -0.155]] etc = [ [[e0tc]], [[e1tc]], [[e2tc]], [[e3tc]] ] ; Which fit to use? idex = filter - 1 ; Get parameter values afit = aeqem(idex) bfit = beqem(idex) xcent = xceqem(idex) ycent = yceqem(idex) eps = [eps0(idex), eps1(idex), eps2(idex), eps3(idex)] ; Corrections to equal emission measure values as quadratic ; in Log10(temperature) IF( KEYWORD_SET(log10t) ) THEN BEGIN dt = log10t - 6. dt2 = dt^2 afit = afit + atc(idex,0) + atc(idex,1)*dt + atc(idex,2)*dt2 bfit = bfit + btc(idex,0) + btc(idex,1)*dt + btc(idex,2)*dt2 xcent = xcent + xctc(idex,0) + xctc(idex,1)*dt + xctc(idex,2)*dt2 ycent = ycent + yctc(idex,0) + yctc(idex,1)*dt + yctc(idex,2)*dt2 eps = eps + etc(idex,0,*) + etc(idex,1,*)*dt + etc(idex,2,*)*dt2 ENDIF ; Polar angle phi = ATAN( yc - ycent, xc - xcent ) ; Distance d = SQRT( ( xc - xcent )^2 + ( yc - ycent )^2 ) / 500. ; Ellipticity epsilon = eps(0) + eps(1) * d + eps(2) * d^2 $ + eps(3) * d^3 ; "Radius" r2 = (x^2 * (1. + epsilon *(1. + COS(2.*phi))) $ + y^2 * (1. + epsilon*(1. - COS(2.*phi))) $ + x * y * 2. * epsilon * SIN(2.*phi)) ; Moffat core of the PSF moffpsf = 1. / (1. + r2/afit^2 )^bfit moffpsf = rebin(moffpsf,orsize,orsize) * 4.^(res) psf = moffpsf endif ; End Moffat PSF section ; Add scattering wings if needed IF not(KEYWORD_SET(moffat_only)) THEN BEGIN ; If not Moffat-only, then ; include scattering wing if not(keyword_set(scat_fraction)) then begin case 1 of ; DEFAULT VALUES filter eq 1: scat_fraction=1.e-8 ; Not used, I hope. filter eq 2: scat_fraction=0.035 ; AL.1 filter eq 3: scat_fraction=0.035 ; ALMG filter eq 6: scat_fraction=0.050 ; Mg, assumption filter eq 5: scat_fraction=0.200 ; AL12 filter eq 4: scat_fraction=0.200 ; Be endcase ; No evidence was found for time-dependent ; scat_fraction...DMcK, 16-May-02 slevel = scat_fraction / 16. ; Empirical guesstimate. ; For Moffat core. For delta core, ; this will get overridden by SCAT_FRACTION below. if keyword_set(log10t) then begin ; Scattering wings parameters vs increasing energy reflam=sxt_mwave(6.84136,index) ; However, this requires an ; index structure to be input. ; 6.84136 is from the 6.94MK reference flare data lamrat=sxt_mwave(log10t,index) ; lamrat=reflam/lamrat ; slevel=slevel*(lamrat)^2 ; Swiped from SXT_SCATIM endif endif else begin scat_fraction=scat_fraction slevel = scat_fraction / 16. ; Empirical guesstimate, see above. endelse ; Level of scattering wing has now been set. ; Insert ray-less scattering wing if SXT_SCATWINGS is not to be used. if keyword_set(no_shadow) then begin r2=(x^2 + y^2) * 0.25 r2(FRE_half, FRE_half) = 1000. ; suppress the center point subset_low = fre_half-orhalf subset_hi = fre_half+orhalf-1 delvarx, x, y, FRE_size, FRE_half scat = slevel * r2^(spower/2.) delvarx, r2 scat = scat(subset_low:subset_hi,subset_low:subset_hi) endif else begin ; Insert SXT_SCATWINGS to include the rays. (DEFAULT) scat=sxt_scatwings(res,slope=spower) * (5.*slevel) / 3.e6 ; The above (5*slevel/3e6) is empirically designed to make ; the Moffat+sxt_scatwings have approximately the right amount ; of scattering. For Delta+sxt_scatwings, that factor will be ; overridden below. scatsz=size(scat) scatsx=scatsz(1) subset_low=scatsx/2 - orhalf subset_hi=scatsx/2 + orhalf -1 scat=scat(subset_low:subset_hi,subset_low:subset_hi) ; End SXT_SCATWINGS insertion. endelse if keyword_set(delta) then begin ; Replace PSF with delta+wings psf = scat psf[orhalf,orhalf] = 0. ctrlmt=64./(2^res) scat_center = psf(orhalf-ctrlmt:orhalf+ctrlmt-1,orhalf-ctrlmt:orhalf+ctrlmt-1) scat_total = total(scat_center,/double) psf = (psf / scat_total) * scat_fraction psf[orhalf,orhalf] = 1. - scat_fraction ; I'm using double here because we're adding small nos. ; I set the middle 128x128 of a full-res to sum up to ; scat_fraction, because the scat_fractions were measured ; with 64x64 F-res flare PFIs. If dim > 128, then the ; outer parts of the wings will contribute a bit to the ; total PSF, like they should, because of the larger fov. ; If half- or qrtr-res, then middle 64x64 or 32x32 portion will ; be total up to scat_fraction. endif else psf = moffpsf > scat ENDIF ; What style of output? IF not(KEYWORD_SET(fourier)) THEN BEGIN ; Normalize to unit total psf = float(psf / TOTAL(psf,/double)) ENDIF ELSE BEGIN ; Normalize to unit average pjunk = STDEV(psf, average) psf = psf / average ; Shift to corner centered psf = SHIFT(psf, -orhalf, -orhalf) ENDELSE RETURN, psf END