PRO 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, $ ; Input date=date, pix=pix, n_params0=n_params0, $ ; Input unc1=unc1,unc2=unc2,t0=t0,sig1=sig1,sig2=sig2, $ ; Input version=version ; Return ;+ ; NAME: ; SXT_TEEM2 ; PURPOSE: ; Low level SXT temperature routine (called by sxt_teem1) ; Temperature from img2/img1 and Em from img1 ; (img1 should from the thinner of the two filters) ; CALLING SEQUENCE: ; sxt_teem2, filt1, img1, filt2, img2, Te , EM ; sxt_teem2, filt1, img1, filt2, img2, Te, EM, d_te, d_em, valid, $ ; t1=t1, t2=t2, interp=interp, gain_ccd=gain_ccd $ ; date=date, pix=pix, n_params0=n_params0, $ ; unc1=unc1, unc2=unc2, version=version ; INPUTS: ; filt1 = Filter B value (1-5) of Img1 ; Img1 = SXT counts (DN) decompressed, registered and background subtracted ; filt2 = Filter B value (1-5) of Img2 ; Img2 = SXT counts (DN) decompressed, registered and background subtracted ; ; OPTIONAL INPUT KEYWORDS: ; t1,t2 = Exposure times in msec. (Default is 1000 ms). ; date = Date to determine fraction of launch entrance filter ; for response function (Default = '31-aug-91') ; pix = Array of pixels to compute observed ratio. ; ccd_gain= camera gain in e-/Dn. (Default is 100). ; interp = If set, use Spline interpolation of SXT response functions ; (default is INTERPOL interpolation). ; unc1,unc2=Array of decompression uncertainties (default is 0) ; n_params0=To override the number of actual parameters. This controls ; what calculations are done: ; 5 for Te, 6 for Te+EM, 7 for Te+EM+dTE+dEM ; OUTPUTS: ; Te = log10(Temp) (invalid temps = 0) ; OPTIONAL OUTPUTS: ; EM = log10(Emission measure). ; d_Te = Statistical uncertainties of Te, (log10 units) ; d_EM = Statistical uncertainties of EM ; (includes uncertainty of d_Te) (log10 units) ; Valid = Array of valid pixels (=0). -1 = bad ratio values. ; OPTIONAL OUTPUT KEYWORDS: ; version = version number of input data base file ; COMMON BLOCKs: ; sxt_teem_db Contains the SXT response functions (sre*genx) ; RESTRICTIONS: ; o Img1 and Img2 must be the same size, decompressed, registered, ; straylight and background corrected. ; o Not many checks are made -- assume this was done by sxt_teem1.pro ; o Will compute ratio of img2/img1 to determine Te. ; The caller should make sure that img1 is the thinner filter. ; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ; 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: ; 18-mar-93, J. Lemen, Written ; 23-mar-93, J. Lemen, Added modified entrance filter code ; 29-mar-93, JRL, Don't return dummy unc1, unc2 arrays (call delvarx) ; 18-nov-93, JRL, Fixed bug with not including Mg3 filter ; 22-dec-93, JRL, Fixed an IDL V3.1 related bug ; 23-mar-94, RDB, Made file name sre*.* (VMS compatibility) ; 25-jun-94, JRL, Call get_yo_dates to get open entrance filter fracton ; Return correctly dimensioned Te, EM even if there are no valid pixels. ; 31-aug-94, JMM, added Klimchuk correction to EM uncertainty ; calculation, includes temperature uncertainty ; in EM uncertainty calculation. ; 15-Apr-02, LWA Modified to pass total signal array and observing ; time from sxt_teem1 for those cases using sxt_prep normalized input. ; 19-Apr-02, LWA Corrected case for t1 not specified. ;- ; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ; COMMON BLOCKs: ; sxt_teem_db Contains the SXT response functions (sre*genx) ; common sxt_teem_db, db, text, header, version1, Filein ; db = array of (n_filter+1,n_temp) ; first colum is temperature(log value) ; rest of colum is number of electrons ; for corresponding filters(not filter id) ; (n_temp is number of temperatures) ; text = some information about the database. ; header = time and date of the creation of database. ; version1 = version number of the database file.(from database ; file) ; filein = Name of the file ; Additional: ; em_assumed = assumed em value.(from database file) ; label_id = filter id number corresponding to col. number ; if label_id(4) = 6 means that id number 4 if at ; col. number 4(fourth filter) in database data. ; This information is from database file. ; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - common sxt_teem_db, db, text, header, version1, Filein ; ------------------------ ; HARD-WIRED Assumptions: ; ------------------------ em_assumed = 44. + 3 ; for /msec and EM=1.e44; HARD-WIRED *** label_id = [1, 2, 3, 6, 5, 4] ; HARD-WIRED *** thick_order = [1, 2, 3, 6, 5, 4] ; HARD-WIRED *** gain_default = 100. ; e/DN, gain at lunch ; if no parameters are present, assume information mode if n_params() lt 5 then begin doc_library,'sxt_teem2' print,format="(75('-'))" fastdoc,'sxt_teem2',/summ ; Give the full parameter list print,format="(75('-'))" return endif ;---------------------------------------------------------------------------- ; **** Step 1: Set up default values of all optional input parameters: ;---------------------------------------------------------------------------- if n_elements(gain_ccd) eq 0 then gain_ccd = gain_default if n_elements(date) eq 0 then date = '31-aug-91' ; Def = Launch date if n_elements(interp) eq 0 then interp=0 ; Def = Linear interpolation if n_elements(n_params0) eq 0 then n_params0 = n_params() if n_elements(t1) ne 0 then begin if t1 ne 1000. then t0 = t1 endif if n_elements(t1) eq 0 then begin t1 = 1000. & t0 = t1 endif if n_elements(t2) eq 0 then t2 = 1000. ; 1000 ms if n_elements(sig1) eq 0 then sig1 = img1 if n_elements(sig2) eq 0 then sig2 = img2 t1 = float(t1) & t2 = float(t2) ; Make sure it is floating Te = 0. & EM = 0. & d_Te= 0. & d_EM= 0. & Valid = -1 ; In case of early return ; **** ; Only compute for img1 and img2 gt 0 ; **** if n_elements(pix) eq 0 then begin pix_new = 1 ; Reset pix before return pix = where((img1 gt 0) and (img2 gt 0)) if pix(0) eq -1 then begin message, ' *** Error : (1) no valid data',/cont & tbeep delvarx,pix return endif endif else pix_new = 0 ; Don't reset at the end ;---------------------------------------------------------------------------- ; **** Step 2: Read the Data Base File ;---------------------------------------------------------------------------- ; Read data from file only for first time call to this function. ; Save the data in common block sxt_teem_db. if(n_elements(db) eq 0) then begin ; sre = electrons vs Te filein = findfile(concat_dir('$DIR_SXT_SENSITIVE','sre*.*')) filein = filein(n_elements(filein)-1) ; Use latest version print,'Reading data base file = ',filein restgen, str=db,file=filein, text=text, header=header break_file, filein, disk, dir, filnam version1 = strmid(filnam,strpos(filnam,'sre')+3,strlen(filnam)) endif version = version1 ; Return to the caller ;---------------------------------------------------------------------------- ; **** Step 3: Perform Checks on Input Data ;---------------------------------------------------------------------------- ; ***** ; determine type of input ; ***** sz1=size(img1) ; check args 1 and 3 sz2=size(img2) if (sz1(sz1(0)+1) eq 1) or (sz2(sz2(0)+1) eq 1) or $ ; Check for byte type (sz1(sz1(0)+2) ne sz2(sz2(0)+2)) then begin $ ; Make sure equal length message,' ** Error **',/cont & tbeep message,' Images must be decompressed and equal length' endif ; ***** ; check that rear filter position is valid ; ***** n_filter = (size(db.elects))(1)/2 ; number of filters in database if((filt1_id gt n_filter) or (filt2_id gt n_filter)) then $ message, 'No such a filter in database' ; ***** ; Set up the data base vectors ; ***** filt1_col = (where(label_id eq filt1_id))(0) ; colum number in database filt2_col = (where(label_id eq filt2_id))(0) ; colum number in database te_db = reform(db.Temp) ; temperature array in database ; Set up Launch values first: i1_db = transpose(db.elects(filt1_col,*)) ; electrons for img1 i2_db = transpose(db.elects(filt2_col,*)) ; electrons for img2 i1_ph = transpose(db.photons(filt1_col,*)) ; photons for img1 i2_ph = transpose(db.photons(filt2_col,*)) ; photons for img2 ; Modify for change in entrance filter: Frac_open = get_yo_dates(date,/ent,/val) ; Get open area fraction of entrance filter if Frac_open gt 0. then begin print,' Open fraction of entrance filter =',frac_open,format='(a,f6.3)' i1_db = i1_db*(1-Frac_open) + Frac_open*transpose(db.elects(6+filt1_col,*)) i2_db = i2_db*(1-Frac_open) + Frac_open*transpose(db.elects(6+filt2_col,*)) i1_ph = i1_ph*(1-Frac_open) + Frac_open*transpose(db.photons(6+filt1_col,*)) i2_ph = i2_ph*(1-Frac_open) + Frac_open*transpose(db.photons(6+filt2_col,*)) endif ;---------------------------------------------------------------------------- ; **** Step 4: Compute Ratio of img1/img2 ;---------------------------------------------------------------------------- ; **** ; Calculate r_db and ratio_r = img1 / img2 ; **** ratio_r = float(img1) * 0. Te = float(img1) * 0. 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(filt1_th gt filt2_th) then begin ; filt1 is thicker message,'** Warning ** ',/cont & tbeep print,' Filter 1 is thicker than Filter 2 = ',filt1_id,filt2_id message,' Solution may not be found',/cont,/noname endif pidb = where(i1_db gt 0.) te_db = te_db(pidb) r_db = i2_db(pidb) / i1_db(pidb) ratio_r(pix) = (img2(pix) * (t1 / t2)) / img1(pix) filt1_str = strtrim(gt_filtb(filt1_id,/str)) filt2_str = strtrim(gt_filtb(filt2_id,/str)) ph_dn1 = gain_ccd * i1_ph(pidb) / i1_db(pidb) ; Number of photons per dn ph_dn2 = gain_ccd * i2_ph(pidb) / i2_db(pidb) ; **** ; Check to see if r_db vs Te_db is double valued or not ; **** xx = max(r_db) & xx = !c if xx ne N_elements(r_db)-1 then begin message,' ** Warning **',/cont & tbeep print, ' Ratio of ',filt2_str,' / ',filt1_str,' is double valued' Te_db = Te_db(0:xx) r_db = r_db(0:xx) endif ; **** ; find the corresponding temperature for each ratio_r from te_db and r_db ; **** v_dbr = where(r_db gt 0.) ; Valid Data Base ratio Te_db = Te_db(v_dbr) & r_db = r_db(v_dbr) ; Only valid data base ph_dn1 = ph_dn1(v_dbr) ; Number of photons per dn ph_dn2 = ph_dn2(v_dbr) ; Number of photons per dn v_odr = lonarr(n_elements(img1)) ; Set up vector of valids v_odr(pix) = 1 v_odr = where((v_odr eq 1) and $ (ratio_r ge min(r_db)) and (ratio_r le max(r_db)),nv_odr) if pix_new eq 1 then delvarx,pix ; Reset the pix array if (nv_odr eq 0) then begin message, ' ** Error **: (2) no valid data',/cont & tbeep if n_params0 ge 6 then EM = Te if n_params0 ge 7 then d_Te = Te if n_params0 ge 8 then d_EM = Te if n_params0 eq 9 then Valid = fix(Te) - 1 return endif else begin ; use dspline (deluxe spline- calls spline or interpol) to interpolate data te(v_odr) = dspline(alog10(r_db), te_db, alog10(ratio_r(v_odr)),interp=interp) endelse ;---------------------------------------------------------------------------- ; **** Step 5: Compute Emission Measure: EM ;---------------------------------------------------------------------------- ; use sig1 for em calculation if n_params0 ge 6 then begin ; Return emmission meas. em = 0. * float(sig1) ; Setup output variable te_db1 = transpose(db.Temp) ; temperature array in database ff_te = Te ff_te(v_odr) = dspline(te_db1(pidb), i1_db(pidb), Te(v_odr),interp=interp) em(v_odr) = alog10(((gain_ccd / t0) * sig1(v_odr)) / ff_te(v_odr)) + em_assumed endif ;---------------------------------------------------------------------------- ; **** Step 6: Compute Uncertainties. ;---------------------------------------------------------------------------- if n_params0 ge 7 then begin ; Return Uncertainties on Te and EM ; Set up some variables: ph1 = Te*0. & ph2 = ph1 & ud1 = ph1 & ud2 = ph1 & pd1 = ph1 & pd2 = ph1 if n_elements(unc1) eq 0 then begin unc1_new = 1 & unc1 = byte(img1) * 0 endif else unc1_new = 0 if n_elements(unc2) eq 0 then begin unc2_new = 1 & unc2 = byte(img2) * 0 endif else unc2_new = 0 ; Compute photons for each filter: pd1(v_odr) = 10.^(dspline(Te_db, alog10(ph_dn1), Te(v_odr),interp=0)) ; photons/dn ph1(v_odr) = (sig1(v_odr) * pd1(v_odr)) > 1.0 ; Number of ph pd2(v_odr) = 10.^(dspline(Te_db, alog10(ph_dn2), Te(v_odr),interp=0)) ; photons/dn ph2(v_odr) = (sig2(v_odr) * pd2(v_odr)) > 1.0 ; Number of ph ; Uncertainty in dn: ud1(v_odr) = sqrt(ph1(v_odr)/(pd1(v_odr))^2 + float(unc1(v_odr))^2 + 1) ud2(v_odr) = sqrt(ph2(v_odr)/(pd2(v_odr))^2 + float(unc2(v_odr))^2 + 1) if unc1_new then delvarx,unc1 if unc2_new then delvarx,unc2 delvarx,ph1,ph2,pd1,pd2 ; Release some memory ; Sigma_R = sigma_Ratio(obs)/Ratio(obs) sigma_R = Te * 0. sigma_R(v_odr) = sqrt((ud1(v_odr)/sig1(v_odr))^2 + $ (ud2(v_odr)/sig2(v_odr))^2) ; Compute d log10(R) / d log10(Te): dlogR_dlogT = deriv(Te_db, alog10(r_db)) ; Userlib routine ; Interpolate to Te values: d_Te = Te * 0. d_Te(v_odr) = alog10(exp(1.)) * $ sigma_R(v_odr) / dspline(Te_db, dlogR_dlogT, Te(v_odr) ,interp=0) ; Compute uncertainty in EM: ;JMM 31-aug-94 d_EM = 0. * sig1 d_EM(v_odr) = ud1(v_odr) / (sig1(v_odr) *alog(10.)) ; Compute d log10(ff) / d log10(Te) * d_Te, the Klimchuk factor dlogff_dlogT = deriv(te_db1(pidb), alog10(i1_db(pidb))) dlffdlT = dspline(te_db1(pidb), dlogff_dlogT, Te(v_odr), interp = 0) d_EM(v_odr) = sqrt(d_EM(v_odr)^2 + (d_Te(v_odr)*dlffdlT)^2) ; Release some memory: ud1 = 0 & ud2 = 0 ENDIF ;---------------------------------------------------------------------------- ; **** Step 7: Assign the Valid array ;---------------------------------------------------------------------------- Valid = fix(Te*0) - 1 ; Fill the array with -1 Valid(v_odr) = 0 ; Valids are = 0 end