;+ ; Name: SPEX_GET_DEM ; ; Purpose: Calculate DEM over the range 1.e6 - 1.e9 K (.086 - 86.2 keV) using parameters and form of DEM requested. Plots ; DEM vs temperature if requested. ; ; Calling sequence: dem = spex_get_dem(summ=summ, par=par, form=form, interval=interval, units=units, plot=plot, temp=temp) ; ; Method: User supplies either a set of parameters and DEM form, or an OSPEX summ structure (contains the results of ; fitting time intervals with one of the multi_therm functinos) and an interval number. The DEM is calculated and returned ; either in units of keV^-1 or K^-1 depending on units keyword, and plotted if requested. ; ; Input keywords: ; par - array of parameters appropriate for the photon model function using the form of the DEM you want (look ; in f_multi_therm_xxx or f_multi+_therm_xxx_defaults for a description of the parameters and good default values. ; form - form of the DEM: 'exp'. 'pow', 'gauss', 'pow_exp', or '2pow' ; summ - (if par and form not supplied), an ospex summ structure containing fitted results using one of the multi_therm funcs ; interval - if summ is provided, this is the # of the interval # used to get parameters, defatuls is 0 ; NOTE: provide EITHER par and form OR summ and interval keywords. ; units - 'K' or 'KEV' (default is KEV) units of DEM to return ; use_temp_minmax - if set, only plot temperatures between min and max defined by parameters a[1] and a[2]. Default=0. ; plot - if set, plot DEM vs temperature in K or keV ; _extra - any plot parameters passed through to plot and legend routines (like charsize, yrange...) ; ; Output keywords: ; temp - array of temperature values used in K or keV depending on units keyword ; leg_text - text of legend with values of DEM parameters ; ; Output - returns array of DEM values in keV-1 or K-1 and a plot if requested ; ; Examples: ; dem=spex_get_dem(par=[1., 0.087, 8.5, 6., 8.00, 1e0, 1.], form='2pow') ; dem=spex_get_dem(par=[1., 0.087, 8.5, 6., 8.00, 1e0, 1.], form='2pow', /plot, units='k', temp=temp) ; ; Written: Kim Tolbert, 25-Mar-2014 ; Modifications: ; 26-Jun-2014, Kim. Added use_temp_minmax keyword ; 15-Dec-2014, Kim. Added linear_gauss option ; ;- function spex_get_dem, par=par, form=form, summ=summ, interval=interval, units=units, use_temp_minmax=use_temp_minmax, plot=plot, $ temp=temp, leg_text=leg_text, _extra=_extra if ~(keyword_set(summ) or (keyword_set(par) and keyword_set(form))) then begin message,'You must pass in a summ structure or par and form keywords.', /cont print," e.g. spex_get_dem, par=[.005,1.,4.,3.,2.], form='exp'" return, -1 endif checkvar, interval, 0 checkvar, units, 'KEV' units = strupcase(units) checkvar, use_temp_minmax, 0 if keyword_set(summ) then begin func_comps = arr2str(summ.spex_summ_fit_function, '+') q = where (strpos(func_comps,'multi_therm') ne -1, count) if count eq 0 then begin message, /cont, 'No multi thermal function found in summary results.' return, -1 endif mth = func_comps[q[0]] form = '' if strpos(mth, 'exp') then form = 'exp' if strpos(mth, 'pow') then form = 'pow' if strpos(mth, 'gauss') then form = 'gauss' if strpos(mth, 'linear_gauss') then form = 'linear_gauss' if strpos(mth, 'pow_exp') then form = 'pow_exp' if strpos(mth, '2pow') then form = '2pow' if form eq '' then begin message, 'DEM dependence is not one of allowed forms: exp, pow, gauss, pow_exp', /cont return, -1 endif a = summ.spex_summ_params[*,interval] endif else begin a = par endelse em = a[0] * 1.d49 ; emission measure resolve_routine,'f_multi_therm', /either, /compile_full_file tmin = a[1] tmax = a[2] leg_text = '' format = '(g8.3)' case form of 'exp': begin func = 'f_mexp' alpha = a[3] leg_text = 'alpha = ' + trim(alpha,format) end 'gauss': begin func = 'f_gauss' alpha = a[3] t_gauss = a[4] leg_text = 'alpha = ' + trim(alpha,format) + ' t_gauss = ' + trim(t_gauss,format) end 'linear_gauss': begin func = 'f_linear_gauss' alpha = a[3] t_gauss = a[4] leg_text = 'alpha = ' + trim(alpha,format) + ' t_gauss = ' + trim(t_gauss,format) end 'pow': begin func = 'f_mpow' alpha = a[3] leg_text = 'alpha = ' + trim(alpha,format) end 'pow_exp': begin func = 'f_mpow_mexp' alpha = a[3] tpeak = a[4] leg_text = 'alpha = ' + trim(alpha,format) + ' tpeak = ' + trim(tpeak,format) + ' keV' end '2pow': begin func = 'f_2pow' alpha = a[3] beta = a[4] t0 = a[5] leg_text = 'alpha = ' + trim(alpha,format) + ' beta = ' + trim(beta,format) + ' t0 = ' + trim(t0,format) + ' keV' end endcase kperkev = 11.6048e6 ; 1 keV = 11.6 MK npts = 1001 logtk = 6. + findgen(1001) / ((npts-1)/3.) ; log T(K) from 6 to 9 K ;logtk = 6. + findgen(1001) / ((npts-1)/2.) ; log T(K) from 6 to 8 K tk = 10^logtk ; T in K tkev = tk / kperkev ; T in keV if use_temp_minmax then begin q = where(tkev ge tmin and tkev le tmax, count) if count le 1 then begin message,/cont, 'No temperatures between you min/max temperatures (a[1] and a[2] parameters). Aborting.' return, -1 endif tkev = tkev[q] tk = tkev * kperkev endif ; units of DEM are keV^(-1). If requested units of K, i.e. DEM in K^(-1), then divide by kperkev dem = em * call_function(func, tkev, alpha, t_gauss=t_gauss, tpeak=tpeak, t0=t0, beta=beta) if units eq 'K' then dem = dem / kperkev temp = units eq 'K' ? tk : tkev if keyword_set(plot) then begin if units eq 'K' then begin xtitle = 'Temperature (K)' ytitle = 'DEM (cm!u-3!n K!u-1!n)' endif else begin xtitle = 'Temperature (keV)' ytitle = 'DEM (cm!u-3!n keV!u-1!n)' endelse plot, temp, dem, /xlog, /ylog, xtitle=xtitle, ytitle=ytitle, _extra=_extra if leg_text ne '' then ssw_legend, leg_text, box=0, _extra=_extra endif return, dem end