Viewing contents of file '../idllib/astron/contrib/varosi/vlibm/allpro/ir_spectrum_fit.pro'
;+
; NAME:
; IR_spectrum_fit
; PURPOSE:
; Fit an infrared emission spectrum with a simple 1D radiative transfer
; model of single temperature dust emission from a homogenous source
; and line of sight absorption by colder dust.
;
; CALLING:
; specfit = IR_spectrum_fit( wvdat, fdat )
;
; INPUTS:
; wvdat = array of wavelengths, in microns.
;
; fdat = array of fluxes at each wavelength,
; calibrated in Jansky's / arcsec^2.
;
; KEYWORDS:
; DUST_KAPPA = structure containing the dust absorption coefficients.
; FIT_PARM_STRUCT = structure containing the fit control parameters,
; default is silicate dust at Galactic Center distance, with
; even mix of graphite & silicate dust in Line-of-sight (LOS).
; SOURCE_NAME = string, optional, name of source.
; BEAM_AREA = optional area of beam in arcsec^2,
; and this area factor is applied to computed luminosities.
; /PLOT : plot the data and overplot fitting result,
; PLOT=2 prints hardcopy.
;
; OUTPUTS:
; Function returns a structure with tags containing
; spectrum fit results: Temperature, mass of dust, etc...
;
; OPTIONAL OUTPUTS:
; wvf = wavelength array.
; flux_obs = flux (at wvf) of fit to (wvdat,fdat).
; flux_src = flux of source with no extinction.
;
; COMMON BLOCKS:
; common IR_spectrum_fit, fitparms, wghts
; EXTERNAL CALLS:
; pro match
; pro non_Lin_Lsq (a version of IDL function curvefit)
; function finterpol
; function rad_trans_1D (called by non_Lin_Lsq)
; PROCEDURE:
; Setup common blocks for function rad_trans_1D,
; call pro non_Lin_Lsq to fit spectrum, put results into structure.
; HISTORY:
; Written: Frank Varosi, HSTX @ NASA/GSFC, 1996.
; FV.1998: added option to keep extinction fixed.
;-
function IR_spectrum_fit, wvdat, fdat, wvf, flux_obs, flux_src, $
PLOTIT=plotit, SHOW=show, INFO=info, PSYM=psym,$
DUST_KAPPA=Dust_Kappa, SOURCE_NAME=src, MAXITER=maxit, $
BEAM_AREA=beam_area, FIT_PARM_STRUCT=fitpario
common IR_spectrum_fit, fitparms, wghts
common IR_spectrum_comp, Fsil_LOS, Fsil_src1, Fsil_src2
common kappa_interp, wki, kigrf, kisil
common kappa_full, wkf, kgrf, ksil
common kappa_V, kvg, kvs
common wave_subs, j9, j12, j19, j20
common RT_1D_fixed, cLOS, Tbb
if N_struct( fitparms ) NE 1 then begin
fitparms = { fitparm10, Fsil_LOS: 0.7, $
Fsil_src: 1.0, $
Fsil_src2: 1.0, $
Threshold: 0.1, $
init_T: 250.0, $
init_T2: 0.0, $
init_mdust_src: 7e-7, $
init_mdust_src2: 0.0, $
init_mdust_LOS: 7e-4, $
constant_ext: 0, $
constant_T: 0, $
distance: 8.5e3, $
Max:1000.0, $
wave_min: 1.0, $
wave_max: 100.0, $
wave_fit_min: 4.0, $
max_iter:30 }
endif
if N_struct( fitpario ) EQ 1 then copy_struct, fitpario, fitparms $
else fitpario = fitparms
Fsil_src1 = ( fitparms.Fsil_src < 1) > 0
Fsil_src2 = ( fitparms.Fsil_src2 < 1) > 0
Fsil_LOS = ( fitparms.Fsil_LOS < 1) > 0
fitparms.threshold = fitparms.threshold > 0
if N_struct( wghts ) LE 0 then begin
wvw = { wvwt, wavelength:0.0, weight:1.0, threshold:0.0 }
wghts = replicate( wvw, 11 )
wghts.wavelength = [ 2.2, 3.2, 4.7, $
7.8, 8.7, 9.8, 10.3, 11.6, 12.4, 18, 20 ]
wghts(3).weight = 10
wghts(4).weight = 10
wghts(5).weight = 100
wghts(6).weight = 10
wghts(9).weight = 10
endif
Nwav = N_elements( wvdat )
if (Nwav LE 0) then return,0
if N_struct( Dust_Kappa ) EQ 1 then begin
wkf = Dust_Kappa.wavelen
m = min( abs( wkf - 0.55 ), jv )
kvg = Dust_Kappa.Kabs(jv,0) + Dust_Kappa.Ksca(jv,0)
kvs = Dust_Kappa.Kabs(jv,1) + Dust_Kappa.Ksca(jv,1)
w = where( wkf GE ( 1 < fitparms.wave_fit_min ) )
wkf = wkf(w)
kgrf = Dust_Kappa.Kabs(w,0) + Dust_Kappa.Ksca(w,0)
ksil = Dust_Kappa.Kabs(w,1) + Dust_Kappa.Ksca(w,1)
m = min( abs( wkf - 9.8 ), j9 )
m = min( abs( wkf - 12.4 ), j12 )
m = min( abs( wkf - 19 ), j19 )
m = min( abs( wkf - 20 ), j20 )
kigrf = exp( finterpol( aLog( kgrf ), aLog(wkf), aLog(wvdat) ) )
kisil = exp( finterpol( aLog( ksil ) ) )
wki = wvdat
endif
Ndat = N_elements( fdat )
if (Ndat LE 0) then return,0
if Ndat ne Nwav then begin
message,"size of wavelength and flux arrays not equal" $
+ string(7b),/INFO
print,Nwav,Ndat
return,0
endif
match, wvdat, wki, md, mki, COUNT=nmatch
if( nmatch LT Nwav ) then begin
message,"stored interpolate wavelengths do not match " $
+"those for data, must restore Dust_Kappa.idl" $
+ string(7b),/INFO
return,0
endif
if N_elements( beam_area ) ne 1 then beam_area = 1.0
LOS = fitparms.distance *!cv.pc
Jy = 1e-26 * 1e7 * 1e-4 ; (Watt/Jansky) * (erg/Watt) * (m/cm)^2
Lumfac = beam_area * 4*!PI * (LOS*Jy) * ( LOS /!cv.Lsun )
parms = [ fitparms.init_mdust_LOS, $
fitparms.init_T, $
fitparms.init_mdust_src ]
if (fitparms.init_T2 gt 0) and $
(fitparms.init_mdust_src2 gt 0) then $
parms = [ parms, fitparms.init_T2, fitparms.init_mdust_src2 ]
if fitparms.constant_ext then begin
cLOS = parms(0)
parmf = parms(1:*)
endif else if fitparms.constant_T then begin
cLOS = parms(0)
Tbb = parms(1)
parms = parms(0:2)
parmf = [parms(2)]
endif else parmf = parms
Npar = N_elements( parmf )
parmsav = parmf
match, wvdat, wghts.wavelength, mwd, mwvs,/ORIG,COUNT=nm
wfth = fltarr( Ndat )
wfth(mwd) = wghts(mwvs).threshold
if N_elements( maxit ) ne 1 then maxit = fitparms.max_iter
if !DEBUG then begin
help,Ndat,Npar,cLOS
print,wvdat,wfth,fdat
endif
wok = where( (fdat GT fitparms.threshold) AND (fdat GT wfth) AND $
(wvdat GT fitparms.wave_fit_min), nok )
if( nok GT Npar ) then begin
match, wvdat(wok), wghts.wavelength, mwd, mwvs,/ORIG,COUNT=nm
wvec = replicate( 1.0, nok )
if (nm GT 0) then wvec(mwd) = wghts(mwvs).weight
ntry = 0
if !DEBUG then print,ntry,parmf
nufnu = fdat(wok) ;; * (!cv.c * 1e4 * 1e-26)/wvdat(wok) * 1e10
;;wvec = 1/nufnu
FIT:
non_Lin_Lsq, wvdat(wok), nufnu, parmf, stdevp, chisq, ff, $
FUNC_NAME="rad_trans_1D", W=wvec, INFO=info, MAX=maxit
if !DEBUG then begin
ffit = Rad_Trans_1D( wvdat(wok), parmf )
print,ntry,parmf
print, wvdat(wok), wvec, nufnu, ffit
endif
if ( min( parmf ) LE 1e-29 ) AND (ntry LT 10) then begin
ntry = ntry+1
parmf = parmsav + parmsav * randomn( seed, Npar ) * 0.1
print,ntry,parmf
goto, FIT
endif
flux_obs = Rad_Trans_1D( wkf, parmf, /FULL, EXT=ext_LOS, /JY )
wvf = wkf
nuw = (!CV.c*1e4)/wkf
LIR_obs = -Lumfac * Trap_Int( nuw, flux_obs,/POW )
flux_src = flux_obs/(ext_LOS>1e-37)
LIR_src = -Lumfac * Trap_Int( nuw, flux_src,/POW )
if N_elements( parmf ) ge 4 then begin
flux_obs2 = Rad_Trans_1D( wkf, parmf, /FULL, /JY, /T2 )
flux_src2 = flux_obs2/(ext_LOS>1e-37)
LIR_src2 = -Lumfac * Trap_Int( nuw, flux_src2,/POW )
endif
endif else begin
parmf(*) = 0
stdevp = parmf
LIR_src = 0.0
LIR_obs = 0.0
chisq = -1.0
endelse
if fitparms.constant_ext then begin
parms(1:*) = parmf
stdevp = [ 0, stdevp ]
endif else if fitparms.constant_T then begin
parms(2) = parmf
stdevp = [ 0, 0, stdevp ]
endif else parms = parmf
Npar = N_elements( parms )
if N_elements( src ) NE 1 then src = ""
ksrc1 = Fsil_src1*ksil + (1-Fsil_src1)*kgrf
ksrc2 = Fsil_src2*ksil + (1-Fsil_src2)*kgrf
kLOS = Fsil_LOS*ksil + (1-Fsil_LOS)*kgrf
kLOSv = Fsil_LOS*kvs + (1-Fsil_LOS)*kvg
sfitp = create_struct( "Fsil_LOS", 0.0, $
"Fsil_src", 0.0, $
"Fsil_src2", 0.0, $
"distance", 0.0, $
"beam_area", beam_area, $
"ChiSq", chisq, $
"constant_ext", 0, $
"constant_T", 0, $
"Mdust_LOS", parms(0), $
"Tau_LOS_V", parms(0) * kLOSv, $
"Tau_LOS_9", parms(0) * kLOS(j9), $
"Tau_LOS_19", parms(0) * kLOS(j19), $
"LIR_obs", LIR_obs, $
"LIR_src_tot", LIR_src )
if Npar ge 5 then begin
sfitp = create_struct( sfitp, $
"LIR_src2", LIR_src2, $
"Mdust_src", parms(2), $
"Mdust_src2", parms(4), $
"Tau_src_20", parms(2) * ksrc1(j20), $
"Tau_src2_20", parms(4) * ksrc1(j20), $
"T", parms(1), $
"T2", parms(3) )
endif else begin
sfitp = create_struct( sfitp, $
"Mdust_src", parms(2), $
"Tau_src_20", parms(2) * ksrc2(j20), $
"T", parms(1) )
endelse
specfit = create_struct( NAME="specfit" + strtrim( Npar,2 ) + "_V28", $
sfitp, $
"Nfitp", nok, $
"parms", parms, $
"stdevp", stdevp, $
"source", src, $
"Rel_RA", 0.0, $
"Rel_DEC", 0.0, $
"Threshold", 0.0, "box_size", 0 )
copy_struct, fitparms, specfit
if keyword_set( plotit ) then plot_specfit, specfit, fitparms, wvdat,fdat, wkf,$
flux_obs, flux_src, flux_obs2, flux_src2, $
HARD=(plotit GT 1), SHOW=show, PSYM=psym
if !DEBUG GT 2 then stop
return, specfit
end