Blame view

src/idl/dustem_run_example.pro 4.97 KB
9ebef1db   Jean-Philippe Bernard   First commit
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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
PRO dustem_run_example,model=model,help=help,plot_it=plot_it,png=png,postcript=postcript

;+
; NAME:
;    dustem_run_example
; PURPOSE:
;    This is an example of how to run dustem SEDs with the dustem wrapper.
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_run_example,model=model,help=help,/postcript,/help,/plot_it
; INPUTS:
;    None
; OPTIONAL INPUT PARAMETERS:
;    None
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    model = Selects one of the dust mixture used by dustem
;           'MC10' from Compiegne et al 2010 (default)
;           'DBP90' from Desert et al 1990
;           'DL01' from Draine & Li 2001
;           'DL07' from Draine & Li 2007
;    postcript = if set, plot is done in DUSTEM/Docs/Figures/Last_dustem_fit.ps
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    None
; RESTRICTIONS:
;    The dustem fortran code must be installed
;    The dustem idl wrapper must be installed
; PROCEDURE:
;    None
; EXAMPLES
;    dustem_run_readme,model='DBP90'
; MODIFICATION HISTORY:
;    Written by J.P. Bernard April 1st 2011
;    see evolution details on the dustem cvs maintained at IRAP
;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-

IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_run_example'
  goto,the_end
ENDIF

loadct,13

IF keyword_set(model) THEN BEGIN
  use_model=strupcase(model)
ENDIF ELSE BEGIN
  use_model='MC10'     ;Default is last dustem model
ENDELSE

;== This just initializes 
dustem_init,model=use_model
!dustem_verbose=1
!dustem_show_plot=1

IF keyword_set(png) THEN BEGIN
  dir_png=!dustem_dat+'/Figures/'
  force_mkdir,dir_png
ENDIF
IF keyword_set(postcript) THEN BEGIN
  dir_ps=!dustem_dat+'/Figures/'
  force_mkdir,dir_ps
ENDIF

IF not keyword_set(dir) THEN BEGIN
;  CASE getenv('DUSTEM_WHICH') OF
  CASE !dustem_which OF
    'DESERT':dir='DESERT_POST_ISO_MORE/'
    'COMPIEGNE': dir='MC_DAT/'
    'VERSTRAETE': dir='d_3.5/'
    'VERSTRAETE': dir='LV_DAT/' 
    'WEB3p8': dir_in=!dustem_soft_dir ; I changed this 
    'RELEASE': dir_in=!dustem_soft_dir ; I changed this 
    ELSE: dir='les_DAT/'
  ENDCASE
ENDIF

;== The following structure contains the inputs to the model
st_model=dustem_read_all(dir_in)
;help,st_model,/str

;== The following line writes inputs to the Fortran
dustem_write_all,st_model,!dustem_dat
;== This runs the Fortran. The ouput structure contains the results
st=dustem_run()
;== Optionaly, this creates an SED to be overploted
;stop
;filters=['IRAS1','IRAS4']
;Nfilt=n_elements(filters)
;sed=dustem_initialize_sed(Nfilt)
;sed.filter=filters
;sed.wave=dustem_filter2wav(filters)
;sed.instru=dustem_filter2instru(filters)
;stop
;stt=dustem_set_data(sed=sed)

;== The following plots the resulting emission SED
win=0L
window,win & win=win+2
xtit=textoidl('\lambda (\mum)')
;=== same plot as in Compiegne et al. 2010
ytit=textoidl('\nuI_\nu^{em} (W/m^2/sr for N_H=1.e20 H/cm^2)')
tit='DUSTEMWrapper Emitted Intensity '+use_model
yr=[1e-11,6.e-7]
xr=[1,5e3]
dustem_plot_nuinu_em,st,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog,title=tit,xtit=xtit,ytit=ytit
IF keyword_set(postcript) THEN BEGIN
  file_ps=dir_ps+'Last_dustem_run_readme_nuinuem.ps'
  previous_device=!d.name
  set_plot,'PS'
  device,filename=file_ps
  dustem_plot_nuinu_em,st,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog,title=tit,xtit=xtit,ytit=ytit
  device,/close
  set_plot,previous_device
  message,'Wrote '+file_ps,/continue
ENDIF

;== The following plots the resulting extinction
window,win & win=win+1
yr=[0,2.5]  ;range of 1/lambda
xr=[1,10]   ;range of sigma
xtit=textoidl('1/\lambda (\mum^{-1})')
ytit=textoidl('\sigma_{ext} (1e-21 cm^2/H)')
tit='DUSTEMWrapper extinction '+use_model
dustem_plot_extinction,st,st_model,xr=xr,yr=yr,/xsty,/ysty,xtit=xtit,ytit=ytit,title=tit
IF keyword_set(postcript) THEN BEGIN
  file_ps=dir_ps+'Last_dustem_run_readme_extinction.ps'
  previous_device=!d.name
  set_plot,'PS'
  device,filename=file_ps
  dustem_plot_extinction,st,st_model,xr=xr,yr=yr,/xsty,/ysty,xtit=xtit,ytit=ytit,title=tit
  device,/close
  set_plot,previous_device
  message,'Wrote '+file_ps,/continue
ENDIF

IF keyword_set(png) THEN BEGIN
  win=10L
  window,win ;& win=win+2
  file_png=dir_png+'Last_dustem_run_readme_nuinuem.png'
  xtit=textoidl('\lambda (\mum)')
  ytit=textoidl('\nuI_\nu^{em} (W/m^2/sr for N_H=1.e20 H/cm^2)')
  tit='DUSTEMWrapper Emitted Intensity '+use_model
  yr=[1e-11,6.e-7]
  xr=[1,5e3]
  dustem_plot_nuinu_em,st,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog,title=tit,xtit=xtit,ytit=ytit
  write_png,file_png,tvrd(/true)
  message,'Wrote '+file_png,/info
  ;=====
  file_png=dir_png+'Last_dustem_run_readme_extinction.png'
  yr=[0,2.5]  ;range of 1/lambda
  xr=[1,10]   ;range of sigma
  xtit=textoidl('1/\lambda (\mum^{-1})')
  ytit=textoidl('\sigma_{ext} (1e-21 cm^2/H)')
  tit='DUSTEMWrapper extinction '+use_model
  dustem_plot_extinction,st,st_model,xr=xr,yr=yr,/xsty,/ysty,xtit=xtit,ytit=ytit,title=tit
  write_png,file_png ,tvrd(/true),/verbose
  message,'Wrote '+file_png,/info
ENDIF

the_end:

END