PRO SXT_TEEM1, arg1, image1, arg3, 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 ; Returned ;+ ; NAME: ; SXT_TEEM1 ; PURPOSE: ; Compute temperature from the ratio of the two input data arrays. ; CALLING SEQUENCE: ; sxt_teem1, index1, Image1, index2, Image2, Te [, EM, gain=gain, /interp] ; sxt_teem1, filt1, Image1, filt2, Image2, Te, EM, t1=t1, t2=t2 ; sxt_teem1, filt1, Image1, index2, Image2, Te, thresh1=thresh1 ; sxt_teem1, index1, Image1, index2, Image2, Te, EM, d_Te, d_EM, Valid ; ; INPUTS: ; index1 = index (structure) of Image1 ; filt1 = Filter B value (1-5) of Image1 ; Image1 = SXT counts as DN ; index2 = index (structure) of Image2 ; filt2 = Filter B value (1-5) of Image2 ; Image2 = SXT counts as DN ; ; 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 ; d_EM = Statistical uncertainties of EM ; Valid = Array of valid pixels (=0). -1 = bad ratio values. ; 1 indicates a saturated pixel (or possibly > 1 if sum>1) ; 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 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. ; When arg1 or arg3 is a filter number, make sure it is a scalar. ;- -------------------------------------------------------------- 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(arg1) ; check args 1 and 3 sz3=size(arg3) if sz0(sz0(0)+1) eq 8 then begin ; Yes -- arg1 is a structure (index) fila1 = gt_filta(arg1) ; Filter1 A t1 = gt_expdur(arg1) ; exposure time in msec date_d= gt_day(arg1,/str) ; Get the date filt1_id=gt_filtb(arg1) ; Filter1 B if fila1 eq 6 then begin ; tbeep,5 ; Ring the bell message,'ND filter for exposure at: '+fmt_tim(arg1),/cont message,' Increased scatter may lead to erroneous result',/cont endif endif else begin t1 = 1000. ; Default is 1000. msec fila1 = 1 ; Default is open if n_elements(arg1) ne 1 then message,'*** ERROR: Filter(1) must be a scalar ***' filt1_id=(gt_filtb(arg1))(0) ; Filter1 B endelse if sz3(sz3(0)+1) eq 8 then begin ; Yes -- arg3 is a structure (index) fila2 = gt_filta(arg3) ; Filter A t2 = gt_expdur(arg3) ; exposure time in msec date_d= gt_day(arg3,/str) ; Get the date filt2_id= gt_filtb(arg3) ; Filter B if fila2 eq 6 then begin ; tbeep,5 ; Ring the bell message,'ND filter for exposure at: '+fmt_tim(arg3),/cont message,' Increased scatter may lead to erroneous result',/cont endif endif else begin t2 = 1000. ; Default is 1000. msec fila2 = 1 ; Default is open if n_elements(arg3) ne 1 then message,'*** ERROR: Filter(2) must be a scalar ***' filt2_id= (gt_filtb(arg3))(0) ; Filter B endelse 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,arg1,image1,nind,img1,uunc1,ssat1,sum=sum,error=error if error then return endif else begin ; Not a byte array img1 = image1 ; If image1 is not byte-type, make guess for decompression error (unc1): if (n_params0 ge 7) then $ if (n_elements(unc1) eq 0) then $ dum=sxt_decomp(sxt_comp(img1),uunc1) else uunc1 = unc1 if n_elements(sat1) eq 0 then ssat1 = byte(img1*0) else ssat1 = sat1 if sum gt 1 then sxt_sumxy,arg1,img1,sum,sum,uunc1,ssat1,/silent endelse if(sz2(sz2(0)+1) eq 1) then begin ; Yes -- This is a byte array sxt_prep,arg3,image2,nind,img2,uunc2,ssat2,sum=sum,error=error if error then return endif else begin ; Not a byte array img2 = image2 ; If image2 is not byte-type, make guess for decompression error (unc2): if (n_params0 ge 7) then $ if (n_elements(unc2) eq 0) then $ dum=sxt_decomp(sxt_comp(img2),uunc2) else uunc2 = unc2 if n_elements(sat2) eq 0 then ssat2 = byte(img2*0) else ssat2 = sat2 if sum gt 1 then sxt_sumxy,arg3,img2,sum,sum,uunc2,ssat2,/silent endelse satpix = ssat1 > ssat2 ; Inclusive saturated pixel array delvarx,ssat1,ssat2 ; Finished with ssat1 and ssat2 ; ***** ; 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 if npix gt 0 then $ 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 if npix2 gt 0 then $ sxt_teem2,filt2_id,total(img2(pix2)),filt1_id,total(img1(pix2)),$ te1,em1,d_te1,d_em1,valid1, $ t1=t2, t2=t1, interp=interp, gain_ccd=gain_ccd, $ date=date, n_params0=n_params0, $ unc1=uunc2, unc2=uunc1 endif else begin ; filt2 is thicker if npix gt 0 then $ 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 if npix2 gt 0 then $ sxt_teem2,filt1_id,total(img1(pix2)),filt2_id,total(img2(pix2)),$ te1,em1,d_te1,d_em1,valid1, $ t1=t1, t2=t2, interp=interp, gain_ccd=gain_ccd, $ date=date, n_params0=n_params0, $ unc1=uunc1, unc2=uunc2 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 if npix2 gt 0 then begin ; Fill in saturated results te(pix2) = te1 if n_params0 ge 6 then em(pix2) = em1 - alog10(npix2); Assign average EM / pixel if n_params0 ge 7 then d_te(pix2) = d_te1 if n_params0 ge 8 then d_em(pix2) = d_em1 ; This might not be the right way if n_params0 ge 9 then Valid(pix2) = satpix(pix2) ; Saturated valids are > 0 endif end