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 unc_p1=unc_p1,unc_p2=unc_p2, $ ; Input corona=corona, photospheric=photospheric, hybrid=hybrid, $ ; Input mewe=mewe, sre_file=sre_file, $ ; 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, sre_file=sre_file ; 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, sre_file=sre_file, $ ; version=version ; INPUTS: ; filt1 = Filter B value (1-6) of Img1 ; Img1 = SXT counts (DN) decompressed, registered and background subtracted ; filt2 = Filter B value (1-6) of Img2 ; Img2 = SXT counts (DN) decompressed, registered and background subtracted ;**** the following keywords are for selecting the response function. **** ; corona = if set, use sre_ch601_corona_chianti.genx ; hybrid = if set, use sre_ch601_hybrid_chianti.genx (default) ; photospheric = if set, sre_ch601_photos_chianti.genx ; mewe = if set, former SXT response, sre950419_02.genx, is used. ; sre_file = specify the response function file name with path. ; ; 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) ; unc_p1,unc_p2= Photon noise (statistical uncertainty) for the img1 and ; img2, calculated in other place and to be passed to ; this program. If not set, calculated from the db with ; no consideration on integrated or averaged DN. ; 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 fraction ; 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. ; 15-Nov-2004 - S.L.Freeland - use sxt_teem_common.pro for common def. ; 7-Jul-2005 LWA Commented out "Open fraction ..." print statement. ; 20-Feb-2010 Aki Takeda, - change the default sre file ; (Mewe to Chianti with coronal abundance). ; - added SRE_FILE keyword to accept user-specified ; response functions. ; - setting hard-wired variable in a separate ; routine, teem_hardwired.pro ; 06-Aug-2010 Aki T, - added UNC_P1, UNC_P2 keyword to pass photon noise ; calculated outside this program. ; 28-Sep-2010 Aki T, - update the path of the default sre file ; (ydb_beta -> $SSW_SXT/response). ; 18-Mar-2011 Aki T, - added /corona, /photospheric, /hibrid ; and /mewe keywords. ; 27-Jan-2018 Aki T, - update response functions to chainti 800 version. ; also changed binding to abundance files. ; /hibrid : sun_coronal_2012_schmelz_ext.abund (from Fludra 1999) ; /photos : sun_photospheric_2011_caffau.abund (from Grevesse 2007) ; /corona : sun_coronal_1992_feldman_ext.abund (no change) ;- ; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ; 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. ; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @sxt_teem_common ;common sxt_teem_db, db, text, header, version1, Filein ; ------------------------ ; HARD-WIRED Assumptions: ; ------------------------ @teem_hardwired ;(current hard wired parameters, checked 22-Feb-2010) ;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 chkarg,'sxt_teem2' ; 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 when the sre* file name is different from the one ; in the common block. Save the new data in common block sxt_teem_db. dir_resp = '$SSW_SXT/response/' if (n_elements(sre_file) eq 0) then begin ; set default sre_file = dir_resp+'sre_ch601_hybrid_chianti.genx' ; (with hybrid abund. & chianti.ioneq ) endif ; if keyword_set(corona) then sre_file = dir_resp+'sre_ch800_corona_chianti.genx' if keyword_set(hybrid) then sre_file = dir_resp+'sre_ch800_hybrid_chianti.genx' if keyword_set(photospheric) then sre_file = dir_resp+'sre_ch800_photos_chianti.genx' ;if keyword_set(corona) then sre_file = dir_resp+'sre_ch601_corona_chianti.genx' ;if keyword_set(hybrid) then sre_file = dir_resp+'sre_ch601_hybrid_chianti.genx' ;if keyword_set(photospheric) then sre_file = dir_resp+'sre_ch601_photos_chianti.genx' if keyword_set(mewe) then sre_file = dir_resp+'sre950419_02.genx' ; ; print,'*** SXT response file '+sre_file+' is used.' ; go_read=0 if (n_elements(filein) eq 0) then go_read=1 if (n_elements(filein) ge 1) then if (filein ne sre_file) then go_read=1 if (go_read eq 1) then begin print,'Reading data base file = ', sre_file restgen, str=db,file=sre_file, text=text, header=header filein=sre_file break_file, filein, disk, dir, filnam version1 = strmid(filnam,strpos(filnam,'sre')+3,strlen(filnam)) if strmid(version1,0,1) eq '_' then $ ; handling new sre file convention. version1 = strmid(version1,strpos(version1,'_',/reverse_search)+1,8) 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 ;----- new KW, unc_p1, unc_p2 handlnig (Aug-2010, AkT) ------ if n_elements(unc_p1) ne 0 then unc_p1sq=unc_p1^2 $ else begin ; Compute photons for the filter 1: 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 unc_p1sq=ph1(v_odr)/(pd1(v_odr))^2 ; Uncertainty in dn endelse if n_elements(unc_p2) ne 0 then unc_p2sq=unc_p2^2 $ else begin ; Compute photons for the filter 2: 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 unc_p2sq=ph2(v_odr)/(pd2(v_odr))^2 ; Uncertainty in dn endelse ud1(v_odr) = sqrt(unc_p1sq + float(unc1(v_odr))^2 + 1) ud2(v_odr) = sqrt(unc_p2sq + float(unc2(v_odr))^2 + 1) ;------------------------------------------------------------ ; 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