PRO xrt_teem4syn_ch1102_slmd_fctm, index1, data1, index2, data2, $ te, em, et, ee, cal = cal, $ bin = bin, trange = trange, no_threshold = no_threshold, $ te_err_threshold = te_err_threshold, photon_noise_threshold = photon_noise_threshold, $ mask1 = mask1, mask2 = mask2, $ lambda_spectrum = lambda_spectrum, te_spectrum = te_spectrum, spectrum = spectrum, $ apac_spectrum = apac_spectrum, eff_area = eff_area, $ corona=corona, photospheric=photospheric, $ frac_open=frac_open, f1ctm_fact=f1ctm_fact,f2ctm_fact=f2ctm_fact ; ========================================================================= ;+ ; PROJECT: ; ; Solar-B / XRT ; ; NAME: ; ; XRT_TEEM (Commented out history checking part to process ; synoptic composite images whose intensities are normalized. ; Modified index (exptime changed to 1.0) is used as input. ; Dummy 2x2 data arrays are used to avoid summation over time. ; Also, LAMBDA_SPECTRUM, TE_SPECTRUM, SPECTRUM are ; re-defined by CHIANTI ver. 10.1 spectra. /hybrid keyword ; was dropped due to Chianti 10.1's new abundance files. ; Taking partially open pre-filter transmission into account) ; ; CATEGORY: ; ; XRT CORONAL TEMPERATURE DIAGNOSTICS ; ; PURPOSE: ; ; Get coronal temperature using filter ratio method. ; Temperature and volume emission measure is derived in log scale. ; ; In default, this program use the solar spectrum ; calculated with CHIANTI database ver. 6.0.1 ; (density: 10^9 [cm^-3], ionization equilibrium: chianti.ioneq, abundance: sun_coronal_ext). ; ; CALLING SEQUENCE: ; ; XRT_TEEM, index1, data1, index2, data2, te, em, et, ee, [bin=bin], [te_err_threshold=te_err_threshold], $ ; [photon_noise_threshold=photon_noise_threshold], [/no_threshold], [/info], $ ; [lambda_spectrum = lambda_spectrum, te_spectrum = te_spectrum, spectrum = spectrum], $ ; [/apac_spectrum], [eff_area = eff_area] ; ; INPUTS: ; ; INDEX1 - (structure array) XRT index. ; DATA1 - (2 or 3-dim float array, [x, y, t]) XRT *non-normalized* level1 data. ; INDEX2 - (structure array) XRT index. ; DATA2 - (2 or 3-dim float array, [x, y, t]) XRT *non-normalized* level1 data. ; FRAC_OPEN - (float) Open fraction of pre-filter (0.0 :no opening -- 1.0 : fully open) ; F1CTAM_FACT - (float) Filter contamination enhancement factor for filter1 (default 1.0 : no correction) ; F2CTAM_FACT - (float) Filter contamination enhancement factor for filter2 (default 1.0 : no correction) ; ; KEYWORDS: ; ; CAL - [Optional input] If set as "cal = 1", the coronal temperature is estimated based on the calibration in Narukage et al. (2011). ; If set as "cal = 2", the coronal temperature is estimated based on the calibration in Narukage et al. (2013). ; Default is "cal = 2" to use the latest calibration result of Narukage et al. (2013). ; BIN - [Optional input] (long) If set, data is binned using this value in spatial. ; After binning process, temperature is derived. ; TE_ERR_THRESHOLD - [Optional input] (float) You can adjust the threshold ratio of temperature error. (default = 0.1 [10%]) ; PHOTON_NOISE_THRESHOLD - [Optional input] (float) You can adjust the threshold ratio of photon noise to signal. (default = 0.1 [10%]) ; /NO_THRESHOLD - [Optional] (Boolean) If set, no threshold is set. ; /INFO - [Optional] (Boolean) If set, information is shown. ; LAMBDA_SPECTRUM - [Optional input] (float array) wavelength for SPECTRUM in an unit of [A]. ; TE_SPECTRUM - [Optional input] (float array) temperature for SPECTRUM in an unit of [log K]. ; SPECTRUM - [Optional input] (2-dim float array, [lambda, te]) photon number spectrum from solar plasma in an unit of [cm^3 s^-1 sr^-1]. ; We recommend that ; (1) The data point of wavelength is from 1 to 400 [A] with 0.1 [A] resolution. ; (2) The data point of temperature is from 10^5.0 to 10^8.0 [K] with 10^0.05 [K] resolution. ; /APAC_SPECTRUM - [Optional] (Boolean) If set, solar spectrum calculated with APAC database is used. (in default, solar spectrum is calculated with CHIANTI database.) ; EFF_AREA - [Optional] (Structure) You can input the output (effective area) from MAKE_XRT_WAVE_RESP.PRO. ; ; OUTPUTS: ; ; TE - (2-dim float array, [x, y]) derived temperature in log scale [log K]. ; EM - (2-dim float array, [x, y]) derived volume emission measure in log scale [log cm^-3]. ; ET - (2-dim float array, [x, y]) error of temperature in log scale [log K]. ; EE - (2-dim float array, [x, y]) error of volume emission measure in log scale [log cm^-3]. ; ; EXAMPLES: ; ; Using this function, you can derive the coronal temperature using filter ratio method. ; IDL> xrt_teem, index1, data1, index2, data2, te, em, et, ee ; ; If you want to bin data in spatial to collect photons (to reduce photon noise), set /bin keyword as following. ; In this case, data is binned as 3x3 in spatial at first. After this, temperature is derived with binned data. ; IDL> xrt_teem, index1, data1, index2, data2, te, em, et, ee, bin = 3 ; ; COMMON BLOCKS: ; ; none ; ; NOTES: ; ; The returned value of pixels where temperature cannot be derived or error is grater than threshold will be 0. ; ; If you input 3-dim XRT data [x, y, t], the data is summed along time axis and then the data becomes 2-dim array [x, y]. ; After this process, temperature is derived using summed data. ; ; The detail of coronal-temperature-diagnostic capability of the Hinode/XRT is described in ; Narukage et al. 2011, Solar Phys., 269, 169. ; http://adsabs.harvard.edu/doi/10.1007/s11207-010-9685-2 ; and ; Narukage et al. 2013, Solar Phys., in press ; http://adsabs.harvard.edu/doi/10.1007/s11207-013-0368-7 ; These two papers are the reference papers of this program. ; ; CONTACT: ; ; Comments, feedback, and bug reports regarding this routine may be ; directed to this email address: ; noriyuki.narukage ~at~ nao.ac.jp ; ; MODIFICATION HISTORY: ; progver = 'v2007-May-18' ;--- (N.Narukage (ISAS/JAXA)) Written. progver = 'v2009-Jul-27' ;--- (N.Narukage (NAOJ)) Updated to consider the contamination ; on focal-plane analysis filters and CCD. progver = 'v2010-Jul-30' ;--- (N.Narukage (NAOJ)) Updated to input the solar spectrum. ; and added the option to select the APAC database. progver = 'v2010-Aug-04' ;--- (N.Narukage (NAOJ)) Updated to input the output from MAKE_XRT_WAVE_RESP.PRO. progver = 'v2010-Aug-12' ;--- (N.Narukage (NAOJ)) Updated to return an error message ; for the inputs of DATA1 and DATA2 normalized by . progver = 'v2011-Feb-17' ;--- (N.Narukage (NAOJ)) Added the information on the reference paper. progver = 'v2013-Nov-04' ;--- (N.Narukage (NAOJ)) Updated to use the calibration result of Narukage et al. (2013). ; Added the keyword of "cal" to select the calibration result. ; And modified for the message from the program. ;--- (20-Oct-2014, Aki Takeda) Modified for the use for irradiance study. ;--- ( 2-Feb-2019, Aki Takeda) Added abundance keyword (corona, hybrid, photos). ; Also avoided pause when printing messages (hprint --> print). ;--- ( 8-Apr-2020, Aki Takeda) Renamed from xrt_teem4syn_ch901.pro to make clear that it accesses ; the solar spectrum files calculated with the latest isothermal.pro (ver.15 + local mod.) ;--- (22-Jun-2020, Aki Takeda) Modified for taking partially open pre-filter transmission. ;--- ( 6-May-2021, Aki Takeda) Modified for taking FxCTM_FACT keyword to consider increasing contamination on filters. ;--- ( 2-Jan-2022, Aki Takeda) Renamed for Chianti 10.0 version. ;--- (18-Aug-2023, Aki Takeda) Renamed for Chianti 10.1 version. /hybrid dropped. ;--- (20-Oct-2025, Aki Takeda) Renamed for Chianti 11.02 version. ;- ; ========================================================================= ; --- setting for message ------------------------------------------------- if not keyword_set(message) then message = ['XRT_TEEM'] ; --- summed data along time axis ----------------------------------------- ; n1 = total(strpos(index1.history, 'XRT_RENORMALIZE') ne -1) ; n2 = total(strpos(index2.history, 'XRT_RENORMALIZE') ne -1) ; if (n1 ne 0) or (n2 ne 0) then begin ; te = 0. ; em = 0. ; et = 0. ; ee = 0. ; if n1 ne 0 and n2 eq 0 then temp = 'DATA1 was' ; if n1 eq 0 and n2 ne 0 then temp = 'DATA2 was' ; if n1 ne 0 and n2 ne 0 then temp = 'DATA1 and DATA2 were' ; message = [message, ''] ; message = [message, '***** ERROR from XRT_TEEM.PRO *******************************************'] ; message = [message, 'DATA1 and DATA2 should be "non-normalized" data,'] ; message = [message, ' but your inputted '+temp+' normalized by XRT_PREP procedure.'] ; message = [message, 'Then XRT_TEEM procedure was stopped.'] ; message = [message, '******************************************* ERROR from XRT_TEEM.PRO *****'] ; return ; endif in1 = index1[0] in2 = index2[0] in1.EXPTIME = total(index1.EXPTIME) in2.EXPTIME = total(index2.EXPTIME) if n_elements(data1[0,0,*]) ge 2 then da1 = total(data1, 3) else da1 = data1 if n_elements(data2[0,0,*]) ge 2 then da2 = total(data2, 3) else da2 = data2 ; --- mask process -------------------------------------------------------- if not keyword_set(mask1) then mask1 = da1 * 0. if not keyword_set(mask2) then mask2 = da2 * 0. if n_elements(mask1[0,0,*]) ge 2 then mk1 = total(mask1, 3) ne 0 else mk1 = mask1 ne 0 if n_elements(mask2[0,0,*]) ge 2 then mk2 = total(mask2, 3) ne 0 else mk2 = mask2 ne 0 ; --- bin option process -------------------------------------------------- if keyword_set(bin) then begin s = size(da1) x = rebin(indgen(s[1]), s[1], s[2]) y = rotate(rebin(indgen(s[2]), s[2], s[1]),1) d1 = fltarr(s[1], s[2]) d2 = fltarr(s[1], s[2]) m1 = fltarr(s[1], s[2]) m2 = fltarr(s[1], s[2]) for i=0, bin - 1 do begin for j=0, bin - 1 do begin mx = fix(x/bin)*bin + i my = fix(y/bin)*bin + j d1 = d1 + da1(mx, my) d2 = d2 + da2(mx, my) m1 = m1 + mk1(mx, my) m2 = m2 + mk2(mx, my) endfor endfor da1 = d1 da2 = d2 mk1 = m1 ne 0 mk2 = m2 ne 0 endif else bin = 1. ; --- data check ---------------------------------------------------------- te = 'Unestimatable' em = 'Unestimatable' et = 'Unestimatable' ee = 'Unestimatable' ; --- set parameter ------------------------------------------------------- at_em = 0. ; --- make flux ----------------------------------------------------------- ;****** below added by AkT(20-Oct-2014) to use chianti ver.713 spectra ******* ;****** below added by AkT(21-Jul-2017) to use chianti ver.800 spectra ******* ;****** below added by AkT( 2-Feb-2019) to added abundance selection ******* ;****** below added by AkT(18-Mar-2020) to use chianti ver.901 spectra ******* ;****** below added by AkT(8-Apr-2020) to use chianti ver.0901 spectra ******* ;****** below added by AkT(2-Jan-2022) to use chianti ver.1000 spectra ******* ;****** below added by AkT(18-Aug-2023) to use chianti ver.1010 spectra ******* ;****** below added by AkT(20-Oct-2025) to use chianti ver.1102 spectra ******* dir_spec = '~takeda/idlpro/xrt/xrt_resp/solar_spectrum/' ;spfile = dir_spec+'solspec_ch713_corona_chianti' ;dir_spec = '$SSW/hinode/xrt/idl/response/solar_spectrum/' spfile = dir_spec+'solspec_ch1102_cr2021chianti_adv.genx' ; default is corona ; if keyword_set(hybrid) then spfile = dir_spec+'solspec_ch1000_hybrid_chianti' if keyword_set(photospheric) then spfile = dir_spec+'solspec_ch1102_ph2021asplund_adv.genx' restgen, te_spectrum, lambda_spectrum, spectrum, file = spfile ;****** above added by AkT(20-Oct-2014) to use chianti ver.713 spectra ******* ;****** above added by AkT(21-Jul-2017) to use chianti ver.800 spectra ******* ;****** above added by AkT( 2-Feb-2019) to added abundance selection ******* ;****** above added by AkT(18-Mar-2020) to use chianti ver.901 spectra ******* ;****** above added by AkT(8-Apr-2020) to use chianti ver.0901 spectra ******* ;****** above added by AkT(2-Jan-2022) to use chianti ver.1000 spectra ******* ;****** above added by AkT(18-Aug-2023) to use chianti ver.1010 spectra ******* ;****** above added by AkT(20-Oct-2025) to use chianti ver.1102 spectra ******* ; flux1 = xrt_flux1102_slmd_fctm(t, cal = cal, index = in1, vem = at_em, /datapoint, lambda_spectrum = lambda_spectrum, te_spectrum = te_spectrum, spectrum = spectrum, apac_spectrum = apac_spectrum, eff_area = eff_area, message = message, frac_open=frac_open, fctm_fact=f1ctm_fact) flux2 = xrt_flux1102_slmd_fctm(t, cal = cal, index = in2, vem = at_em, /datapoint, lambda_spectrum = lambda_spectrum, te_spectrum = te_spectrum, spectrum = spectrum, apac_spectrum = apac_spectrum, eff_area = eff_area, message = message, frac_open=frac_open, fctm_fact=f2ctm_fact) temp = min(abs(t-6.3), pos) temp = flux1 / flux2 if temp[pos] gt 1. then rev_ratio = 1 else rev_ratio = 0 ; check the ratio at 2MK ;stop ; --- trange process ------------------------------------------------------ if keyword_set(trange) then begin n = where(t ge min(trange) and t le max(trange)) if n[0] ne -1 then begin t = t[n] flux1 = flux1[n] flux2 = flux2[n] endif else begin message = [message, ''] message = [message, '***** WARNING from XRT_TEEM.PRO *****************************************'] message = [message, 'The inputted "TRANGE" keyword is out of range.'] message = [message, 'Then your inputted "TRANGE" is ignored.'] message = [message, '***************************************** WARNING from XRT_TEEM.PRO *****'] endelse endif ; --- make ratio ---------------------------------------------------------- data_ratio = ( da1 / in1.EXPTIME ) / ( da2 / in2.EXPTIME ) model_ratio = flux1 / flux2 if rev_ratio eq 1 then begin data_ratio = ( da2 / in2.EXPTIME ) / ( da1 / in1.EXPTIME ) model_ratio = flux2 / flux1 endif ; --- derive Te ------------------------------------------------------- s = size(data_ratio) ok_num = fltarr(s[1], s[2]) ok_cnt = fltarr(s[1], s[2]) for i=1, n_elements(model_ratio)-1 do begin n = where( data_ratio ge min(model_ratio[i-1:i]) and data_ratio le max(model_ratio[i-1:i]) ) if n[0] ne -1 then begin ok_num[n] = i ok_cnt[n] ++ endif endfor ok_num = ok_num * (ok_cnt eq 1) a = abs(model_ratio[ok_num] - data_ratio) b = abs(model_ratio[(ok_num-1)>0] - data_ratio) te = (t[ok_num] * b + t[(ok_num-1)>0] * a) / (a + b) OK_pixel = (ok_cnt eq 1) * (da1 ge 0) * (da2 ge 0) * (mk1 eq 0) * (mk2 eq 0) NG_pixel = where( OK_pixel eq 0 ) if NG_pixel[0] ne -1 then te[NG_pixel] = 0. ;stop ;wdef,10,500,500 ;plot,t,model_ratio ; ; --- derive EM ----------------------------------------------------------------- DN = INTERPOL(flux1, t, te) em = alog10( da1 / (DN * in1.EXPTIME) ) + at_em - alog10(bin^2.) if NG_pixel[0] ne -1 then em[NG_pixel] = 0. ; --- derive Te error, EM error ------------------------------------------------- dlnR_dlnT = abs( DERIV(alog(10.^t), alog(model_ratio)) ) dlnR_dlnT = INTERPOL(dlnR_dlnT, t, te) K1 = xrt_cvfact(temp, cal = cal, index = in1, /datapoint, /error, lambda_spectrum = lambda_spectrum, te_spectrum = te_spectrum, spectrum = spectrum, eff_area = eff_area, message = message) K1 = INTERPOL(K1, temp, te) K2 = xrt_cvfact(temp, cal = cal, index = in2, /datapoint, /error, lambda_spectrum = lambda_spectrum, te_spectrum = te_spectrum, spectrum = spectrum, eff_area = eff_area, message = message) K2 = INTERPOL(K2, temp, te) dlnf1_dlnT = DERIV(alog(10.^t), alog(flux1)) dlnf1_dlnT = INTERPOL(dlnf1_dlnT, t, te) dlnf2_dlnT = DERIV(alog(10.^t), alog(flux2)) dlnf2_dlnT = INTERPOL(dlnf2_dlnT, t, te) et = alog10( 1./dlnR_dlnT * sqrt(K1/da1 + K2/da2) ) + te if NG_pixel[0] ne -1 then et[NG_pixel] = 0. ee = alog10( 1./dlnR_dlnT * sqrt(dlnf2_dlnT^2.*K1/da1 + dlnf1_dlnT^2.*K2/da2) ) + em if NG_pixel[0] ne -1 then ee[NG_pixel] = 0. ;stop ; --- threshold process --------------------------------------------- if not(keyword_set(te_err_threshold)) then te_err_threshold = 0.5 ; changed 2020 if not(keyword_set(photon_noise_threshold)) then photon_noise_threshold = 0.2 ; changed 2020 if not(keyword_set(no_threshold)) then begin OK_pixel = OK_pixel * ( (et - te) le alog10(te_err_threshold) ) * ( sqrt(K1/(da1>0)) le photon_noise_threshold ) * ( sqrt(K2/(da2>0)) le photon_noise_threshold ) NG_pixel = where( OK_pixel eq 0 ) if NG_pixel[0] ne -1 then te[NG_pixel] = 0. if NG_pixel[0] ne -1 then em[NG_pixel] = 0. if NG_pixel[0] ne -1 then et[NG_pixel] = 0. if NG_pixel[0] ne -1 then ee[NG_pixel] = 0. message = [message, ''] message = [message, '----- NOTE from XRT_TEEM.PRO --------------------------------------------'] message = [message, '+ Examined Te range : '+ string(min(t), format='(f4.2)') + ' - ' + string(max(t), format='(f4.2)')+ ' [log K]'] message = [message, '+ Applied thresholds :'] message = [message, ' - Te error < ' + string(te_err_threshold*100. , format='(f6.2)') + '%'] message = [message, ' - Photon noise < ' + string(photon_noise_threshold*100., format='(f6.2)') + '%'] message = [message, '-------------------------------------------- NOTE from XRT_TEEM.PRO -----'] endif else begin message = [message, ''] message = [message, '----- NOTE from XRT_TEEM.PRO --------------------------------------------'] message = [message, '+ Examined Te range : '+ string(min(t), format='(f4.2)') + ' - ' + string(max(t), format='(f4.2)')+ ' [log K]'] message = [message, '-------------------------------------------- NOTE from XRT_TEEM.PRO -----'] endelse ; --- show message ------------------------------------------------------------ if n_elements(message) ge 2 then begin if message[0] eq 'XRT_TEEM' then print, [message[1:*], ''] ; if message[0] eq 'XRT_TEEM' then hprint, [message[1:*], ''] endif end