fit_sed_example.pro
2.62 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
PRO fit_sed_example
;== Initialize DustemWrap
dustem_init
;== set DustemWrap behavior
!dustem_verbose=1
!dustem_show_plot=1
;== Read an SED
dir=!dustem_wrap_soft_dir+'/Data/SEDs/'
file=dir+'Gal_composite_spectrum.xcat'
spec=read_xcat(file,/silent)
;=== SET THE OBSERVATION STRUCTURE
;dustem_set_data,spec
result1 = dustem_set_data(sed=spec)
;== define which model parameters to vary
;=== pd is the parameter definition
pd=['(*!dustem_params).grain.G0'];, '(*!dustem_params).grain.grains(1).mdust_o_mh','(*!dustem_params).grain.grains(0).mdust_o_mh'];,'(*!dustem_params).grain.grains(2).mdust_o_mh','(*!dustem_params).grain.grains(3).mdust_o_mh','(*!dustem_params).grain.grains(4).mdust_o_mh'];,'dustem_create_continuum_2']
;pd=[$
; '(*!dustem_params).G0', $ ;G0 radiation field intensity
; '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
; '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
; '(*!dustem_params).grains(2).mdust_o_mh',$ ;amCBEx mass fraction
; '(*!dustem_params).grains(3).mdust_o_mh',$ ;amCBEx mass fraction
; '(*!dustem_params).grains(4).mdust_o_mh',$ ;aSil mass fraction
; 'dustem_create_continuum_2'] ;Intensity of NIR continuum
Npar=n_elements(pd)
;=== iv gives the corresponding initial values
iv = [1.0];, 7.8e-4, 7.8e-4];,1.65e-4,1.45e-3,7.8e-3];,0.001]
;=== The following indicate if the parameters are upper or lower limited
ulimed=replicate(1,Npar)
llimed=replicate(1,Npar)
;=== The following gives corresponding lower limits
llims=[0.1] ;,0.0,0.0]
ulims=[100.]
;== SET THE FITTED PARAMETERS
dustem_init_parinfo,pd,iv,up_limited=ulimed, lo_limited=llimed,up_limits=ulims,lo_limits=llims
;dustem_init_parinfo,pd,iv,up_limited=ulimed, $
; lo_limited=llimed,up_limits=ulims,lo_limits=llims
;== Run the SED fitting
tol=1.e-14 ;Tolerance used
Nitermax=20 ;maximum number of iteration
;This is the criterium which will stop the fit procedure
loadct,13
;This is to ajust plot range
!y.range=[1e-4,10]
res=dustem_mpfit_data(tol=tol,Nitermax=Nitermax)
;=== Plot best fit
yr=[1e-4,10] & xr=[1,2000]
tit='DUSTEM Example'
ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
xtit=textoidl('\lambda (\mum)')
errors=(*(*!dustem_fit).current_param_errors)
chi2=(*!dustem_fit).chi2
rchi2=(*!dustem_fit).rchi2
loadct,13
st=dustem_read_all_res(!dustem_res,/silent)
dustem_plot_nuinu_em,st,yr=[1e-11,6.e-7],/ysty,xr=[1,5e3],/xsty
;dustem_sed_plot,(*(*!dustem_fit).current_param_values), $
; ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty, $
; res=res,errors=errors,chi2=chi2,rchi2=rchi2
END