function sxt_flux2, LogTe, filter, em=em, t1=t1, date=date, version=version, $ DN=DN, gain_ccd=gain_ccd, photons=photons, frac_open=frac_open, $ noverbose=noverbose, file_db=file_db, sre_file=sre_file ;+ ; NAME: ; sxt_flux ; PURPOSE: ; Specify Te + filter + EM ==> compute DN or photon flux. ; Specify Te + filter + DN ==> compute EM ; ; CALLING SEQUENCE: ; Flux = sxt_flux(LogTe,filter_number) ; DN/s for log10(EM) = 44 ; Flux = sxt_flux(LogTe,index) ; DN/get_expdur(index) for gt_filtb(index) ; Flux = sxt_flux(LogTe,index,/photons) ; Photon/get_expdur(index) ; EM = sxt_flux(LogTe,index,DN=DN) ; Supply DN and return Emission measure ; ; INPUTS: ; LogTe = log10(Temp) (May be a vector: min = 5.5, max = 8.0) ; filter_number= Filter B value (range is 1 to 6) or SXT index structure. ; Note: If FILTER_NUMBER is an SXT index structure, then this ; will be used for Filter B, t1 (unless overridden with T1=T1) ; and date (unless overridden with DATE=DATE or ; FRAC_OPEN=FRAC_OPEN keywords). ; ; Either LogTe or filter_number may be vectors. However, if both are vectors ; then their lengths must be the same. In this case each element of the ; output vector, Flux(i), of SXT_FLUX will computed for LogTe(i) and ; filter_number(i). ; ; OUTPUTS: ; This function returns DN or Photon fluxes. ; If DN=DN keyword is set, will return Emission Measures. ; ; OPTIONAL INPUT KEYWORDS: ; em = log10(Emission measure). Default=44 (implies ==> 10^44 cm-3) ; DN = SXT DN counts. If specified, sxt_flux returns the log10(EM) ; instead of flux and the /photon and EM= keywords are ignored. ; t1 = Exposure time in msec (default = 1000 msec). ; Overrides exposure time specified by FILTER (index). ; date = Time in any format (including structure). Used to ; determine entrance filter transmission. ; Overrides date specified by FILTER (index). ; frac_open = Fraction of open entrance filter. 0 = launch, 1 = no filter. ; Overrides DATE and FILTER values. ; gain_ccd = camera gain in e-/Dn. (Default is 100). ; photons = If set, return photon flux (rather than DN). ; noverbose = If set, suppress some of the informational messages ; file_db = If provided, read the line spectral data from the specified ; file instead of the default Mewe file. ; sre_file = Same as FILE_DB (new SXT_TEEM like keyword). ; ; OPTIONAL OUTPUT KEYWORDS: ; version = version number of SXT electron vs Te data base file ; ; RESTRICTIONS: ; DATE or FRAC_OPEN (if specified) must be length 1. If the date is ; obtained from FILTER (= SXT index), the first time in the structure ; will be used to determine the entrance filter fraction. ; ; If the dates of interest span the time of 13-nov-92 16:50 UT, then ; separate calls must be made for data before and after this date. ; ; MODIFICATION HISTORY: ; 5-jul-93, J. R. Lemen, Written. ; 21-jul-93, JRL Modified to allow date to come in in external format ; 22-dec-93, JRL, Fixed an IDL V3.1 related bug ; 25-jan-94, JRL, Call get_yo_dates to get date/value of entrance filter ; 31-aug-94, JRL, Make sure sre data base has equally spaced temperatures ; Use interp1d to do linear interpolation. ; 14-mar-96, JRL, Include Loren Acton's header modifications. ; 31-jul-97, JRL, Added file_db keyword ; 15-nov-2004, S.L.Freeland - isolate sxt_teem common block in own routine ; 2-jul-2010, Aki Takeda - Update default sre-file, ; Mewe --> Chianti (sre_ch601_corona_chianti.genx) ;- @sxt_teem_common ;common sxt_teem_db, db, text, header, version1, Filein, Filein_last ; ------------------------ ; HARD-WIRED Assumptions: ; ------------------------ @teem_hardwired ;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 ;interp=0 ; HARD-WIRED *** Linear interpolation ; if no parameters are present, assume information mode if n_params() lt 2 then begin doc_library,'sxt_flux' print,format="(75('-'))" fastdoc,'sxt_flux',/summ ; Give the full parameter list print,format="(75('-'))" return,-1 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(gain_ccd) gt 1 then begin print,' **** Error in SXT_Flux (1)*******************************' & tbeep print,' gain_ccd must be length of 1' help,gain_ccd print,' *********************************************************' return,-1 endif D_date = '31-aug-91' ; Def = Day after Launch date D_t1 = 1000. ; Def integration time = 1000 ms if n_elements(EM) eq 0 then J_EM = 44. $; Default Emission measure (cm^3) else J_EM = float(EM) ;---------------------------------------------------------------------------- ; **** 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. ; Determine if we need to read the file or not: if n_elements(db) eq 0 then qread = 1 else qread = 0 case 1 of keyword_set(file_db) : filein = file_db keyword_set(sre_file) : filein = sre_file else : begin dir_chianti_sre='/disk/hl1/ylegacy/ydb_beta/sxt_response/' filein = file_search(dir_chianti_sre+'sre_ch601_corona_chianti.genx') ; (with coronal abun. & chianti.ioneq ) end endcase if not keyword_set(Filein_last) then Filein_last = Filein if Filein_last ne Filein then qread = 1 ; Read the file: if qread then begin ; sre = electrons vs Te print,'Reading data base file = ',filein restgen, str=db,file=filein, text=text, header=header Filein_last = Filein ; -- Coded added 31-aug-94: Make sure Temp vector is regularly spaced --- std = stdev(db.Temp(1:*)-db.Temp,dT) if abs(std/dT) ge 1.e-5 then begin message,'Aligning SXT response function to equally-spaced Temp grid',/info Ntemp = n_elements(db.Temp) Temp = db.Temp(0) + findgen(Ntemp) * dT ; Use the average T step size elects = transpose(db.elects)*0. ; Set up the variables photons = transpose(db.photons) ; Transpose for convenience for i=0,n_elements(elects(0,*)) -1 do elects(0,i) = dspline(db.Temp,db.elects(i,*), Temp,interp=0) for i=0,n_elements(photons(0,*))-1 do photons(0,i) = dspline(db.Temp,db.photons(i,*),Temp,interp=0) db.Temp = Temp ; Reassign to the data structure db.elects = transpose(elects) db.photons = transpose(photons) endif 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 for filter ; ***** if (n_elements(LogTe) eq 0) or (n_elements(filter) eq 0) then begin print,' **** Error in SXT_Flux (2)*******************************' & tbeep print,'sxt_flux: *** LogTe or Filter is not defined ' help,LogTe,Filter print,' *********************************************************' return,-1 endif sz1=size(LogTe) ; check args 1 and 3 sz2=size(filter) if (sz2(sz2(0)+1) eq 8) then begin ; Yes = filter = index j_filter = gt_filtb(filter) j_date = anytim2ints(filter(0)) ; Get the first date (JRL Mod 20-jul-93) j_t1 = gt_expdur(filter) endif else j_filter = filter if (min(j_filter) lt 1) or (max(j_filter) gt 6) then begin print,' **** Error in SXT_Flux (3)*********************************' & tbeep print,' Valid range for filter is between 1 and 6 (inclusive)' print,' ***********************************************************' return,-1 endif ; ***** ; Set up t1 ; ***** if n_elements(t1) gt 0 then j_t1 = float(t1) ; Make sure it is floating if n_elements(j_t1) eq 0 then j_t1 = d_t1 ; Use the default value ; ***** ; The next section of code ensures that the lengths of LogTe, Filter, Em, T1 all match. ; ***** if n_elements(DN) eq 0 then em_mode = 0 else begin J_DN = DN & em_mode = 1 if not keyword_set(noverbose) then $ print,'sxt_flux: Returning EM values (instead of flux)' endelse J_Te = LogTe ; Save the original value of LogTe NTe = N_elements(J_Te) & arg_problem = 0 ; Assume no problem NFt = N_elements(J_filter) NT1 = N_elements(J_T1) Nxx = N_elements(J_em) ; Are we supposed to return EM (emission measure values)? if keyword_set(photons) then j_photons = photons else j_photons = 0 if em_mode then begin Nxx = N_elements(J_DN) J_EM = em_assumed j_photons = 0 endif aa = [NTe,NFt,NT1,Nxx] if n_elements(uniq(aa(sort(aa)))) gt 2 then arg_problem=1 else begin if n_elements(uniq(aa(sort(aa)))) eq 2 then begin MaxN = max(aa) MinN = min(aa) if MinN ne 1 then arg_problem=1 else begin if NTe ne MaxN then J_Te = replicate(J_Te(0), MaxN) if NFt ne MaxN then J_filter = replicate(J_filter(0), MaxN) if NT1 ne MaxN then J_T1 = replicate(J_T1(0), MaxN) if Nxx ne MaxN then $ if em_mode then J_DN = replicate(J_DN(0), MaxN) else $ J_EM = replicate(J_EM(0), MaxN) endelse endif endelse ; n_elements(uniq(aa(sort(aa)))) ... if arg_problem then begin print,' **** Error in SXT_Flux (4)*******************************' & tbeep print, ' Length of Filter (and if specified, t1, em)' print, ' must match length of LogTe ' print, ' -- OR Filter (and other specified parameters) must have a length of 1' if em_mode then help,te,filter,t1,DN else help,te,filter,t1,EM print,' *********************************************************' return,-1 endif flux = J_Te * 0. ; Set up the output array ; ***** ; Set up frac_open ==> use frac_open=, then date=, then gt_date(index) ; ***** if n_elements(frac_open) eq 0 then begin if n_elements(date) ne 0 then j_date = anytim2ints(date) ; (JRL Mod 20-jul-93) if n_elements(j_date) eq 0 then j_date = d_date if n_elements(j_date) gt 1 then begin print,' **** Warning in SXT_Flux (1)*******************************' & tbeep print,' Date must be a scalar' j_date = j_date(0) print,' Will use the first value = ',fmt_tim(j_date) endif ; Get the open filter area (Changed after 13-Nov-92) J_Frac_open = get_yo_dates(j_date,/entr,/value) if not keyword_set(noverbose) and (J_Frac_open gt 0.) then $ print,'sxt_flux: Using open frac of entrance filter (frac_open) = ',j_frac_open,format='(a,f7.4)' endif else begin j_frac_open = Frac_open ; Use the user-supplied value if n_elements(j_frac_open) gt 1 then begin print,' **** Warning in SXT_Flux (2)*******************************' & tbeep print,' Frac_open must be a scalar' j_frac_open = j_frac_open(0) print,' Will use the first value supplied = ',j_frac_open endif if (j_frac_open gt 1.0) or (j_frac_open lt 0.) then begin print,' **** Warning in SXT_Flux (3)*******************************' & tbeep print,' Frac_open must be: 0. <= Frac_open <= 1.' print,' Value supplied = ',j_frac_open print,' Setting frac_open to 0.' j_frac_open = 0. endif endelse ; ***** ; Set up the data base vectors ; ***** te_db = reform(db.Temp) ; temperature array in database filt1_col = label_id([J_filter-1])-1 ; colum number in database filt1_unq = filt1_col(uniq(filt1_col,sort(filt1_col))) ; Get unique cases ;---------------------------------------------------------------------------- ; **** Step 4: Compute the Fluxes ;---------------------------------------------------------------------------- Tmin = min(te_db) & Tmax = max(te_db) if not keyword_set(noverbose) then begin if (min(J_Te) lt Tmin) or (max(J_Te) gt Tmax) then begin print,' **** Warning in SXT_Flux (4)*******************************' & tbeep print,' Valid range of Te must be between:',Tmin,' and',Tmax,format='(2(a,f5.2))' print,' A Value of 0 is returned for out of range values' endif endif for i=0,n_elements(filt1_unq)-1 do begin ; Set up Launch values first: if keyword_set(j_photons) then begin i1_db = transpose(db.photons(filt1_unq(i),*)) ; photons (Launch) i2_db = transpose(db.photons(6+filt1_unq(i),*)) ; photons (No entrance) endif else begin i1_db = transpose(db.elects(filt1_unq(i),*)) ; electrons (Launch) i2_db = transpose(db.elects(6+filt1_unq(i),*)) ; electrons (No entrance) endelse ; Modify for change in entrance filter: ;;print,'Frac_open=',J_frac_open ;*** if J_Frac_open le 1.0 then i1_db = i1_db*(1-J_Frac_open) + J_Frac_open*i2_db if not keyword_set(J_photons) then i1_db = i1_db / gain_ccd(0) ii = where((filt1_col eq filt1_unq(i)) and (J_Te ge Tmin) and (J_Te le Tmax),npts) if npts gt 0 then Flux(ii) = interp1d(Te_db, i1_db, J_Te(ii))* $ (j_t1(ii)/1.)*10.^(J_em(ii)-em_assumed) ; Assumed time=1.msec endfor ; IF EM is requested, compute it now: if em_mode then begin ii = where(Flux gt 0.,npts) if npts gt 0. then Flux(ii) = alog10(DN(ii) / Flux(ii)) + em_assumed endif return,Flux end