PRO xrt_teem_c710, index1, data1, index2, data2, te, em, et, ee, $ 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, $ ratio_arr=ratio_arr, $ ; output keyword for debugging hybrid=hybrid, photospheric=photospheric, corona=corona, $ ; abundance switch photon_noise1=photon_noise1, photon_noise2=photon_noise2, $ ; output photon noise cr2012=cr2012 ; new abundance switch ; ========================================================================= ;+ ; PROJECT: ; ; Solar-B / XRT ; ; NAME: ; ; XRT_TEEM ; ; 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. ; ; KEYWORDS: ; ; 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 ; This is the reference paper 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 = 'v2011-Jul-25' ;--- (A.Takeda (MSU)) Added the abundance switch to the XRT_FLUX call. progver = 'v2011-Sep-27' ;--- (A.Takeda (MSU)) Added the keywords (optional output) to extract ; the amount of photon noise. progver = 'v2013-Feb-20' ;--- (A.Takeda (MSU)) Modified to access to the spectrum files ; calculated with Chianti ver.7.1 (by using XRT_FLUX710.PRO). progver = 'v2013-Feb-25' ;--- (A.Takeda (MSU)) Added a new keyword, /CR2012 to access to ; a new spectrum file. ; ;- ; ========================================================================= ; --- 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 message = 'DATA1 was' if n1 eq 0 and n2 ne 0 then message = 'DATA2 was' if n1 ne 0 and n2 ne 0 then message = 'DATA1 and DATA2 were' print print, '***** ERROR from XRT_TEEM.PRO *******************************************' print, 'DATA1 and DATA2 should be "non-normalized" data,' print, ' but your inputted '+message+' normalized by XRT_PREP procedure.' print, 'Then XRT_TEEM procedure was stopped.' print, '******************************************* ERROR from XRT_TEEM.PRO *****' print 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 ----------------------------------------------------------- flux1 = xrt_flux710(t, index = in1, vem = at_em, /datapoint, lambda_spectrum = lambda_spectrum, $ te_spectrum = te_spectrum, spectrum = spectrum, apac_spectrum = apac_spectrum, $ eff_area = eff_area, hybrid=hybrid, photospheric=photospheric, corona=corona, $ cr2012=cr2012) flux2 = xrt_flux710(t, index = in2, vem = at_em, /datapoint, lambda_spectrum = lambda_spectrum, $ te_spectrum = te_spectrum, spectrum = spectrum, apac_spectrum = apac_spectrum, $ eff_area = eff_area, hybrid=hybrid, photospheric=photospheric, corona=corona, $ cr2012=cr2012) 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 ; --- 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 print print, '***** WARNING from XRT_TEEM.PRO *****************************************' print, 'The inputted "TRANGE" keyword is out of range.' print, 'Then your inputted "TRANGE" is ignored.' print, '***************************************** WARNING from XRT_TEEM.PRO *****' print 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 ratio_arr=data_ratio ; ; --- 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. ; --- 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, index = in1, /datapoint, /error, lambda_spectrum = lambda_spectrum, te_spectrum = te_spectrum, spectrum = spectrum, eff_area = eff_area) K1 = INTERPOL(K1, temp, te) K2 = xrt_cvfact(temp, index = in2, /datapoint, /error, lambda_spectrum = lambda_spectrum, te_spectrum = te_spectrum, spectrum = spectrum, eff_area = eff_area) K2 = INTERPOL(K2, temp, te) photon_noise1=sqrt(K2*da1) photon_noise2=sqrt(K2*da2) 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. ; --- threshold process --------------------------------------------- if not(keyword_set(te_err_threshold)) then te_err_threshold = 0.5 if not(keyword_set(photon_noise_threshold)) then photon_noise_threshold = 0.2 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. print print, '----- NOTE from XRT_TEEM.PRO --------------------------------------------' print, '+ Examined Te range : '+ string(min(t), format='(f4.2)') + ' - ' + string(max(t), format='(f4.2)')+ ' [log K]' print, '+ Applied thresholds :' print, ' - Te error < ' + string(te_err_threshold*100. , format='(f6.2)') + '%' print, ' - Photon noise < ' + string(photon_noise_threshold*100., format='(f6.2)') + '%' print, '-------------------------------------------- NOTE from XRT_TEEM.PRO -----' print endif else begin print print, '----- NOTE from XRT_TEEM.PRO --------------------------------------------' print, '+ Examined Te range : '+ string(min(t), format='(f4.2)') + ' - ' + string(max(t), format='(f4.2)')+ ' [log K]' print, '-------------------------------------------- NOTE from XRT_TEEM.PRO -----' print endelse end