OV-GSO-DC / dustem-wrapper_idl

Commit e93c23143ca70b7420f9ca46bebee8cf816cadf0

Authored by Annie Hughes
1 parent c598f318
Exists in master

cleaned up for tuto

Showing 2 changed files with 230 additions and 206 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_run_example.pro
  Diff comments   View file @e93c231
1   -PRO dustem_run_example,model=model $
2   - ,png=png $
3   - ,postscript=postscript $
4   - ,help=help
  1 +PRO dustem_run_example, model $
  2 + ,show=show $
  3 + ,postscript=postscript $
  4 + ,help=help
5 5  
6 6 ;+
7 7 ; NAME:
... ... @@ -11,17 +11,9 @@ PRO dustem_run_example,model=model $
11 11 ; CATEGORY:
12 12 ; DustEMWrap, Distributed, High-Level, User Example
13 13 ; CALLING SEQUENCE:
14   -; dustem_run_example,model=model,postscript=postscript,png=png,help=help
  14 +; dustem_run_example,model,show=show,postscript=postscript,help=help
15 15 ; INPUTS:
16   -; None
17   -; OPTIONAL INPUT PARAMETERS:
18   -; None
19   -; OUTPUTS:
20   -; None
21   -; OPTIONAL OUTPUT PARAMETERS:
22   -; None
23   -; ACCEPTED KEY-WORDS:
24   -; model = specifies the interstellar dust mixture used by DustEM
  16 +; model = specifies the interstellar dust mixture used by DustEM
25 17 ; 'MC10' model from Compiegne et al 2010 (default)
26 18 ; 'DBP90' model from Desert et al 1990
27 19 ; 'DL01' model from Draine & Li 2001
... ... @@ -33,11 +25,29 @@ PRO dustem_run_example,model=model $
33 25 ; polarisation. See Tables 2 and 3 of that paper for details.
34 26 ; 'G17_ModelB' model B from Guillet et al (2018)
35 27 ; 'G17_ModelC' model C from Guillet et al (2018)
36   -; 'G17_ModelD' model A from Guillet et al (2018)
37   -; postscript = if set, final plot is saved as postscript in the
38   -; current working directory
39   -; png = if set, final plot is saved as png in the
40   -; current working directory
  28 +; OPTIONAL INPUT PARAMETERS:
  29 +; show = vector specifying the type of plot(s) to show. Possible
  30 +; values (for all dust models) are
  31 +; ["emis", "extuv", "extir", "alb", "sdist","rfield"]
  32 +; Additionally, for the G17_MODEL?s with polarisation, you can specify
  33 +; [ "polext", "polsed", "align"]
  34 +; "emis" -- Stokes I SED in emission predicted by the model
  35 +; "extuv" -- Extinction in the UV predicted by the model
  36 +; "extir" -- Extinction in the IR predicted by the model
  37 +; "alb" -- Albedo of the grain types
  38 +; "sdist" -- Size distribution of the grain types
  39 +; "rfield" -- Dust-heating radiation field
  40 +; "polsed" -- Polarised intensity SED predicted by the model
  41 +; "polext" -- Polarised extinction predicted by the model
  42 +; "align" -- Grain alignment fraction predicted by the model
  43 +; "all" -- All plots appropriate to the model
  44 +; OUTPUTS:
  45 +; Plots
  46 +; OPTIONAL OUTPUT PARAMETERS:
  47 +; None
  48 +; ACCEPTED KEY-WORDS:
  49 +; postscript = on/off if set, plots are saved as postscript in the
  50 +; current working directory (and not shown on screen)
41 51 ; help = if set, print this help
42 52 ; COMMON BLOCKS:
43 53 ; None
... ... @@ -49,7 +59,7 @@ PRO dustem_run_example,model=model $
49 59 ; PROCEDURE:
50 60 ; None
51 61 ; EXAMPLES
52   -; dustem_run_example,model='DBP90',png='test_dbp90.png'
  62 +; dustem_run_example,'DBP90',show="all",/post
53 63 ; MODIFICATION HISTORY:
54 64 ; Written JPB Apr-2011
55 65 ; Evolution details on the DustEMWrap gitlab.
... ... @@ -61,135 +71,128 @@ IF keyword_set(help) THEN BEGIN
61 71 goto,the_end
62 72 ENDIF
63 73  
64   -loadct,13
65   -
66   -IF keyword_set(model) THEN BEGIN
67   - use_model=strupcase(model)
  74 +IF keyword_set(show) THEN BEGIN
  75 + use_show=show
68 76 ENDIF ELSE BEGIN
69   - use_model='MC10' ;Default is the Compiegne et al 2010 model
  77 + use_show=['all'] ;Default is to show everything
70 78 ENDELSE
71 79  
  80 +known_mdls=['MC10','DBP90','DL01','WD01_RV5P5B','DL07','J13','G17_MODELA','G17_MODELB','G17_MODELC','G17_MODELD']
  81 +pol_mdls=['G17_MODELA','G17_MODELB','G17_MODELC','G17_MODELD']
  82 +test_model = where(known_mdls eq strupcase(model),ct)
  83 +if ct eq 0 then begin
  84 + message,'ISM dust model '+model+' unknown',/continue
  85 + message,'Known models are MC10,DBP90,DL01,WD01_RV5P5B,DL07,J13,G17_MODELA,G17_MODELB,G17_MODELC,G17_MODELD',/continue
  86 + stop
  87 +end
  88 +test_pol_model = where(pol_mdls eq model,polct)
  89 +
72 90 ;== INITIALISE DUSTEM
73   -dustem_init,model=use_model
  91 +dustem_init,model=model
74 92 !dustem_verbose=1
75 93 !dustem_show_plot=1
76   -
77   -CASE !dustem_which OF
78   - 'DESERT': begin
79   - dir='DESERT_POST_ISO_MORE/'
80   - message,'DESERT version disabled for this routine',/info
81   - end
82   - 'COMPIEGNE': begin
83   - dir='MC_DAT/'
84   - message,'COMPIEGNE version disabled for this routine',/info
85   - end
86   - 'VERSTRAETE': begin
87   - dir='LV_DAT/'
88   - message,'VERSTRAETE version disabled for this routine',/info
89   - end
90   - 'RELEASE': dir_in=!dustem_soft_dir
91   - ELSE: dir_in=!dustem_soft_dir
92   - ;dir='les_DAT/'
93   -ENDCASE
  94 +!dustem_which='RELEASE'
  95 +dir_in=!dustem_soft_dir
94 96  
95 97 ;== Fill the following structure contains with default inputs to the model
96 98 st_model=dustem_read_all(dir_in)
97   -;help,st_model,/str
98   -
99   -;=== You can modify the model inputs by changing the values in the
100   -;=== st_model structure
101   -
102   -;=== e.g. modify the default value of the ISRF intensity used by DustEM
103   -;st_model.G0=2.
104 99  
105   -;===== e.g. modify keywords related to spinning dust. Does not work with DBP90 !
106   -;st_model.grains[0].type_keywords='logn-chrg-spin'
107   -;st_model.gas.NH=1.e2
108   -;st_model.gas.Tgas=1.e5
  100 +;=== You could modify the model inputs here by directly editing the values in the
  101 +;=== st_model structure. To inspect contents, type IDL> help,st_model,/str
  102 +;=== e.g.
  103 +;=== st_model.G0=2.0 ; change the default value of the ISRF intensity used by DustEM
  104 +;=== st_model.isrf.isrf=sqrt(st_model.isrf.isrf) ; change the shape of ISRF
  105 +;=== st_model.grains[0].mdust_O_mh=5.e-4 ; change the abundance of grain[0]
  106 +;=== st_model.grains[1].type_keywords='plaw' ; change the size distribution description of grain[1]
109 107  
110   -;== The following line saves the modified inputs for use by the Fortran
  108 +;== The following line saves the modified inputs for use by the fortran
111 109 dustem_write_all,st_model,!dustem_dat
112 110  
113 111 ;== The following line runs the Fortran. The ouput structure st contains the results
114 112 st=dustem_run()
115 113  
116   -;== Optionaly, this creates an SED to be overploted
117   -;stop
118   -;filters=['IRAS1','IRAS4']
119   -;Nfilt=n_elements(filters)
120   -;sed=dustem_initialize_sed(Nfilt)
121   -;sed.filter=filters
122   -;sed.wave=dustem_filter2wav(filters)
123   -;sed.instru=dustem_filter2instru(filters)
124   -;stop
125   -;stt=dustem_set_data(sed=sed)
126   -
127   -;== The following lines plot the resulting SED (emission)
128   -win=0L
129   -window,win & win=win+2
130   -xtit=textoidl('\lambda (\mum)')
131   -;=== same plot as in Compiegne et al. 2010
132   -ytit=textoidl('\nuI_\nu^{em} (W/m^2/sr for N_H=1.e20 H/cm^2)')
133   -tit='DustEMWrap Emitted Intensity '+use_model
134   -yr=[1e-11,6.e-7]
135   -xr=[1,5e3]
136   -dustem_plot_nuinu_em,st,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog,title=tit,xtit=xtit,ytit=ytit
137   -IF keyword_set(postscript) THEN BEGIN
138   - dir_ps='./'
139   - file_ps=dir_ps+postscript
140   - previous_device=!d.name
141   - set_plot,'PS'
142   - device,filename=file_ps
143   - dustem_plot_nuinu_em,st,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog,title=tit,xtit=xtit,ytit=ytit
144   - device,/close
145   - set_plot,previous_device
146   - message,'Wrote '+file_ps,/continue
147   -ENDIF
148   -
149   -;== The following plots the resulting SED (extinction)
150   -window,win & win=win+1
151   -yr=[0,2.5] ;range of 1/lambda
152   -xr=[1,10] ;range of sigma
153   -xtit=textoidl('1/\lambda (\mum^{-1})')
154   -ytit=textoidl('\sigma_{ext} (1e-21 cm^2/H)')
155   -tit='DustEMWrap extinction '+use_model
156   -dustem_plot_extinction,st,st_model,xr=xr,yr=yr,/xsty,/ysty,xtit=xtit,ytit=ytit,title=tit
157   -IF keyword_set(postscript) THEN BEGIN
158   - file_ps=dir_ps+file_basename(file_ps,'.ps')+'_extinction.ps'
159   - previous_device=!d.name
160   - set_plot,'PS'
161   - device,filename=file_ps
162   - dustem_plot_extinction,st,st_model,xr=xr,yr=yr,/xsty,/ysty,xtit=xtit,ytit=ytit,title=tit
163   - device,/close
164   - set_plot,previous_device
165   - message,'Wrote '+file_ps,/continue
166   -ENDIF
167   -
168   -IF keyword_set(png) THEN BEGIN
169   - win=10L
170   - window,win ;& win=win+2
171   - dir_png='./'
172   - file_png=dir_png+png
173   - xtit=textoidl('\lambda (\mum)')
174   - ytit=textoidl('\nuI_\nu^{em} (W/m^2/sr for N_H=1.e20 H/cm^2)')
175   - tit='DustEMWrap Emitted Intensity '+use_model
176   - yr=[1e-11,6.e-7]
177   - xr=[1,5e3]
178   - dustem_plot_nuinu_em,st,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog,title=tit,xtit=xtit,ytit=ytit
179   - write_png,file_png,tvrd(/true)
180   - message,'Wrote '+file_png,/info
181   - ;=====
182   - file_png=dir_png+file_basename(file_png,'.png')+'_extinction.png'
183   -
184   - yr=[0,2.5] ;range of 1/lambda
185   - xr=[1,10] ;range of sigma
186   - xtit=textoidl('1/\lambda (\mum^{-1})')
187   - ytit=textoidl('\sigma_{ext} (1e-21 cm^2/H)')
188   - tit='DustEMWrap extinction '+use_model
189   - dustem_plot_extinction,st,st_model,xr=xr,yr=yr,/xsty,/ysty,xtit=xtit,ytit=ytit,title=tit
190   - write_png,file_png ,tvrd(/true),/verbose
191   - message,'Wrote '+file_png,/info
192   -ENDIF
  114 +known_plots=["emis", "extuv", "extir", "alb", "sdist"]
  115 +if polct gt 0 then known_plots=["emis", "extuv", "extir", "alb", "sdist","polext", "polsed", "align"]
  116 +if strupcase(use_show) eq "ALL" then use_show=[known_plots,"rfield"]
  117 +
  118 +;== The following lines generate plots that show the outputs of the fortran.
  119 +;=== For this release (V2.0), we use the old method (dustem_show_fortran) to
  120 +;=== access model quantities.
  121 +;=== This method will be deprecated in a future release and plots will be
  122 +;=== constructed directly from the st and st_model structures.
  123 +match,known_plots,use_show,a,b
  124 +if total(a,/nan) ne -1 and keyword_set(postscript) then $
  125 + dustem_show_fortran,model=model,st=st,show=known_plots[a],hard=1,tit="DustEM_run_example "+model+" :"
  126 +if total(a,/nan) ne -1 and not keyword_set(postscript) then $
  127 + dustem_show_fortran,model=model,st=st,show=known_plots[a],tit="DustEM_run_example "+model+" :"
  128 +
  129 +;== The following lines generate a plot of the ISRF using the
  130 +;== st_model structure. For illustration purposes.
  131 +use_win=12
  132 +rfield_plot=["rfield"]
  133 +match,rfield_plot,use_show,a,b
  134 +if total(a,/nan) ne -1 then begin
  135 + xtit=textoidl('\lambda (\mum)')
  136 + ytit=textoidl('log ISRF [4 \pi I_\nu (erg/cm^2/s/Hz)]')
  137 + tit='DustEM_run_example ISRF'
  138 + yr=[min(st_model.isrf.isrf),5.*max(st_model.isrf.isrf)]
  139 + xr=[min(st_model.isrf.lambisrf),max(st_model.isrf.lambisrf)]
  140 + IF not keyword_set(postscript) THEN BEGIN
  141 + window,use_win,xs=600,ys=400,tit='DUSTEM_SHOW_FORTRAN: ISRF '+strtrim(use_win,2) & use_win=use_win+2
  142 + cgplot,st_model.isrf.lambisrf,st_model.isrf.isrf $
  143 + ,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog $
  144 + ,title=tit,xtit=xtit,ytit=ytit,color=cgcolor('black'),/nodata
  145 + oplot,st_model.isrf.lambisrf,st_model.isrf.isrf,col=cgcolor('black')
  146 + END ELSE BEGIN
  147 + file_ps="isrf.ps"
  148 + previous_device=!d.name
  149 + !y.thick = 5
  150 + !x.thick = 5
  151 + !p.thick = 5
  152 + !p.charsize = 1.3
  153 + !p.charthick = 5
  154 + !x.ticklen = 0.04
  155 + set_plot,'PS'
  156 + device,filename=file_ps,/portrait,/color
  157 + cgplot,st_model.isrf.lambisrf,st_model.isrf.isrf $
  158 + ,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog $
  159 + ,title=tit,xtit=xtit,ytit=ytit,color=cgcolor('black'),/nodata
  160 + oplot,st_model.isrf.lambisrf,st_model.isrf.isrf ,col=cgcolor('black')
  161 + device,/close
  162 + !y.thick = 0
  163 + !x.thick = 0
  164 + !p.thick = 0
  165 + !p.charsize = 1
  166 + !p.charthick = 0
  167 + !x.ticklen = 0.02
  168 + set_plot,previous_device
  169 + END
  170 +end
  171 +
  172 +goto,the_end
  173 +
  174 +;;== The following lines plot the resulting SED (emission)
  175 +;;== using st/st_model structure and a dedicated plotting function
  176 +;;== Prototype for future releases.
  177 +;; win=0L
  178 +;; window,win & win=win+2
  179 +;; xtit=textoidl('\lambda (\mum)')
  180 +;; ytit=textoidl('\nuI_\nu^{em} (W/m^2/sr for N_H=1.e20 H/cm^2)')
  181 +;; tit='DustEMWrap Emission: '+model
  182 +;; yr=[1e-11,6.e-7]
  183 +;; xr=[1,5e3]
  184 +;; dustem_plot_nuinu_em,st,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog,title=tit,xtit=xtit,ytit=ytit,post=postscript
  185 +
  186 +;;== The following plots the resulting SED (extinction)
  187 +;;== using st/st_model structure and a dedicated plotting function
  188 +;;== Prototype for future releases.
  189 +;; window,win & win=win+1
  190 +;; yr=[0,2.5] ;range of 1/lambda
  191 +;; xr=[1,10] ;range of sigma
  192 +;; xtit=textoidl('1/\lambda (\mum^{-1})')
  193 +;; ytit=textoidl('\sigma_{ext} (1e-21 cm^2/H)')
  194 +;; tit='DustEMWrap extinction: '+model
  195 +;; dustem_plot_extinction,st,st_model,xr=xr,yr=yr,/xsty,/ysty,xtit=xtit,ytit=ytit,title=tit,post=postscript
193 196  
194 197 the_end:
195 198  
... ...
src/idl/dustem_show_fortran.pro
  Diff comments   View file @e93c231
1 1 PRO dustem_show_fortran, model=model $
  2 + , st=st $
  3 + , restart=restart $
2 4 , data_dir=data_dir $
3 5 , dustem_dir=dustem_dir $
4 6 , nh=nh, fsed=fsed, sw=sw $
... ... @@ -58,11 +60,10 @@ PRO dustem_show_fortran, model=model $
58 60 ; ''emis'', ''extuv'', ''extir'', ''alb'', ''sdist'', ''polext'', ''polsed'', ''align''
59 61 ; Default is [''emis'']. SHOW=0 no plot
60 62 ; WN : array of window numbers for plots
61   -; HARD : array (as show) or keyword (/HARD) for hardcopy or filename of ps files, emission'
62   -; extinction albedo and size dist.
63   -; Default = [''show.ps'']
64   -;
65   -; OPTIONAL INPUT PARAMETERS:
  63 +; HARD : /HARD for postscript of requested plots (if set, plots will not
  64 +; show on screen)
  65 +;
  66 + ; OPTIONAL INPUT PARAMETERS:
66 67 ; None
67 68 ; OUTPUTS:
68 69 ; OPTIONAL OUTPUT PARAMETERS:
... ... @@ -119,11 +120,12 @@ ENDIF ELSE BEGIN
119 120 ENDELSE
120 121  
121 122 ; run the fortran via the wrapper once
122   -dustem_init,model=use_model
123   -st_model=dustem_read_all(!dustem_soft_dir)
124   -dustem_write_all,st_model,!dustem_dat
125   -st=dustem_run()
126   -
  123 +if keyword_set(restart) then begin
  124 + dustem_init,model=use_model
  125 + st_model=dustem_read_all(!dustem_soft_dir)
  126 + dustem_write_all,st_model,!dustem_dat
  127 + st=dustem_run()
  128 +endif
127 129  
128 130 data_path = !dustem_wrap_soft_dir+'/Data/'
129 131 fortran_path = !dustem_dat
... ... @@ -347,9 +349,10 @@ if keyword_set(data_dir) then data_path = data_dir
347 349 !x.ticklen = 0.04
348 350 set_plot,'ps'
349 351 ; device,file=hard(0),xs=26, ys=20,/portrait,/color
350   - device,file=hard(0),/portrait,/color
  352 + device,file="emis.ps",/portrait,/color
351 353 ihard = ihard + 1
352   - endif else begin
  354 + endif else begin
  355 + cleanplot
353 356 window, wn(0), xs=900,ys=600,tit='DUSTEM_SHOW_FORTRAN: SED '+strtrim(wn(0),2)
354 357 endelse
355 358  
... ... @@ -357,8 +360,8 @@ if keyword_set(data_dir) then data_path = data_dir
357 360 ytit = '!4m!3 I!d!4m!3!n (erg s!u-1!n cm!u-2!n sr!u-1!n)'
358 361 fine = 1
359 362 ; plot_oo, INDGEN(1),/NODAT,/xs,/ys,XR=xr,YR=yr,XTIT=xtit,YTIT= ytit,tit=tit
360   - plot_oo, INDGEN(1),/NODAT,xs=9,/ys,XR=xr,YR=yr,XTIT=xtit,YTIT= ytit
361   - axis, /data, xr=1d-5*clight/xr, xax=1,/xlo,xs=1,xtit='Frequency (GHz)'
  363 + cgplot, INDGEN(1),/NODAT,/xs,/ys,/xlo,/ylo,XR=xr,YR=yr,XTIT=xtit,YTIT= ytit,tit=tit+" SED",color=cgcolor('black')
  364 +; axis, /data, xr=1d-5*clight/xr, xax=1,/xlo,xs=1,xtit='Frequency (GHz)',color=cgcolor('black'),ticksize_font=1.5
362 365 if DATA NE 0 then begin
363 366 if HARD(0) NE '0' then begin
364 367 oploterror,xg,yg,0.,eg,psym=5
... ... @@ -381,7 +384,7 @@ if keyword_set(data_dir) then data_path = data_dir
381 384 endelse
382 385 endif
383 386 if DATA EQ 1 then begin
384   - OPLOT, LAMBDAf, SP1, lines=1 ; FIR blackbody
  387 + OPLOT, LAMBDAf, SP1, lines=1,color=cgcolor('black') ; FIR blackbody
385 388 sp1_firas = INTERPOL( sp1, lambdaf, wfirasf )
386 389 firas_chi2 = DUSTEM_CHI2( firasf, sp1_firas, 3, err=ufirasf )
387 390 print,'Grey body: Td = ',strtrim(bbpar(0),2),' and beta = ',strtrim(bbpar(1),2)
... ... @@ -390,7 +393,7 @@ if keyword_set(data_dir) then data_path = data_dir
390 393 if keyword_set( CMB ) then begin
391 394 lambda_cmb = 1.d2^( 1.d0 + dindgen(500)*0.01 )
392 395 sp_cmb = 1d3 * lambda_cmb * DUSTEM_BLACKBODY(LAMBDA_CMB, 2.728, unit='w') ; SED in erg/s/cm2/s/sr
393   - OPLOT, LAMBDA_CMB, SP_CMB, lines=0 ; CMB
  396 + OPLOT, LAMBDA_CMB, SP_CMB, color=cgcolor('black'),lines=0 ; CMB
394 397 endif
395 398 if keyword_set( COSMO ) AND DATA NE 0 then begin
396 399 oploterror, x_pep_ame, y_pep_ame, y_pep_ame*0d, e_pep_ame, ps=4
... ... @@ -473,15 +476,15 @@ if keyword_set(data_dir) then data_path = data_dir
473 476 xpr = xpr*[dx,dx^4]
474 477 ; overlay model SED in erg/s/cm2/sr
475 478 for i=0,ntype-1 do begin
476   - oplot, x, sed(*,i), lin=ls(i+1)
477   - oplot,xpr,yps*[1.,1.], lin=ls(i+1)
478   - xyouts,/data,xpr(1)*1.1,yps,gtype(i),chars=1
  479 + oplot, x, sed(*,i), lin=ls(i+1),color=cgcolor('black')
  480 + oplot,xpr,yps*[1.,1.], lin=ls(i+1),color=cgcolor('black')
  481 + xyouts,/data,xpr(1)*1.1,yps,gtype(i),chars=1.3,color=cgcolor('black')
479 482 yps = yps/dy
480 483 endfor
481   - oplot, x, sed(*,ntype),lin=ls(0),col=2
  484 + oplot, x, sed(*,ntype),lin=ls(0),col=cgcolor('blue')
482 485 s1 = STRMID(STRTRIM(ROUND(1e1*g0)/1e1,2), 0, 4)
483 486 s2 = STRMID(STRTRIM(ROUND(1e2*alog10(nh))/1e2,2), 0, 5)
484   - xyouts,/norm,0.6,0.85,'G!d0!n='+s1+' N!dH!n=10!u'+s2+'!ncm!u-2!n'
  487 + xyouts,/norm,0.6,0.85,'G!d0!n='+s1+' N!dH!n=10!u'+s2+'!ncm!u-2!n',color=cgcolor('black'),chars=1.3
485 488 ; overlay model polarized SED in erg/s/cm2/sr
486 489 ;IF (C_POLSED EQ 1) THEN oplot, x, sed_p(*,ntype),lin=ls(0),col=3
487 490 endif
... ... @@ -550,11 +553,11 @@ if keyword_set(data_dir) then data_path = data_dir
550 553 ; then other instruments
551 554 if ctg GT 0 then begin
552 555 for k=0,ctg-1 do begin
553   - oplot, smdat.(itg(k)).x, smdat.(itg(k)).ym(*,ntype), ps=6, syms=1.5, col=3
  556 + oplot, smdat.(itg(k)).x, smdat.(itg(k)).ym(*,ntype), ps=6, syms=1.5, col=cgcolor('red')
554 557 endfor
555 558 endif
556 559 endif
557   -
  560 +
558 561 ;
559 562 ; get extinction
560 563 ;
... ... @@ -667,9 +670,10 @@ if keyword_set(data_dir) then data_path = data_dir
667 670 ip = WHERE( STRPOS(show,'extuv') GE 0, cp )
668 671 IF (SHOW(0) NE '0') AND (CP EQ 1) THEN BEGIN
669 672 if HARD(0) NE '0' then begin
670   - device, file=hard(ihard), /portrait,/color
  673 + device, file="extuv.ps", /portrait,/color
671 674 ihard = ihard + 1
672 675 endif else begin
  676 + cleanplot
673 677 window, wn(1), xs=600,ys=400, tit='DUSTEM_SHOW_FORTRAN: Extinction UV '+strtrim(wn(1),2)
674 678 endelse
675 679 xtit = 'x (!4l!3m!u-1!n)'
... ... @@ -688,11 +692,11 @@ if keyword_set(data_dir) then data_path = data_dir
688 692 yscl = 1d0
689 693 ytit='!4r!3!dext!n (10!u-21!n cm!u2!n per H)'
690 694 endelse
691   - plot, [1d], [1d], lin=ls(0), xr=xr,/xs,xtit=xtit, yr=yr,/ys,ytit=ytit, tit=tit
692   - oplot, x, sext(*,ntype)*yscl, lin=ls(0), col=2
693   - oplot, 1d/x_sm, e_sm*1d21*yscl, ps=4
  695 + cgplot, [1d], [1d], lin=ls(0), xr=xr,/xs,xtit=xtit, yr=yr,/ys,ytit=ytit, tit=tit+" Ext UV",col=cgcolor('black')
  696 + oplot, x, sext(*,ntype)*yscl, lin=ls(0), col=cgcolor('blue')
  697 + oplot, 1d/x_sm, e_sm*1d21*yscl, ps=4,col=cgcolor('red')
694 698 ; oploterror, 1d/x_sm, e_sm*1d21, smdat.i_sm79.err*1d21, psym = 4, /nohat
695   - for i=0,ntype-1 do oplot, x, (sext(*,i))*yscl, lin=ls(i+1)
  699 + for i=0,ntype-1 do oplot, x, (sext(*,i))*yscl, lin=ls(i+1),col=cgcolor('black')
696 700 xyouts,/norm,0.2,0.8,'R!dV!n='+STRMID(STRTRIM(ROUND(1e1*rv)/1e1,2), 0, 4)
697 701 ENDIF
698 702  
... ... @@ -700,9 +704,10 @@ if keyword_set(data_dir) then data_path = data_dir
700 704 ip = WHERE( STRPOS(show,'extir') GE 0, cp )
701 705 IF (SHOW(0) NE '0') AND (CP EQ 1) THEN BEGIN
702 706 if HARD(0) NE '0' then begin
703   - device, file=hard(ihard), /portrait,/color
  707 + device, file="extir.ps", /portrait,/color
704 708 ihard = ihard + 1
705 709 endif else begin
  710 + cleanplot
706 711 window, wn(2), xs=600,ys=400, tit='DUSTEM_SHOW_FORTRAN: Extinction IR '+strtrim(wn(2),2)
707 712 endelse
708 713 xtit = textoidl('\lambda (\mum)')
... ... @@ -710,11 +715,12 @@ if keyword_set(data_dir) then data_path = data_dir
710 715 yr = [1e-6,1.]
711 716 if yscl NE 1d0 THEN ytit='Normalized !4r!3!dext!n' ELSE $
712 717 ytit='!4r!3!dext!n (10!u-21!n cm!u2!n per H)'
713   - plot_oo, [1d], [1d], xr=xr,/xs,xtit=xtit, yr=yr,/ys,ytit=ytit, tit=tit
714   - oplot, wlm, sext(*,ntype)*yscl, lin=ls(0), col=2
715   - oplot, x_sm, e_sm*1d21*yscl, ps=4
  718 +; plot_oo, [1d], [1d], xr=xr,/xs,xtit=xtit, yr=yr,/ys,ytit=ytit, tit=tit
  719 + cgplot, [1d], [1d], xr=xr,/xs,xtit=xtit, yr=yr,/ys,ytit=ytit, tit=tit+" Ext IR",col=cgcolor('black'),/xlo,/ylo
  720 + oplot, wlm, sext(*,ntype)*yscl, lin=ls(0), col=cgcolor('blue')
  721 + oplot, x_sm, e_sm*1d21*yscl, ps=4,col=cgcolor('red')
716 722 ; oploterror, x_sm, e_sm*1d21, smdat.i_sm79.err*1d21, psym = 4, /nohat
717   - for i=0,ntype-1 do oplot, wlm, sext(*,i)*yscl, lin=ls(i+1)
  723 + for i=0,ntype-1 do oplot, wlm, sext(*,i)*yscl, lin=ls(i+1),col=cgcolor('black')
718 724 xyouts,/norm,0.8,0.85,'R!dV!n='+STRMID(STRTRIM(ROUND(1e1*rv)/1e1,2), 0, 4)
719 725 ENDIF
720 726  
... ... @@ -724,18 +730,19 @@ if keyword_set(data_dir) then data_path = data_dir
724 730 ip = WHERE( STRPOS(show,'alb') GE 0, cp )
725 731 IF (SHOW(0) NE '0') AND (CP EQ 1) THEN BEGIN
726 732 if HARD(0) NE '0' then begin
727   - device, file=hard(ihard), /portrait,/color
  733 + device, file="alb.ps", /portrait,/color
728 734 ihard = ihard + 1
729 735 endif else begin
730   - window, wn(3), xs=500,ys=350, tit='DUSTEM_SHOW_FORTRAN: Albedo '+strtrim(wn(3),2)
  736 + cleanplot
  737 + window, wn(3), xs=600,ys=400, tit='DUSTEM_SHOW_FORTRAN: Albedo '+strtrim(wn(3),2)
731 738 endelse
732 739 xtit = 'x (!4l!3m!u-1!n)'
733 740 xr = [0,11]
734 741 ytit='Albedo'
735 742 yr=[ 0, 1]
736   - plot, x, alb(*,ntype), xr=xr,/xs,xtit=xtit, yr=yr,/ys,ytit=ytit, tit=tit
737   - oplot, x, alb(*,ntype), col=2
738   - for i=0,ntype-1 do oplot, x, alb(*,i), lin=ls(i+1)
  743 + cgplot, x, alb(*,ntype), xr=xr,/xs,xtit=xtit, yr=yr,/ys,ytit=ytit, tit=tit+" Albedo",/nodata,col=cgcolor('black')
  744 + oplot, x, alb(*,ntype), col=cgcolor('blue')
  745 + for i=0,ntype-1 do oplot, x, alb(*,i), lin=ls(i+1),col=cgcolor('black')
739 746 ENDIF
740 747  
741 748  
... ... @@ -791,19 +798,21 @@ if keyword_set(data_dir) then data_path = data_dir
791 798 endif
792 799  
793 800 if HARD(0) NE '0' then begin
794   - device, file=hard(nplots-1), /portrait,/color
  801 + device, file="sdist.ps", /portrait,/color
795 802 endif else begin
796   - window, wn(5), xs=500,ys=350, xpos=0,ypos=0,tit='DUSTEM_SHOW_FORTRAN: Size Distribution '+strtrim(wn(5),2)
  803 + cleanplot
  804 + window, wn(5), xs=600,ys=400, xpos=0,ypos=0,tit='DUSTEM_SHOW_FORTRAN: Size Distribution '+strtrim(wn(5),2)
797 805 endelse
798 806 yscl = 1.D29
799 807 xr=[0.1,1e4]
800 808 yr=[-3,3]
801 809 xtit='a (nm)'
802 810 ytit='log( 10!u29!n n!dH!u-1!n a!u4!n dn/da (cm!u3!n/H)) '
803   - tit='DustEM'
804   - plot_oi, ax(0:nsize(0)-1,0)*1.e7, alog10(ava(0:nsize(0)-1,0)*yscl), line=ls(1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit='DUSTEM'
  811 +; tit='DustEM'
  812 +; plot_oi, ax(0:nsize(0)-1,0)*1.e7, alog10(ava(0:nsize(0)-1,0)*yscl), line=ls(1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit='DUSTEM'
  813 + cgplot, ax(0:nsize(0)-1,0)*1.e7, alog10(ava(0:nsize(0)-1,0)*yscl), line=ls(1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit=tit+" Size Dist",col=cgcolor('black'),/xlo
805 814 FOR i=1,ntype-1 DO begin
806   - oplot, ax(0:nsize(i)-1,i)*1.e7, alog10(ava(0:nsize(i)-1,i)*yscl), line=ls(i+1)
  815 + oplot, ax(0:nsize(i)-1,i)*1.e7, alog10(ava(0:nsize(i)-1,i)*yscl), line=ls(i+1),col=cgcolor('black')
807 816 ENDFOR
808 817  
809 818 ENDIF
... ... @@ -814,9 +823,10 @@ if keyword_set(data_dir) then data_path = data_dir
814 823 ip = WHERE( STRPOS(show,'polsed') GE 0, c_polsed )
815 824 IF (SHOW(0) NE '0') AND (C_POLsed EQ 1) THEN BEGIN
816 825 if HARD(0) NE '0' then begin
817   - device, file=hard(ihard), /portrait,/color
  826 + device, file="polsed.ps", /portrait,/color
818 827 ihard = ihard + 1
819 828 endif else begin
  829 + cleanplot
820 830 window, wn(6), xs=600,ys=400, tit='DUSTEM_SHOW_FORTRAN: Polarized SED '+strtrim(wn(6),2)
821 831 endelse
822 832  
... ... @@ -849,15 +859,16 @@ if keyword_set(data_dir) then data_path = data_dir
849 859 xtit = 'Wavelength (!4l!3m)'
850 860 ytit = '!4m!3 P!d!4m!3!n (erg s!u-1!n cm!u-2!n sr!u-1!n)'
851 861 fine = 1
852   - plot_oo, INDGEN(1),/NODAT,xs=9,/ys,XR=xr,YR=yr,XTIT=xtit,YTIT= ytit
853   - axis, /data, xr=1d-5*clight/xr, xax=1,/xlo,xs=1,xtit='Frequency (GHz)'
  862 +; plot_oo, INDGEN(1),/NODAT,xs=9,/ys,XR=xr,YR=yr,XTIT=xtit,YTIT= ytit
  863 + cgplot, INDGEN(1),/NODAT,/xs,/ys,XR=xr,YR=yr,XTIT=xtit,YTIT= ytit,tit=tit+" PolSED",col=cgcolor('black'),/xlo,/ylo
  864 +; axis, /data, xr=1d-5*clight/xr, xax=1,/xlo,xs=1,xtit='Frequency (GHz)',col=cgcolor('black')
854 865 for i=0,ntype-1 do begin
855   - oplot, x, sed_p(*,i), lin=ls(i+1)
856   - oplot,xpr,yps*[1.,1.], lin=ls(i+1)
857   - xyouts,/data,xpr(1)*1.1,yps,gtype(i),chars=1
  866 + oplot, x, sed_p(*,i), lin=ls(i+1),col=cgcolor('black')
  867 + oplot,xpr,yps*[1.,1.], lin=ls(i+1),col=cgcolor('black')
  868 + xyouts,/data,xpr(1)*1.1,yps,gtype(i),chars=1.3,col=cgcolor('black')
858 869 yps = yps/dy
859 870 endfor
860   - oplot,x,sed_p(*,ntype),lin=ls(0),col=3
  871 + oplot,x,sed_p(*,ntype),lin=ls(0),col=cgcolor('blue')
861 872 ENDIF
862 873  
863 874  
... ... @@ -905,9 +916,10 @@ if keyword_set(data_dir) then data_path = data_dir
905 916 smdat.m_ext.sca_p = spsca
906 917  
907 918 if HARD(0) NE '0' then begin
908   - device, file=hard(ihard), /portrait,/color
  919 + device, file="serkowski.ps", /portrait,/color
909 920 ihard = ihard + 1
910 921 endif else begin
  922 + cleanplot
911 923 window, wn(7), xs=600,ys=400, xpos=0,ypos=0,tit='DUSTEM_SHOW_FORTRAN: Serkowski '+strtrim(wn(7),2)
912 924 endelse
913 925  
... ... @@ -918,10 +930,12 @@ if keyword_set(data_dir) then data_path = data_dir
918 930 yr=[0.1,5]
919 931 xtit = textoidl('1/\lambda (\mum^{-1})')
920 932 ytit = textoidl('\sigma_{pol} (10^{-23}cm^2/H)')
921   - tit='DUSTEM'
922   - plot_oo, 1./x, yscl2*(spabs(*,0)+spsca(*,0)), line=ls(1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit=tit
923   - FOR i=1,ntype-1 DO oplot, 1./x, yscl2*(spabs(*,i)+spsca(*,i)), line=ls(i+1)
924   - oplot, 1./x, yscl2*(spabs(*,ntype)+spsca(*,ntype)), col=3
  933 +; tit='DUSTEM: Serkowski'
  934 +; plot_oo, 1./x, yscl2*(spabs(*,0)+spsca(*,0)), line=ls(1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit=tit
  935 + cgplot, 1./x, yscl2*(spabs(*,0)+spsca(*,0)), line=ls(1),xr=xr,/xs $
  936 + ,xtit=xtit,/ys,yr=yr,ytit=ytit,tit=tit+" Serkwoski",col=cgcolor('black'),/xlo,/ylo
  937 + FOR i=1,ntype-1 DO oplot, 1./x, yscl2*(spabs(*,i)+spsca(*,i)), line=ls(i+1),col=cgcolor('black')
  938 + oplot, 1./x, yscl2*(spabs(*,ntype)+spsca(*,ntype)), col=cgcolor('blue')
925 939 ix = WHERE(x_sm LE 7d, csm) ; keep wave below 7 microns
926 940 if csm GT 0 then begin
927 941 xx = 1./x_sm(ix)
... ... @@ -930,7 +944,7 @@ if keyword_set(data_dir) then data_path = data_dir
930 944 ix = SORT(xx)
931 945 xx = xx(ix) & xe = xe(ix)
932 946 fpol = DUSTEM_SERKOWSKI(xx)
933   - oplot,xx,yscl*fpol/5.8e21*3.1,ps=4
  947 + oplot,xx,yscl*fpol/5.8e21*3.1,ps=4,col=cgcolor('red')
934 948 ; oplot,xx,25*fpol*3.1/5.8e21/xe, ps=5
935 949 ENDELSE
936 950  
... ... @@ -942,22 +956,24 @@ if keyword_set(data_dir) then data_path = data_dir
942 956 IF (SHOW(0) NE '0') AND (C_POLSED EQ 1) THEN BEGIN
943 957 ; plot fractional polarisation
944 958 if HARD(0) NE '0' then begin
945   - device, file=hard(ihard), /portrait,/color
  959 + device, file="polfrac.ps", /portrait,/color
946 960 ihard = ihard + 1
947   - endif else begin
  961 + endif else begin
  962 + cleanplot
948 963 window, wn(8), xs=600,ys=400, xpos=0,ypos=0,tit='DUSTEM_SHOW_FORTRAN: P/I '+strtrim(wn(8),2)
949 964 endelse
950 965 xr=[10,1.e4]
951 966 yr=[0.,0.4]
952 967 xtit = textoidl('\lambda (\mum)')
953 968 ytit = textoidl('P/I')
954   - tit=''
955   - plot_oi, x, sed_p(*,0)/sed(*,0), line=ls(1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit=tit
956   - axis, /data, xr=1d-5*clight/xr, xax=1,/xlo,xs=1,xtit='Frequency (GHz)'
  969 +; tit=''
  970 +; plot_oi, x, sed_p(*,0)/sed(*,0), line=ls(1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit=tit
  971 + cgplot, x, sed_p(*,0)/sed(*,0), line=ls(1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit=tit+ "PolFrac",/xlo
  972 +; axis, /data, xr=1d-5*clight/xr, xax=1,/xlo,xs=1,xtit='Frequency (GHz)'
957 973 FOR i=1,ntype-1 DO begin
958   - oplot, x, sed_p(*,i)/sed(*,i), line=ls(i+1)
  974 + oplot, x, sed_p(*,i)/sed(*,i), line=ls(i+1),col=cgcolor('black')
959 975 ENDFOR
960   - oplot, x, sed_p(*,ntype)/sed(*,ntype),col=3
  976 + oplot, x, sed_p(*,ntype)/sed(*,ntype),col=cgcolor('blue')
961 977 ENDIF
962 978  
963 979 ;
... ... @@ -966,9 +982,10 @@ if keyword_set(data_dir) then data_path = data_dir
966 982 IF (SHOW(0) NE '0') AND (C_POLEXT EQ 1) THEN BEGIN
967 983 ; plot fractional polarisation
968 984 if HARD(0) NE '0' then begin
969   - device, file=hard(ihard), /portrait,/color
  985 + device, file="polext.ps", /portrait,/color
970 986 ihard = ihard + 1
971   - endif else begin
  987 + endif else begin
  988 + cleanplot
972 989 window, wn(9), xs=600,ys=400, xpos=0,ypos=0,tit='DUSTEM_SHOW_FORTRAN: p/tau '+strtrim(wn(9),2)
973 990 endelse
974 991 yscl = 1.
... ... @@ -976,12 +993,13 @@ if keyword_set(data_dir) then data_path = data_dir
976 993 yr=[0.,0.5]
977 994 xtit = textoidl('\lambda (\mum)')
978 995 ytit = textoidl('p/\tau')
979   - tit=''
980   - plot_oi, x, yscl*(spabs(*,0)+spsca(*,0))/sext(*,0), line=ls(1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit=tit
  996 +; tit=''
  997 +; plot_oi, x, yscl*(spabs(*,0)+spsca(*,0))/sext(*,0), line=ls(1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit=tit
  998 + cgplot, x, yscl*(spabs(*,0)+spsca(*,0))/sext(*,0), line=ls(1),xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit=tit+" Pol Ext",/xlo
981 999 FOR i=1,ntype-1 DO begin
982   - oplot, x, yscl*(spabs(*,i)+spsca(*,i))/sext(*,i), line=ls(i+1)
  1000 + oplot, x, yscl*(spabs(*,i)+spsca(*,i))/sext(*,i), line=ls(i+1), col=cgcolor('black')
983 1001 ENDFOR
984   - oplot, x, yscl*(spabs(*,ntype)+spsca(*,ntype))/sext(*,ntype), col=3
  1002 + oplot, x, yscl*(spabs(*,ntype)+spsca(*,ntype))/sext(*,ntype), col=cgcolor('blue')
985 1003 ENDIF
986 1004  
987 1005 ;
... ... @@ -1003,9 +1021,10 @@ if keyword_set(data_dir) then data_path = data_dir
1003 1021 readf, uu, anisG0
1004 1022  
1005 1023 if HARD(0) NE '0' then begin
1006   - device, file=hard(ihard), /portrait,/color
  1024 + device, file="align.ps", /portrait,/color
1007 1025 ihard = ihard + 1
1008   - endif else begin
  1026 + endif else begin
  1027 + cleanplot
1009 1028 window, wn(10), xs=600,ys=400, xpos=0,ypos=0,tit='DUSTEM_SHOW_FORTRAN: f_align '+strtrim(wn(10),2)
1010 1029 endelse
1011 1030  
... ... @@ -1027,7 +1046,8 @@ if keyword_set(data_dir) then data_path = data_dir
1027 1046 yr = [0,1.05]
1028 1047 xtit=textoidl('Grain radius (\mum)')
1029 1048 ytit='Alignment efficiency'
1030   - plot_oi,radius,fpol,xr=xr,yr=yr,/ys,xtit=xtit,ytit=ytit
  1049 +; plot_oi,radius,fpol,xr=xr,yr=yr,/ys,xtit=xtit,ytit=ytit
  1050 + cgplot,radius,fpol,xr=xr,yr=yr,/ys,xtit=xtit,ytit=ytit,tit=tit+" Align Eff",/xlo
1031 1051 endif
1032 1052  
1033 1053 ENDIF
... ... @@ -1051,7 +1071,8 @@ if keyword_set(data_dir) then data_path = data_dir
1051 1071 if HARD(0) NE '0' then begin
1052 1072 device, /close
1053 1073 set_plot,'x'
1054   - print,'(W): hardcopies in ',hard
  1074 +; print,'(W): hardcopies in ',hard
  1075 + print,'(W): Generated postscript files in current directory'
1055 1076 !y.thick = 0
1056 1077 !x.thick = 0
1057 1078 !p.thick = 0
... ...