PRO SXT_TEEM1, index1, image1, index2, image2, te, em, d_te, d_em, Valid, $ t1=tt1, t2=tt2, interp=interp, gain=gain, $ ; Input thresh1=thresh1, thresh2=thresh2, date=date, $ ; Input average=average, sum=sum, satval=satval, $ ; Input sat1=sat1, sat2=sat2, unc1=unc1, unc2=unc2, $ ; Input n_params0=n_params1, version=version, $ ; Input mewe=mewe, corona=corona, hybrid=hybrid, $ ; Input photospheric=photospheric, sre_file=sre_file ;+ ; NAME: ; SXT_TEEM1 ; PURPOSE: ; Compute Te, EM, d_Te, d_EM from the ratio of the two input data arrays. ; NOTE: To do this without use of index structures use sxt_teem2.pro. ; CALLING SEQUENCE: ; sxt_teem, index1, Image1, index2, Image2, Te ; sxt_teem, index1, Image1, index2, Image2, Te, EM, gain=gain, subs=subs ; sxt_teem, index1, Image1, index2, Image2, Te, EM, d_Te, d_EM, Valid ; sxt_teem, index1, Image1, index2, Image2, Te, EM, d_Te, d_EM, sum=sum ; INPUTS: ; index1 = index (structure) of Image1 ; Image1 = SXT image (DN), if byte type then sxt_prep is called. ; index2 = index (structure) of Image2 ; Image2 = SXT image (DN) through a different filter, if byte type then sxt_prep is called. ; NOTE: The Image variables can be either the raw compressed ; SXT images, simple scalars, or images that have already ; been processed through SXT_PREP. In the latter case ; the index structure created by SXT_PREP must be supplied. ; It is O.K. to input exposure-normalized PREP'd images ; although the d_te and d_em results will not be identical ; as when using raw compressed data because of round-off ; error. ; When comparing results using raw compressed vs. PREP'd ; input care must be taken to assure that the registration ; reference image is the same in both cases. ; OPTIONAL INPUT KEYWORDS: ; t1,t2 = Exposure times in msec. These values will override ; those in index1 and index2 if index1 and index2 are supplied. ; If index1 and t1 are not supplied, default to 1000 msec. ; gain = camera gain in e-/Dn. ; If not supplied, default= (~100). Used in EM calculation. ; interp = If set, use Spline interpolation of SXT response functions ; Default is linear interpolation. ; (NOTE: SPLINE [IDL user library] fails on large images.) ; sat1,sat2=Array of saturated pixels (if Image is not byte-type). ; unc1,unc2=Array of decompression uncertainties (if Image not byte-type). ; sum = Sum over sum X sum pixels before computing Te. ; satval = value of sat1 or sat2 to treat as a saturated pixel (def=1). ; average = If set, return average Te and total EM for entire image. ; thresh1,thresh2 = Minimum threshold (in DN) of background subtracted ; and rebinned (sum=X) image1, image2 ; date = Time in any format (including structure). Used to ; determine entrance filter transmission (not needed if ; index1 or index3 is an SXT index structure). ; OUTPUTS: ; Te = log10(Temp) (invalid temps = 0) ; OPTIONAL OUTPUTS: ; EM = log10(Emission measure). ; d_Te = Statistical uncertainties of Te, log10 format. ; d_EM = Statistical uncertainties of EM, log10 format. ; Valid = Array of valid pixels (=0). -1 = bad ratio values. ; OPTIONAL OUTPUT KEYWORDS: ; version = version number of input data base file ; COMMON BLOCKs: ; None. (In sxt_teem2, sxt_teem_db contains SXT response curves: sre*genx) ; RESTRICTIONS: ; o Image1 and Image2 must be the same size and 1-2 or 2-d (not a cube). ; PROCEDURE: ; o Will compute ratio of thicker/thinner to determine Te. ; o Calls sxt_prep with /reg, /dc_orbit_correct and /float if imageX is byte type ; o Use thinner image for EM determination. ; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ; Double Valued Functions: ; The following filter ratios are doubled valued (launch entr. filters) ; Al12/Noback, Al12/Al.1, AL12/AlMg, Al12/Mg3 ; Mg3 /Noback, Mg3 /Al.1, Mg3 /AlMg ; The temperatures returned may be a LOWER limit to the actual value. ; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ; MODIFICATION HISTORY: ; 10-Jul-91, Fe-Mei Lee Chou. Modified Jan 1992. SXT_TE was a function. ; 23-jan-93, G. Linford, Changed to a procedure, parameter changes. ; 28-jan-93, J. Lemen, Extensive mods., GAL & KTS minor typo correction. ; 5-feb-93, J. Lemen, J. McTiernan: Corrections and minor mods. ; 15-feb-93, J. Lemen, gt_expdur already contains 10% trans. ; 2-mar-93, J. Lemen, Add call to sxt_prep. Update calling arguments ; 11-mar-93, J. Lemen, Fixed a typo ; 23-mar-93, J. Lemen, Call sxt_teem1 to do the hard work. ; 24-mar-93, J. Lemen, Fixed some typo's ; 29-mar-93, JRL, Enable sum= option even if sxt_prep is not called ; internally. thresh1 now applied after background subtraction. ; Renamed to sxt_teem1 (sxt_teem1 renamed to sxt_teem2) ; 1-apr-93, JRL, Fixed a typo (error only if img2 was byte type). ; 18-nov-93, JRL, Fixed some data checking code. ; 22-dec-93, JRL, Fixed a IDL V3.1 related bug ; 25-jan-94, JRL, Add a warning message if the ND filter is used ; 11-feb-94, JRL, Fixed the warning message if the ND filter is used ; 25-jun-94, JRL, If all data is below the threshold, simply return 0s. ; If all data is saturated, return macro pixel value anyway. ; 15-apr-02, LWA, Added /float, /register, /dc_orbit_correct, and ; /dc_interpolate switches to sxt_prep. ; Deleted return of bogus Te, etc., values in saturated areas. ; Corrected code to work with sxt_prep normalized data. ; Renamed "arg" variables to "index". ; Updated header. ; 9-nov-2011, Aki Takeda, Added /mewe, /corona, /hybrid, /photospheric, ; and sre_file keywords. ;- -------------------------------------------------------------- ;on_error,1 ; Return to caller if there is problem thick_order = [1, 2, 3, 6, 5, 4] ; HARD-WIRED *** ; = Noback, Al.1, AlMg, Mg3, Al12, Be119 ;---------------------------------------------------------------------------- ; **** Step 1: Set up default values of all optional input parameters ;---------------------------------------------------------------------------- gain_default = 100. ; e/DN, gain at lunch if n_elements(gain) ne 0 then gain_ccd = gain else gain_ccd = gain_default if n_elements(sum) eq 0 then sum = 1 ; Summation mode if n_elements(satval) eq 0 then satval = 1 ; Threshold for sat. pixel if n_elements(interp) eq 0 then interp=0 ; Default is Linear interpolation if n_elements(thresh1) eq 0 then thresh1 = 0 ; Minimum threshold of image1 if n_elements(thresh2) eq 0 then thresh2 = 0 if n_elements(n_params1) eq 0 then n_params0 = n_params() else $ n_params0 = n_params1 ; Control calculations ; if no parameters are present, assume information mode if n_params0 lt 5 then begin doc_library,'sxt_teem1' print,format="(75('-'))" fastdoc,'sxt_teem1',/summ ; Give the full parameter list print,format="(75('-'))" return endif ;---------------------------------------------------------------------------- ; **** Step 2: Perform Checks on Input Data ;---------------------------------------------------------------------------- ; ***** ; Make sure image1 and image2 are scalar, 1-d or 2-d ; AND check image1 and image2 have the same dimensions: ; ***** sz1 = size([image1]) sz2 = size([image2]) if ((sz1(0) gt 2) or (sz2(0) gt 2)) or $ (sz1(0) ne sz2(0)) or $ (sz1(1) ne sz2(1)) or $ (n_elements(image1) ne n_elements(image2)) then begin message,' ** Error ** image1 and image2 must be 1-d or 2-d',/cont help,image1,image2 & tbeep message,' Size/Length of image1 and image2 must match' endif ; ***** ; determine type of input ; ***** date_d = '31-aug-91' ; Default date sz0=size(index1) ; check args 1 and 3 sz3=size(index2) fila1 = gt_filta(index1) ; Filter1 A t1 = gt_expdur(index1) ; exposure time in msec date_d= gt_day(index1,/str) ; Get the date filt1_id=gt_filtb(index1) ; Filter1 B if fila1 eq 6 then begin ; tbeep,5 ; Ring the bell message,'ND filter for exposure at: '+fmt_tim(index1),/cont message,' Increased scatter may lead to erroneous result',/cont endif fila2 = gt_filta(index2) ; Filter A t2 = gt_expdur(index2) ; exposure time in msec date_d= gt_day(index2,/str) ; Get the date filt2_id= gt_filtb(index2) ; Filter B if fila2 eq 6 then begin message,'ND filter for exposure at: '+fmt_tim(index2),/cont message,' Increased scatter may lead to erroneous result',/cont endif if n_elements(tt1) ne 0 then t1 = tt1 ; Keyword overrides if n_elements(tt2) ne 0 then t2 = tt2 ; Keyword overrides if n_elements(date) eq 0 then date=date_d ; ***** ; Data decompression and background subtraction ; ***** ; if data BYTE-type, assume it is compressed. ; if data is BYTE-type and the index is supplied, subtract background if(sz1(sz1(0)+1) eq 1) then begin ; Yes -- This is a byte array sxt_prep,[index1,index2],[[[image1]],[[image2]]],$ nind,img,uunc,ssat,sum=sum,error=error,/float,/reg,$ /dc_orbit_correct,/dc_interpolate img1 = reform(img(*,*,0)) uunc1 = reform(uunc(*,*,0)) sig1 = img1 img2 = reform(img(*,*,1)) uunc2 = reform(uunc(*,*,1)) sig2 = img2 satpix = reform(ssat(*,*,0)) + reform(ssat(*,*,1)) gt 0 ;Inclusive saturated pixel array delvarx,img,uunci,ssat ;Finished with these. endif else begin ; Not a byte array img1 = image1 img2 = image2 ; If image1 is not byte-type, make guess for decompression error (unc1): ; LWA 4/12/02, added exposure correction for normalized data. if gt_expdur(index1) eq 1000. then begin case 1 of gt_res(index1) eq 0 : pixarea=1 gt_res(index1) eq 1 : pixarea=4 gt_res(index1) eq 2 : pixarea=16 endcase sig1=img1*pixarea*gt_expdur(index1,/original)/1000. endif else begin sig1=img1 ;To pass on to sxt_teem2. endelse if (n_elements(unc1) eq 0) then begin dum=sxt_decomp(sxt_comp(sig1),uunc1) endif else begin uunc1 = unc1 endelse if n_elements(sat1) eq 0 then ssat1 = byte(img1*0) else ssat1 = sat1 if sum gt 1 then sxt_sumxy,index1,img1,sum,sum,uunc1,ssat1,/silent ; If image2 is not byte-type, make guess for decompression error (unc2): if gt_expdur(index2) eq 1000. then begin case 1 of gt_res(index2) eq 0 : pixarea=1 gt_res(index2) eq 1 : pixarea=4 gt_res(index2) eq 2 : pixarea=16 endcase sig2=img2*pixarea*gt_expdur(index2,/original)/1000. endif else begin sig2=img2 ;To pass on to sxt_teem2. endelse if (n_elements(unc2) eq 0) then begin dum=sxt_decomp(sxt_comp(img2),uunc2) endif else begin uunc2 = unc2 endelse if n_elements(sat2) eq 0 then ssat2 = byte(img2*0) else ssat2 = sat2 if sum gt 1 then sxt_sumxy,index2,img2,sum,sum,uunc2,ssat2,/silent satpix = ssat1 > ssat2 endelse ; ***** ; check that front filter position is valid ; ***** if ((fila1 ne 1) and (fila1 ne 6)) or ((fila2 ne 1) and (fila2 ne 6)) then begin message,' ** Error ** Check the following: ',/cont if ((fila1 ne 1) and (fila1 ne 6)) then invalid1 = 'Invalid' else invalid1='' if ((fila2 ne 1) and (fila2 ne 6)) then invalid2 = 'Invalid' else invalid2='' print,'Filter 1 B/A= ',filt1_id,'/',fila1,invalid1,format='(a,i2,a,i2,2x,a)' print,'Filter 2 B/A= ',filt2_id,'/',fila2,invalid2,format='(a,i2,a,i2,2x,a)' message,' Invalid Front Filter (A)' endif ; ***** ; check that rear filter position is valid ; ***** if (filt1_id gt 6) or (filt2_id gt 6) or $ (filt1_id lt 1) or (filt2_id lt 1) then begin message,'** Error ** No such filter in database',/cont message,string(' Check the following: filt1=',filt1_id,' filt2=',filt2_id, $ format='(a,i3,a,i3)') endif if (filt1_id eq filt2_id) then message,string('** Error ** Filter1 = Filter2 = ',filt1_id,format='(a,i3)') ; ***** ; Find the order of thickness of the corresponding filter ; ***** filt1_th = thick_order(filt1_id-1) ; thickness order for filter 1 filt2_th = thick_order(filt2_id-1) ; thickness order for filter 2 if keyword_set(average) then begin pix = where((img1 ge thresh1) and (img2 ge thresh2),npix) npix2= 0 endif else begin pix = where((img1 ge thresh1) and (img2 ge thresh2) and $ (img1 gt 0) and (img2 gt 0) and $ (satpix lt satval),npix) ; Pixels to calculate ratio(thick/thin) ; Saturated pixel array: pix2= where(satpix ge satval,npix2) ; Saturated pixels endelse if npix eq 0 then begin szz = sz1 szz(szz(0)+1) = 4 ; Floating array Te = make_array(size=szz) ; Return a floating array if n_params() ge 6 then EM = Te ; Return EM if requested if n_params() ge 7 then d_Te = Te ; Return d_Te if requested if n_params() ge 8 then d_EM = Te ; Return d_EM if requested if n_params() eq 9 then begin ; Return Valid if requested szz(szz(0)+1) = 2 ; Set to I*2 valid = make_array(size=szz) - 1 ; Set to -1 endif endif if (npix+npix2 eq 0) then begin message,' ** Error ** no valid data',/cont return endif ; ***** ; If average keyword set, total the images. ; ***** if keyword_set(average) then begin ; Includes all saturated pixels img1 = total(img1(pix)) uunc1 = sqrt(float(uunc1(pix))^2) / n_elements(uunc1(pix)) img2 = total(img2(pix)) uunc2 = sqrt(float(uunc2(pix))^2) / n_elements(uunc2(pix)) ppix = 0 ; Only one pixel to calculate endif else begin ppix = pix ; Pixels to calculate ratio of thick/thin endelse ;---------------------------------------------------------------------------- ; **** Step 3: Call SXT_TEEM2 to calculate Te, Em, d_Te, d_EM ;---------------------------------------------------------------------------- if(filt1_th gt filt2_th) then begin ; filt1 is thicker t0=gt_expdur(index2,/original) sxt_teem2,filt2_id,img2,filt1_id,img1,te,em,d_te,d_em,valid, $ t1=t2, t2=t1, interp=interp, gain_ccd=gain_ccd, $ date=date, pix=ppix, n_params0=n_params0, $ unc1=uunc2, unc2=uunc1, version=version, $ sig1=sig1, sig2=sig2, t0=t0, $ mewe=mewe, corona=corona, hybrid=hybrid, $ photospheric=photospheric, sre_file=sre_file endif else begin ; filt2 is thicker t0=gt_expdur(index1,/original) sxt_teem2,filt1_id,img1,filt2_id,img2,te,em,d_te,d_em,valid, $ t1=t1, t2=t2, interp=interp, gain_ccd=gain_ccd, $ date=date, pix=ppix, n_params0=n_params0, $ unc1=uunc1, unc2=uunc2, version=version, $ sig1=sig1, sig2=sig2, t0=t0, $ mewe=mewe, corona=corona, hybrid=hybrid, $ photospheric=photospheric, sre_file=sre_file endelse ;---------------------------------------------------------------------------- ; **** Step 4: Assign the results of the saturated pixel ;---------------------------------------------------------------------------- if n_params0 ge 9 then begin ; Not necessary if Valid was not requested if keyword_set(average) then begin if valid ne -1 then begin valid = 0*fix(image1) -1 ; Assign -1's initially valid(pix) = satpix(pix) ; Show where the saturated pixels are endif else valid = 0*fix(image1) -1 ; Return all -1's endif else begin jj = where(valid eq 0,njj) if njj gt 0 then valid(jj) = satpix(jj) ; Fill in saturated pixel info endelse endif ; The saturated pixel was treated as a "macro" pixel ; LWA 4/12/02, delete all this and return zeroes in saturated areas. end