fit_test.pro
2.49 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_test
;;;;; On initialise sans polarisation. Il faut au préalable choisir un modèle (copier GRAIN_MC10.DAT par exemple dans GRAIN.DAT )
dustem_init,pol=0
rchi2_weight = { $
ext: 1., $
sed: 1.}
;;=== Set weight for reduced chi2 ponderation between emission, extinction, polarized extinction and emission ===
; Calculate total according to fileds in !dustem_data structure
total = 0.
for i = 0, n_tags(rchi2_weight)-1 do total = total + rchi2_weight.(i)
; values are relative to total
for i = 0, n_tags(rchi2_weight)-1 do rchi2_weight.(i) /= total
;stop
;== set DustemWrap behavior
!dustem_verbose=1
!dustem_show_plot=1
;;;;; On lit les fichiers de données (Emission et Extinction)
dir=!dustem_wrap_soft_dir+'/Data/SEDs/'
file=dir+'Gal_composite_spectrum.xcat'
spec=read_xcat(file,/silent)
file_ext=!dustem_wrap_soft_dir+'/Data/EXTs/Mathis1990_DISM.xcat'
extinction=read_xcat(file_ext,/silent)
;;;; on remplit les structures
result = dustem_set_data(sed=spec, ext=extinction)
;;; On définit les paramètres libres du modèle et les valeurs
;;; initiales
pd=['(*!dustem_params).grain.G0', '(*!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']
iv = [1.0, 3.E-04,5.E-04,1.4500E-03,2.E-03,0.001]
;=== paramètre fixé
;fpd=['(*!dustem_params).grain.grains(1).mdust_o_mh']
;fiv=[3.e-4];,1.e-16]
;dustem_init_fixed_params,fpd,fiv
;;; On définit les valeurs min max possibles des paramètres : ici
;;; pas de valeurs max
Npar=n_elements(pd)
ulimed=[0,0,0,0,0,0] ;pas de valeurs max
llimed=[1,1,1,1,1,0] ; on impose une valeur min pour 4 des 5 paramètres
llims=[0.1,1.e-6,1.e-6,1.e-6,1.e-6,0]
;ulims=[]
dustem_init_parinfo,pd,iv,up_limited=ulimed, lo_limited=llimed,up_limits=ulims,lo_limits=llims
tol=1.e-14 ;Tolerance used
Nitermax=20 ;maximum number of iteration
loadct,13
;This is to ajust plot range
!y.range=[1e-4,10]
res=dustem_mpfit_data(tol=tol,Nitermax=Nitermax)
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
;;;; On lit les fichiers de sortie .RES pour affichier la SED obtenue
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
END