FUNCTION calc_xrt_temp_resp, specfile=specfile, trans_file=trans_file, $ electrons=electrons, ccdport=ccdport, $ qabort=qabort, qstop=qstop ; ========================================================================= ;+ ; PROJECT: ; Solar-B / XRT ; ; NAME: ; ; CALC_XRT_TEMP_RESP ; ; CATEGORY: ; ; Instrument response ; ; PURPOSE: ; ; Produce the temperature response for each XRT x-ray channel, ; assuming a spectral emission model. See Note #1. ; WARNING: This routine has been deprecated. See Note #2. ; ; CALLING SEQUENCE: ; ; Result = CALC_XRT_TEMP_RESP( [,specfile=specfile] [,/electrons] ; [,ccdport=ccdport] [,qabort=qabort] ; [,/qstop] ) ; ; INPUTS: ; ; No parameters. ; ; KEYWORDS: ; ; SPECFILE - [Optional] (string scalar) ; If SPECFILE is not given, a default file will be used. ; A "genx" type file that contains an XRT structure ; of type="spectrum". These can be generated with ; . (See that code for a ; detailed description of the contents.) The "spectrum" ; structure contains information about a spectral ; emission model as a function of wavelength and ; temperature. ; /ELECTRONS - [Optional] (Boolean) Cause the temperature response ; output to have units of 'el cm^5 s^-1 pix^-1'. ; The default is 'DN cm^5 s^-1 pix^-1'. ; CCDPORT - [Optional] (string scalar) The XRT flight camera ; has two read-out ports, which have slightly ; different gains (electrons per DN). This keyword ; may be used to select for the port, giving a string ; value of 'LEFT' or 'RIGHT'. For XRT data stored in ; an INDEX/DATA pair, you can find out the port with ; IDL> print, index.readport ; However, it is generally quite safe to rely on the ; default value, especially when representing the ; XRT temperature response directly (i.e., this keyword ; is generally not necessary). ; (Also see GET_XRT_GAIN.pro.) ; /QSTOP - [Optional] (Boolean) For debugging. ; ; OUTPUTS: ; ; Result - [Mandatory] (structure array, [Nchannels]) ; The temperature response for each XRT x-ray channel. ; This function is an integral of a spectral emission ; model with the channel's spectral response function. ; Units = 'DN cm^5 s^-1 pix^-1'. (Since a "column" ; emission measure has units of 'cm^-5', the units of ; the T-response may also be written as ; 'DN s^-1 pix^-1 EM^-1'.) See Note #1. ; Here is a description of the "temp_resp" structure. ; TYPE STRING 'temp_resp' ; CHANNEL_NAME STRING ; TEMP FLOAT Array[50] ; TEMP_UNITS STRING 'log K' ; TEMP_RESP FLOAT Array[50] ; TEMP_RESP_UNITS STRING 'DN cm^5 s^-1 pix^-1' ; LENGTH LONG ; DATA_FILES STRING Array[6] ; HISTORY STRING Array[5] ; COMMENTS STRING Array[5] ; TYPE = 'temp_resp', common to all XRT temp-response ; structures. ; LENGTH = N_temperatures (may be less than 50); also ; applies to TEMP_RESP. ; QABORT - [Optional] (Boolean) Indicates that the program ; exited gracefully without completing. (Might be ; useful for calling programs.) ; 0: Program ran to completion. ; 1: Program aborted before completion. ; ; EXAMPLES: ; ; Basic usage: ; IDL> t_resp = calc_xrt_temp_resp() ; IDL> help, t_resp ; IDL> help, t_resp, /st ; ; COMMON BLOCKS: ; ; none ; ; NOTES: ; ; 1) A discussion of units. ; It is easier to understand these units in the larger context. ; For a spectral model S, such as is returned by ; : S(wave, T) ~ [ph cm^3 s^-1 sr^-1 A^-1]. ; For a spectral response R, such as ; returns: ; R(wave) ~ [el cm^2 sr ph^-1 pix^-1], which says something about ; how many electrons are generated in the CCD for a photon of ; a given wavelength. Then: ; S(wave, T) * R(wave) ~ [el A^-1 s^-1 pix^-1 (cm^-5)^-1] ; ~ [el A^-1 s^-1 pix^-1 EM^-1], ; where EM ~ [cm^-5] is the "line of sight" (or "column") ; emission measure. For a channel's temperature response F{T): ; F(T) = integrate{ S(wave, T) * R(wave) * d(wave)} ; ~ [el s^-1 pix^-1 EM^-1] ; Note that one must assume a spectral model to calculate a ; temperature response. ; ; Going just a bit further, for an EM at T = T0, the signal G ; observed is ; G|T0 = F(T0) * EM|T0 ~ [el s^-1 pix^-1]. ; For an emission measure continuously distributed over a range ; of temperatures T, one uses the differential emission measure: ; DEM(T) ~ [cm^-5 K^-1]. ; Now the net signal G is ; G = integrate{F(T) * DEM(T) * dT} ~ [el s^-1 pix^-1]. ; ; 2) This routine has been deprecated. It is not compatible with the ; most recent calibrations, and does not account for the CCD ; contamination layer. Use instead. ; ; CONTACT: ; ; Comments, feedback, and bug reports regarding this routine may be ; directed to this email address: ; xrt_manager ~at~ head.cfa.harvard.edu ; ; MODIFICATION HISTORY: ; progver = 'v2007-May-17' ;--- (M.Weber) Written. (Based on heavily on ; earlier work, though.) progver = 'v2008-Oct-01' ;--- (M.Weber) This routine has been ; deprecated. See Note #2. ; ;- ; ========================================================================= ; === Initial setup ============================== ;=== Initialize program constants. prognam = 'CALC_XRT_TEMP_RESP' prognul = ' ' ;=== Set some keyword Booleans. qspecfile = keyword_set(specfile) qelec = keyword_set(electrons) qport = keyword_set(ccdport) qstop = keyword_set(qstop) qabort = 0B ; === Check inputs =============================== ;; SPECFILE is evaluated inside . ; === Get spec model and XRT response ============ ;=== Get spectral model, S(wavelength, temperature) spec = get_xrt_spec_genx(infile=specfile, qabort=qabort) if qabort then begin return, 0 endif ntmp = spec.tlength swave0 = spec.wave[0:spec.wlength-1] ;; dsw0 is the differential for the integral. dsw0 = (shift(swave0,-1) - shift(swave0,1)) / 2. dsw0[0] = (swave0[1] - swave0[0]) / 2. dsw0[spec.wlength-1] = (swave0[spec.wlength-1]-swave0[spec.wlength-2])/2. spec0 = spec.spec[0:spec.wlength-1,0:ntmp-1] ;=== Get spectral response, R(wavelength) if keyword_set(trans_file) then begin print, 'calc_tresp: ' + trans_file end sp_resp = calc_xrt_spec_resp(trans_file=trans_file, qabort=qabort) if qabort then begin return, 0 endif nchn = n_elements(sp_resp) ; === Get the camera gain ======================== ;=== If going to convert to DN... if not qelec then begin gain = -1d if qport then gain = get_xrt_gain(ccdport, qabort=qabort) if ((gain[0] eq -1) or (qport eq 0)) then begin gain = get_xrt_gain('default') qabort = 0B endif ;=== else if leaving in electrons... endif else gain = 1.0d ; === Generate temp_resp structure =============== outvar = {type:'temp_resp', channel_name:'', temp:fltarr(50), $ temp_units:'log K', temp_resp:fltarr(50), $ temp_resp_units:'', length:1L, data_files:strarr(6), $ history:strarr(5), comments:strarr(5) } outvar.temp[0:spec.tlength-1] = spec.temp[0:spec.tlength-1] if qelec then outvar.temp_resp_units = 'el cm^5 s^-1 pix^-1' $ else outvar.temp_resp_units = 'DN cm^5 s^-1 pix^-1' outvar.length = spec.tlength outvar.data_files[1:5] = spec.data_files[0:4] outvar.history = spec.history outvar.comments = spec.comments outvar = replicate(temporary(outvar), nchn) outvar.channel_name = sp_resp.channel_name outvar.data_files[0] = sp_resp.data_files[0] ; === Loop: integrate for temperature resp ======= for ichn = 0,(nchn-1) do begin rwave0 = sp_resp[ichn].wave[0:sp_resp[ichn].length-1] resp0 = sp_resp[ichn].spec_resp[0:sp_resp[ichn].length-1] ;=== Find swave0 within limits of rwave0. wss = where((swave0 ge min(rwave0)) and (swave0 le max(rwave0)), $ nwave) ;=== Interpolate the response function onto swave0, so it may ;=== may be integrated with the spectrum. resp1 = interpol(resp0, rwave0, swave0[wss], /quad) ;=== This is a fancy matrix multiplication that does the ;=== integration over wavelength, for every choice of temperature. outvar[ichn].temp_resp[0:ntmp-1] = $ (transpose(spec0[wss,*]) # (resp1 * dsw0[wss])) / gain endfor ; === Finish ===================================== if qstop then stop return, outvar END ;======================================================================