OV-GSO-DC / dustem-wrapper_idl

Commit b5ccb7069e98fea359ca79d7a37e4867cc72cba3

Authored by Jean-Philippe Bernard
1 parent 906f3c27
Exists in master

improved to fit polarization seds

Showing 36 changed files with 1613 additions and 513 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_cc.pro
  Diff comments   View file @b5ccb70
... ... @@ -89,7 +89,10 @@ cc(*)=1.d0
89 89  
90 90 FOR i=0L,Nbands-1 DO BEGIN
91 91 iinstru=where(inst_tags EQ instrus(i),countinstr)
92   - IF countinstr EQ 0 THEN stop
  92 + IF countinstr EQ 0 THEN BEGIN
  93 + message,'Instrument '+instrus(i)+' not found',/continue
  94 + stop
  95 + ENDIF
93 96 st=(*!dustem_filters).(iinstru)
94 97 ii=(where(st.filter_names EQ filter_names(i)))(0)
95 98 spec0=interpol(spec2,wavein2,st.central_wavelengths(ii))
... ... @@ -158,13 +161,18 @@ FOR i=0L,Nbands-1 DO BEGIN
158 161 ; den=integral(*(st.use_wavelengths(ii)),*(st.use_transmissions(ii))/(*(st.use_wavelengths(ii))^2.5),st.use_wmin(ii),st.use_wmax(ii),/double)
159 162 ; cc(i)=num/den/spec0/(st.central_wavelengths(ii))^2.5
160 163  
161   - alpha_new=0.
  164 + alpha_new=0.
162 165 spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
163 166 num=integral(*(st.use_wavelengths(ii)),spec_int*(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii))),st.use_wmin(ii),st.use_wmax(ii),/double)
164 167 den=integral(*(st.use_wavelengths(ii)),*(st.use_transmissions(ii))/(*(st.use_wavelengths(ii))^(alpha_new+2.)),st.use_wmin(ii),st.use_wmax(ii),/double)
165 168 cc(i)=num/den/spec0*(st.central_wavelengths(ii))^(-(alpha_new+1.))
166 169  
167 170 END
  171 + ELSE:BEGIN
  172 + message,strupcase(fluxconvs(i))+' not recognized',/info
  173 + cc[i]=1.
  174 + ;stop
  175 + END
168 176 ENDCASE
169 177 !dustem_previous_cc=ptr_new(cc)
170 178 ENDIF ELSE BEGIN ;When no new cc calculation is required
... ...
src/idl/dustem_compute_polsed.pro
  Diff comments   View file @b5ccb70
... ... @@ -44,11 +44,11 @@ ENDIF
44 44 IF not keyword_set(stp) THEN BEGIN
45 45 ;st=dustem_run(p_dim)
46 46 dustem_create_functions,p_dim/(*(*!dustem_fit).param_init_values),res=res
47   - st=dustem_run(p_dim)
  47 + stp=dustem_run(p_dim)
48 48 ENDIF
49 49 If keyword_set(result) THEN result=st
50   -fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/stp.wav)*1.e20/1.e7
51   -spec = stp.em_tot * fact
  50 +fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/stp.polsed.wav)*1.e20/1.e7
  51 +spec = stp.polsed.em_tot * fact
52 52 ;stop
53 53 if not isarray(spec) THEN stop
54 54 ;stop
... ... @@ -68,7 +68,7 @@ ind_polsed=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_polsed)
68 68 IF count_polsed NE 0 THEN BEGIN
69 69 filter_names=((*!dustem_data.polsed).filt_names)(ind_polsed)
70 70 ;FOR ii=0L,n_elements(filter_names)-1 DO BEGIN
71   - spolsed=dustem_cc(stp.wav,spec,filter_names,cc=cc)
  71 + spolsed=dustem_cc(stp.sed.wav,spec,filter_names,cc=cc)
72 72 ;ENDFOR
73 73 dustem_polsed(ind_polsed)=spolsed
74 74 ENDIF
... ...
src/idl/dustem_compute_sed.pro
  Diff comments   View file @b5ccb70
1   -FUNCTION dustem_compute_sed,p_dim,st,cont=cont,_extra=extra,out_st=out_st,wave0=wave0
  1 +FUNCTION dustem_compute_sed,p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron,_extra=extra,out_st=out_st,wave0=wave0
2 2  
3 3 ;+
4 4 ; NAME:
... ... @@ -47,8 +47,11 @@ ENDIF
47 47 IF not keyword_set(st) THEN BEGIN
48 48 ;st=dustem_run(p_dim)
49 49 IF !dustem_idl_continuum THEN cont=dustem_create_continuum()
50   - dustem_create_functions,p_dim/(*(*!dustem_fit).param_init_values),cont=cont,res=res
  50 + IF !dustem_idl_freefree THEN freefree=dustem_create_freefree()
  51 + IF !dustem_idl_synchrotron THEN synchrotron=dustem_create_synchrotron()
  52 + dustem_create_functions,p_dim/(*(*!dustem_fit).param_init_values),cont=cont,freefree=freefree,synchrotron=synchrotron,res=res
51 53 st=dustem_run(p_dim)
  54 + ;stop
52 55 ENDIF
53 56  
54 57 ;If keyword_set(result) THEN result=st
... ... @@ -58,6 +61,8 @@ if not keyword_set(wave0) then wave0 = 0.
58 61 fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/st.sed.wav)*1.e20/1.e7
59 62 spec = st.sed.em_tot * fact
60 63 IF !dustem_idl_continuum THEN spec=spec+cont
  64 +IF !dustem_idl_freefree THEN spec=spec+freefree
  65 +IF !dustem_idl_synchrotron THEN spec=spec+synchrotron
61 66 if not isarray(spec) THEN stop
62 67 ;COMPUTE THE MODEL SED
63 68 dustem_sed = (*!dustem_data.sed).values * 0.
... ...
src/idl/dustem_create_continuum.pro
  Diff comments   View file @b5ccb70
1   -FUNCTION dustem_create_continuum,key=key,val=val
2   -
3   - temp=1000.
4   - ampl=1.d-2
5   -
6   - IF keyword_set(key) THEN BEGIN
7   - a=where(key EQ 1,count1)
8   - b=where(key EQ 2,count2)
9   - IF count1 NE 0 then temp=(val(a))(0)
10   - IF count2 NE 0 then ampl=(val(b))(0)
11   -; IF count1 NE 0 then temp=1d0*10.^((val(a))(0))
12   -; IF count2 NE 0 then ampl=1d0*10.^((val(b))(0))
13   -; print, "key=",key
14   -; print, "val=",val
15   - ENDIF
  1 +FUNCTION dustem_create_continuum,key=key,val=val,help=help
  2 +
  3 +;+
  4 +; NAME:
  5 +; dustem_create_continuum
  6 +; PURPOSE:
  7 +; Make products for a given flight
  8 +; CATEGORY:
  9 +; DUSTEM Wrapper
  10 +; CALLING SEQUENCE:
  11 +; cont=dustem_create_continuum([,key=][,val=])
  12 +; INPUTS:
  13 +;
  14 +; OPTIONAL INPUT PARAMETERS:
  15 +; key = input parameter number
  16 +; val = input parameter value
  17 +; OUTPUTS:
  18 +; cont = continuum spectrum (on dustem wavelengths)
  19 +; OPTIONAL OUTPUT PARAMETERS:
  20 +; None
  21 +; ACCEPTED KEY-WORDS:
  22 +; help = if set, print this help
  23 +; COMMON BLOCKS:
  24 +; None
  25 +; SIDE EFFECTS:
  26 +; None
  27 +; RESTRICTIONS:
  28 +; None
  29 +; PROCEDURE:
  30 +; This is a dustem pluggin
  31 +; EXAMPLES
  32 +;
  33 +; MODIFICATION HISTORY:
  34 +; Written by JPB
  35 +;-
  36 +
  37 +IF keyword_set(help) THEN BEGIN
  38 + doc_library,'dustem_create_continuum'
  39 + output=0.
  40 + goto,the_end
  41 +ENDIF
  42 +
  43 +;default values of input parameters
  44 +temp=1000.
  45 +ampl=1.d-2
  46 +
  47 +IF keyword_set(key) THEN BEGIN
  48 + a=where(key EQ 1,count1)
  49 + b=where(key EQ 2,count2)
  50 + IF count1 NE 0 then temp=(val(a))(0)
  51 + IF count2 NE 0 then ampl=(val(b))(0)
  52 +ENDIF
16 53  
17   -; stop
18   -; IF getenv('DUSTEM_WHICH') EQ 'DESERT' THEN BEGIN
19   - IF !dustem_which EQ 'DESERT' THEN BEGIN
20   - lambir=((*!dustem_params).gemissiv.lambir)
21   - ENDIF ELSE BEGIN
22   - lambir=((*!dustem_params).lambda.lambda)
23   - ENDELSE
  54 +IF !dustem_which EQ 'DESERT' THEN BEGIN
  55 + lambir=((*!dustem_params).gemissiv.lambir)
  56 +ENDIF ELSE BEGIN
  57 + lambir=((*!dustem_params).lambda.lambda)
  58 +ENDELSE
  59 +
  60 +;=== normalize spectrum to the requested amplitude.
  61 +norm = max(dustem_planck_function(temp,lambir))
  62 +output = ampl*dustem_planck_function(temp,lambir)/norm
24 63  
25   - norm = max(dustem_planck_function(temp,lambir))
26   - output = ampl*dustem_planck_function(temp,lambir)/norm
27   -
28   - return,output
  64 +the_end:
  65 +RETURN,output
29 66  
30 67 END
... ...
src/idl/dustem_create_freefree.pro 0 → 100644
  Diff comments   View file @b5ccb70
... ... @@ -0,0 +1,107 @@
  1 +FUNCTION dustem_create_freefree,key=key,val=val,help=help
  2 +
  3 +;+
  4 +; NAME:
  5 +; dustem_create_freefree
  6 +; PURPOSE:
  7 +; DUSTEM pluggin to compute free-free emission
  8 +; CATEGORY:
  9 +; DUSTEM Wrapper
  10 +; CALLING SEQUENCE:
  11 +; freefree=dustem_create_freefree([,key=][,val=])
  12 +; INPUTS:
  13 +; None
  14 +; OPTIONAL INPUT PARAMETERS:
  15 +; key = input parameter number
  16 +; val = input parameter value
  17 +; OUTPUTS:
  18 +; freefree = free-free spectrum (on dustem wavelengths)
  19 +; OPTIONAL OUTPUT PARAMETERS:
  20 +; None
  21 +; ACCEPTED KEY-WORDS:
  22 +; help = if set, print this help
  23 +; COMMON BLOCKS:
  24 +; None
  25 +; SIDE EFFECTS:
  26 +; None
  27 +; RESTRICTIONS:
  28 +; None
  29 +; PROCEDURE:
  30 +; This is a dustem pluggin
  31 +; EXAMPLES
  32 +;
  33 +; MODIFICATION HISTORY:
  34 +; Written by JPB
  35 +;-
  36 +
  37 +IF keyword_set(help) THEN BEGIN
  38 + doc_library,'dustem_create_freefree'
  39 + output=0.
  40 + goto,the_end
  41 +ENDIF
  42 +
  43 +;default values of input parameters
  44 +Tgas=10000. ;default gas temperature
  45 +Amplitude=1. ;Amplitude
  46 +
  47 +IF keyword_set(key) THEN BEGIN
  48 + a=where(key EQ 1,count1)
  49 + b=where(key EQ 2,count2)
  50 + IF count1 NE 0 then Tgas=(val(a))(0)
  51 + IF count2 NE 0 then Amplitude=(val(b))(0)
  52 +ENDIF
  53 +
  54 +IF !dustem_which EQ 'DESERT' THEN BEGIN
  55 + lambir=((*!dustem_params).gemissiv.lambir)
  56 +ENDIF ELSE BEGIN
  57 + lambir=((*!dustem_params).lambda.lambda)
  58 +ENDELSE
  59 +
  60 +cmic=3.e14
  61 +nu=cmic/lambir ;Hz
  62 +mjy=1 ;output is in MJy/sr
  63 +lambir_ref=10000.
  64 +
  65 +;stop
  66 +;use_method='Deschenes2008'
  67 +;use_method='WallsGabaud1998'
  68 +;use_method='Halpha_scaled'
  69 +use_method='Dickinson2003_norm'
  70 +
  71 +CASE use_method OF
  72 + 'WallsGabaud1998':BEGIN
  73 + em=1. ;This is a stupid value for the Emission measure. Result will be rescaled based on I_halpha_R anyway
  74 + output=intensity_free_free(nu,Tgas,em,I_halpha_R=I_halpha_R,nHe=nHe,ergs=ergs,mjy=mjy)
  75 + fact=Amplitude/I_halpha_R
  76 + ;=== normalize to the requested Halpha value.
  77 + output=output*fact
  78 + END
  79 + 'WallsGabaud1998_scaled':BEGIN
  80 + em=1. ;This is a stupid value for the Emission measure. Result will be rescaled based on I_halpha_R anyway
  81 + output=intensity_free_free(nu,Tgas,em,I_halpha_R=I_halpha_R,nHe=nHe,ergs=ergs,mjy=mjy)
  82 + norm=interpol(output,lambir,lambir_ref)
  83 + output=output/norm*Amplitude
  84 + END
  85 + 'Deschenes2008':BEGIN
  86 + beta_freefree=2.+1./(10.48+1.5*alog(Tgas/8.e3)-alog(nu/1.e9))
  87 + output=nu^(-1.*beta_freefree)
  88 + norm=interpol(output,lambir,lambir_ref)
  89 + output=output/norm*Amplitude
  90 + END
  91 + 'Dickinson2003_norm':BEGIN
  92 + ;stop
  93 + em=1.
  94 + T4=Tgas/1.e4
  95 + nu_ghz=nu/1.e9
  96 + a=0.366*(nu_ghz/1.e9)^(-0.15)*(alog(4.995*1e-2*(nu_ghz)^(-1.))+1.5*alog(Tgas))
  97 + Tb=8.369e3*a*(nu_ghz)^(-2.)*T4^0.667*10^(0.029*T4)*(1.+0.08) ;in mK for 1 Rayleigh
  98 + convert_mk_mjy, lambir, Tb, I_Mjy, /RJ
  99 + norm=interpol(I_Mjy,lambir,lambir_ref)
  100 + output=I_Mjy/norm*Amplitude
  101 + END
  102 +ENDCASE
  103 +
  104 +the_end:
  105 +RETURN,output
  106 +
  107 +END
... ...
src/idl/dustem_create_functions.pro
  Diff comments   View file @b5ccb70
1   -PRO dustem_create_functions,p_min,cont=cont,res=res
  1 +PRO dustem_create_functions,p_min,cont=cont,freefree=freefree,synchrotron=synchrotron,res=res
2 2  
3 3 p_dim = p_min * (*(*!dustem_fit).param_init_values)
4 4  
... ... @@ -18,14 +18,26 @@ FOR i=0L,n_elements(p_min)-1 DO BEGIN
18 18 key_ionfrac=index
19 19 val_ionfrac=value
20 20 ENDIF
21   - IF ftn NE 'dustem_create_continuum' THEN BEGIN
  21 + IF ftn NE 'dustem_create_continuum' and ftn NE 'dustem_create_freefree' and ftn NE 'dustem_create_synchrotron' THEN BEGIN
22 22 str='res='+ftn+'(key=index,val=value)'
23 23 toto=execute(str)
24 24 IF !dustem_verbose NE 0 THEN message,str,/info
25 25 ENDIF ELSE BEGIN
26   - str='cont='+ftn+'(key=index,val=value)'
27   - toto=execute(str)
28   - IF !dustem_verbose NE 0 THEN message,str,/info
  26 + IF ftn EQ 'dustem_create_continuum' THEN BEGIN
  27 + str='cont='+ftn+'(key=index,val=value)'
  28 + toto=execute(str)
  29 + IF !dustem_verbose NE 0 THEN message,str,/info
  30 + ENDIF
  31 + IF ftn EQ 'dustem_create_freefree' THEN BEGIN
  32 + str='freefree='+ftn+'(key=index,val=value)'
  33 + toto=execute(str)
  34 + IF !dustem_verbose NE 0 THEN message,str,/info
  35 + ENDIF
  36 + IF ftn EQ 'dustem_create_synchrotron' THEN BEGIN
  37 + str='synchrotron='+ftn+'(key=index,val=value)'
  38 + toto=execute(str)
  39 + IF !dustem_verbose NE 0 THEN message,str,/info
  40 + ENDIF
29 41 ENDELSE
30 42 f=f+1 & i=i+count-1
31 43 ENDIF
... ...
src/idl/dustem_create_synchrotron.pro 0 → 100644
  Diff comments   View file @b5ccb70
... ... @@ -0,0 +1,74 @@
  1 +FUNCTION dustem_create_synchrotron,key=key,val=val,help=help
  2 +
  3 +;+
  4 +; NAME:
  5 +; dustem_create_synchrotron
  6 +; PURPOSE:
  7 +; DUSTEM pluggin to compute synchrotron emission
  8 +; CATEGORY:
  9 +; DUSTEM Wrapper
  10 +; CALLING SEQUENCE:
  11 +; synch=dustem_create_synchrotron([,key=][,val=])
  12 +; INPUTS:
  13 +; None
  14 +; OPTIONAL INPUT PARAMETERS:
  15 +; key = input parameter number
  16 +; val = input parameter value
  17 +; OUTPUTS:
  18 +; synch = synch spectrum (on dustem wavelengths)
  19 +; OPTIONAL OUTPUT PARAMETERS:
  20 +; None
  21 +; ACCEPTED KEY-WORDS:
  22 +; help = if set, print this help
  23 +; COMMON BLOCKS:
  24 +; None
  25 +; SIDE EFFECTS:
  26 +; None
  27 +; RESTRICTIONS:
  28 +; None
  29 +; PROCEDURE:
  30 +; This is a dustem pluggin
  31 +; EXAMPLES
  32 +;
  33 +; MODIFICATION HISTORY:
  34 +; Written by JPB
  35 +;-
  36 +
  37 +IF keyword_set(help) THEN BEGIN
  38 + doc_library,'dustem_create_synchrotron'
  39 + output=0.
  40 + goto,the_end
  41 +ENDIF
  42 +
  43 +;default values of input parameters
  44 +s=3 ;power spectrum of CR E distribution
  45 +A=1. ;Synchrotron radiation amplitude at 10 mm
  46 +
  47 +IF keyword_set(key) THEN BEGIN
  48 + a=where(key EQ 1,count1)
  49 + b=where(key EQ 2,count2)
  50 + IF count1 NE 0 then s=(val(a))(0)
  51 + IF count2 NE 0 then A=(val(b))(0)
  52 +ENDIF
  53 +
  54 +IF !dustem_which EQ 'DESERT' THEN BEGIN
  55 + lambir=((*!dustem_params).gemissiv.lambir)
  56 +ENDIF ELSE BEGIN
  57 + lambir=((*!dustem_params).lambda.lambda)
  58 +ENDELSE
  59 +
  60 +cmic=3.e14
  61 +nu=cmic/lambir
  62 +lambir_ref=10000.
  63 +
  64 +;see Deschenes et al 2008, eq 6
  65 +synch=nu^(-(s+3.)/2.)
  66 +norm=interpol(synch,lambir,lambir_ref)
  67 +synch=synch/norm*A
  68 +
  69 +output=synch
  70 +
  71 +the_end:
  72 +RETURN,output
  73 +
  74 +END
... ...
src/idl/dustem_fit_modified_bb_readme.pro
  Diff comments   View file @b5ccb70
... ... @@ -116,7 +116,7 @@ Nitermax=100 ;maximum number of iteration. This is the criterium w
116 116 loadct,13
117 117 !y.range=[1e-4,10] ;This is to ajust plot range from outside the routine
118 118 t1=systime(0,/sec)
119   -res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,function_name='dustem_greybody_mpfit')
  119 +res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,function_name='dustem_greybody_mpfit',/xlog,/ylog)
120 120 t2=systime(0,/sec)
121 121  
122 122 stop
... ... @@ -140,7 +140,7 @@ rchi2=(*!dustem_fit).rchi2
140 140 loadct,13
141 141 dustem_sed_plot,(*(*!dustem_fit).current_param_values), $
142 142 ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty, $
143   - res=res,errors=errors,chi2=chi2,rchi2=rchi2,function_name='dustem_greybody_mpfit'
  143 + res=res,errors=errors,chi2=chi2,rchi2=rchi2,function_name='dustem_greybody_mpfit',/xlog,/ylog
144 144  
145 145 ;stop
146 146  
... ... @@ -175,7 +175,7 @@ IF keyword_set(postcript) THEN BEGIN
175 175 device,filename=ps_file,/color
176 176 ENDIF
177 177 dustem_sed_plot,*(*!dustem_fit).current_param_values, $
178   - ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,function_name='dustem_greybody_mpfit'
  178 + ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,function_name='dustem_greybody_mpfit',/xlog,/ylog
179 179 IF keyword_set(postcript) THEN BEGIN
180 180 device,/close
181 181 set_plot,'X'
... ...
src/idl/dustem_fit_sed_polsed_readme.pro 0 → 100644
  Diff comments   View file @b5ccb70
... ... @@ -0,0 +1,215 @@
  1 +PRO dustem_fit_sed_polsed_readme,postcript=postcript,mode=mode,help=help
  2 +
  3 +;This Readme describes how to fit total intensity and polarization SEDs jointly
  4 +;Obviously, the dustem package must have been installed succesfully
  5 +;see DustemWrap documentation for install instructions.
  6 +
  7 +;+
  8 +; NAME:
  9 +; dustem_fit_sed_polsed_readme
  10 +; PURPOSE:
  11 +; This is an example of how to fit total intensity and polarization SEDs jointly
  12 +; It is meant to be an example to follow when writing your own
  13 +; programs using the dustem IDL wrapper.
  14 +; CATEGORY:
  15 +; Dustem
  16 +; CALLING SEQUENCE:
  17 +; dustem_fit_sed_polsed_readme,postcript=postcript,mode=mode,help=help
  18 +; INPUTS:
  19 +; None
  20 +; OPTIONAL INPUT PARAMETERS:
  21 +; None
  22 +; OUTPUTS:
  23 +; None
  24 +; OPTIONAL OUTPUT PARAMETERS:
  25 +; None
  26 +; ACCEPTED KEY-WORDS:
  27 +; postcript = if set plot is done in DUSTEM/Docs/Figures/Last_dustem_fit.ps
  28 +; help = If set print this help
  29 +; COMMON BLOCKS:
  30 +; None
  31 +; SIDE EFFECTS:
  32 +; None
  33 +; RESTRICTIONS:
  34 +; The dustem idl wrapper must be installed
  35 +; PROCEDURE:
  36 +; None
  37 +; EXAMPLES
  38 +; dustem_fit_sed_polsed_readme
  39 +; MODIFICATION HISTORY:
  40 +; Written by J.P. Bernard June 29th 2011
  41 +; see evolution details on the dustem cvs maintained at CESR
  42 +; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
  43 +;-
  44 +
  45 +
  46 +IF keyword_set(help) THEN BEGIN
  47 + doc_library,'dustem_fit_sed_polsed_readme'
  48 + goto,the_end
  49 +ENDIF
  50 +
  51 +;=== initialise dustem
  52 +;dustem_init,mode='COMPIEGNE_ETAL2010',/pol
  53 +use_mode='G17_MODELA'
  54 +dustem_init,mode=use_mode,/pol
  55 +;dustem_init,/pol
  56 +!dustem_verbose=1
  57 +!dustem_show_plot=1
  58 +;!run_pol=1
  59 +help,!run_pol
  60 +
  61 +;stop
  62 +
  63 +;=== make a fake one out for given filters
  64 +;filters=['DIRBE1','DIRBE2','DIRBE3','DIRBE4','IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3']
  65 +filters=['IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3']
  66 +Nfilt=n_elements(filters)
  67 +sed=dustem_initialize_sed(Nfilt)
  68 +sed.filter=filters
  69 +sed.wave=dustem_filter2wav(filters)
  70 +sed.instru=dustem_filter2instru(filters)
  71 +polsed=sed
  72 +;st=dustem_set_data(sed=sed,polsed=polsed)
  73 +st=dustem_set_data(ext=ext,polsed=polsed,polext=polext,sed=sed,polfrac=polfrac,rchi2_weight=rchi2_weight,f_HI=f_HI)
  74 +
  75 +;=== Set which parameters you want to fit
  76 +pd = [ $
  77 + '(*!dustem_params).G0', $ ;G0
  78 + '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
  79 + '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
  80 + '(*!dustem_params).grains(2).mdust_o_mh' $ ;amCBEx
  81 + ]
  82 +
  83 +p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04]
  84 +
  85 +;=== set initial value of parameters you want to fit
  86 +;iv=p_truth ;exact solution
  87 +iv=p_truth+[0.,3.e-5,1.e-5,-2.e-5] ;shifted
  88 +
  89 +ulimed=[0 , 0 ,0 ,0 ]
  90 +llimed=[1 , 1 ,1 ,1 ]
  91 +llims =[0. ,0. ,0. ,0.]
  92 +
  93 +;== SET THE FITTED PARAMETERS (now done in dustem_init_parinfo)
  94 +dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
  95 +
  96 +;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-chrg-spin'
  97 +;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-chrg-spin'
  98 +;(*!dustem_params).grains[4].TYPE_KEYWORDS='pol-plaw-ed'
  99 +
  100 +help,*!dustem_fit
  101 +
  102 +help,!run_pol
  103 +
  104 +;stop
  105 +
  106 +sed.spec=dustem_compute_sed(p_truth,sst,_extra=extra) ;sst is not set on purpose, for dustem_compute_sed to compute the sed using fortran code
  107 +sed.error=(sed.spec*0.1);>0.1
  108 +
  109 +polsed.spec=dustem_compute_polsed(p_truth,sstp,_extra=extra) ;sst is not set on purpose, for dustem_compute_sed to compute the sed using fortran code
  110 +polsed.error=(polsed.spec*0.1)
  111 +;sed.error[*]=100.
  112 +
  113 +;=== SET THE OBSERVATION STRUCTURE
  114 +
  115 +;ind=where(sed.wave GE 5.)
  116 +ind=where(sed.wave GE 0.1)
  117 +
  118 +st=dustem_set_data(polsed=polsed[ind],sed=sed[ind])
  119 +
  120 +;stop
  121 +
  122 +;=== RUN fit
  123 +;tol=1.e-20
  124 +tol=1.e-10
  125 +xtol=1.e-10
  126 +Nitermax=5 ;maximum number of iteration. This is the criterium which will stop the fit procedure
  127 +
  128 +xrange=[1.,2.e4]
  129 +yrange=[1e-4,1.e3]
  130 +
  131 +loadct,13
  132 +;!y.range=[1e-8,100] ;This is to ajust plot range from outside the routine
  133 +t1=systime(0,/sec)
  134 +res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,xrange=xrange,/xstyle,yrange=yrange,/ysty,/ylog,/xlog,xtit='wavelength [mic]',ytit='Brightness []')
  135 +t2=systime(0,/sec)
  136 +
  137 +print,(res-p_truth)/p_truth*100.
  138 +
  139 +help,(*!dustem_fit),/STRUCTURE
  140 +
  141 +;=== SAVE FIT RESULTS
  142 +;file_out=getenv('DUSTEM_RES')+'DUSTEM_fit_example.sav'
  143 +file_out=!dustem_res+'DUSTEM_polsed_fit_example.sav'
  144 +dustem_save_system_variables,file_out
  145 +message,'Saved '+file_out,/continue
  146 +
  147 +stop
  148 +
  149 +;======================================
  150 +;====You can exit IDL here and re-enter
  151 +;======================================
  152 +file='/tmp/ramdisk8MB/DUSTEM_polsed_fit_example.sav'
  153 +dustem_restore_system_variables,file
  154 +
  155 +;=== Plot best fit
  156 +tit='DUSTEM spinning Example'
  157 +ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
  158 +xtit=textoidl('\lambda (\mum)')
  159 +;errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values)
  160 +;errors=(*(*!dustem_fit).current_param_errors)
  161 +;chi2=(*!dustem_fit).chi2
  162 +;rchi2=(*!dustem_fit).rchi2
  163 +
  164 +;loadct,13
  165 +;window,win & win=win+1
  166 +;dustem_sed_plot,(*(*!dustem_fit).current_param_values), $
  167 +; ytit=ytit,xtit=xtit,title=tit,yr=yrange,xr=xrange,/ysty,/xsty, $
  168 +; res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog
  169 +win=1
  170 +window,win & win=win+1
  171 +xrange=[1.,2.e4]
  172 +yrange=[1e-4,1.e3]
  173 +
  174 +;=== RESTORE FIT RESULTS
  175 +;=== initialise dustem
  176 +;use_mode='G17_MODELA'
  177 +;dustem_init,mode=use_mode,/pol ;CAUTION: MUST use the same mode as earlier !!
  178 +
  179 +;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-chrg-spin'
  180 +;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-chrg-spin'
  181 +
  182 +;=== recover best fit values
  183 +res=*(*!dustem_fit).current_param_values
  184 +chi2=(*!dustem_fit).chi2
  185 +rchi2=(*!dustem_fit).rchi2
  186 +;errors=*(*!dustem_fit).current_param_errors
  187 +errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values)
  188 +
  189 +;=== Plot best fit
  190 +tit='DUSTEM polarization dust Example (Saved)'
  191 +ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
  192 +xtit=textoidl('\lambda (\mum)')
  193 +
  194 +loadct,13
  195 +IF keyword_set(postcript) THEN BEGIN
  196 + set_plot,'PS'
  197 + ps_file=!dustem_wrap_soft_dir+'/Docs/Figures/'+'Last_dustem_spinning_fit.ps'
  198 + device,filename=ps_file,/color
  199 +ENDIF
  200 +dustem_sed_plot,*(*!dustem_fit).current_param_values,ytit=ytit,xtit=xtit,title=tit,yr=yrange,xr=xrange,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog,/pol
  201 +
  202 +
  203 +;dustem_sed_plot,*(*!dustem_fit).current_param_values, $
  204 +; ytit=ytit,xtit=xtit,title=tit,yr=yrange,xr=xrange,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog
  205 +IF keyword_set(postcript) THEN BEGIN
  206 + device,/close
  207 + set_plot,'X'
  208 + message,'wrote '+ps_file,/info
  209 +ENDIF
  210 +
  211 +print,'dustem_mpfit_sed executed in ',t2-t1,' sec'
  212 +
  213 +the_end:
  214 +
  215 +END
... ...
src/idl/dustem_fit_sed_readme.pro
  Diff comments   View file @b5ccb70
... ... @@ -237,7 +237,7 @@ Nitermax=2 ;maximum number of iteration. This is the criterium whi
237 237 loadct,13
238 238 !y.range=[1e-4,10] ;This is to ajust plot range from outside the routine
239 239 t1=systime(0,/sec)
240   -res=dustem_mpfit_data(tol=tol,Nitermax=Nitermax,gtol=gtol)
  240 +res=dustem_mpfit_data(tol=tol,Nitermax=Nitermax,gtol=gtol,/xlog,/ylog)
241 241 ;res=dustem_mpfit_sed(tol=tol,Nitermax=Nitermax) ;,func_name=func_name,cf_min=cf_min)
242 242 t2=systime(0,/sec)
243 243  
... ... @@ -274,7 +274,7 @@ rchi2=(*!dustem_fit).rchi2
274 274 loadct,13
275 275 dustem_sed_plot,(*(*!dustem_fit).current_param_values), $
276 276 ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty, $
277   - res=res,errors=errors,chi2=chi2,rchi2=rchi2
  277 + res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog
278 278  
279 279 ;stop
280 280  
... ... @@ -319,7 +319,7 @@ IF keyword_set(postcript) THEN BEGIN
319 319 device,filename=ps_file,/color
320 320 ENDIF
321 321 ;dustem_sed_plot,(*!dustem_fit_params),ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2
322   -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
  322 +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,/xlog,/ylog
323 323 IF keyword_set(postcript) THEN BEGIN
324 324 device,/close
325 325 set_plot,'X'
... ...
src/idl/dustem_fit_sed_spinning_readme.pro
  Diff comments   View file @b5ccb70
1 1 PRO dustem_fit_sed_spinning_readme,postcript=postcript,mode=mode,help=help
2 2  
  3 +;This is the version with spinning dust activated. Not very useful for the proposed SED
  4 +
3 5 ;This Readme describes how to fit SEDs with dustem with spinning dust turned on.
4 6 ;It runs on the SED stored into the file sample_SED.xcat
5 7 ;in this directory. Remember that the goal here is not necessarily to
... ... @@ -213,7 +215,7 @@ Nitermax=2 ;maximum number of iteration. This is the criterium whi
213 215 loadct,13
214 216 !y.range=[1e-4,10] ;This is to ajust plot range from outside the routine
215 217 t1=systime(0,/sec)
216   -res=dustem_mpfit_data(tol=tol,Nitermax=Nitermax,gtol=gtol)
  218 +res=dustem_mpfit_data(tol=tol,Nitermax=Nitermax,gtol=gtol,/xlog,/ylog,xrange=[1,20000])
217 219 t2=systime(0,/sec)
218 220  
219 221  
... ... @@ -226,7 +228,7 @@ dustem_save_sed_fit,file_out
226 228  
227 229 ;=== Plot best fit
228 230 yr=[1e-4,10]
229   -xr=[1,2000]
  231 +xr=[1,20000]
230 232 tit='DUSTEM Example'
231 233 ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
232 234 xtit=textoidl('\lambda (\mum)')
... ... @@ -237,7 +239,7 @@ rchi2=(*!dustem_fit).rchi2
237 239 loadct,13
238 240 dustem_sed_plot,(*(*!dustem_fit).current_param_values), $
239 241 ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty, $
240   - res=res,errors=errors,chi2=chi2,rchi2=rchi2
  242 + res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog
241 243  
242 244 ;stop
243 245  
... ... @@ -269,7 +271,7 @@ errors=*(*!dustem_fit).current_param_errors
269 271  
270 272 ;=== Plot best fit
271 273 yr=[1e-4,10]
272   -xr=[1,2000]
  274 +xr=[1,20000]
273 275 tit='DUSTEM Example (Saved)'
274 276 ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
275 277 xtit=textoidl('\lambda (\mum)')
... ... @@ -280,7 +282,7 @@ IF keyword_set(postcript) THEN BEGIN
280 282 ps_file=!dustem_wrap_soft_dir+'/Docs/Figures/'+'Last_dustem_fit.ps'
281 283 device,filename=ps_file,/color
282 284 ENDIF
283   -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
  285 +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,/xlog,/ylog
284 286 IF keyword_set(postcript) THEN BEGIN
285 287 device,/close
286 288 set_plot,'X'
... ...
src/idl/dustem_fit_spin_test_readme.pro
  Diff comments   View file @b5ccb70
... ... @@ -20,7 +20,7 @@ PRO dustem_fit_spin_test_readme,postcript=postcript,mode=mode,help=help
20 20 ; CATEGORY:
21 21 ; Dustem
22 22 ; CALLING SEQUENCE:
23   -; dustem_fit_modified_bb_readme,postcript=postcript,help=help
  23 +; dustem_fit_spin_test_readme,postcript=postcript,help=help
24 24 ; INPUTS:
25 25 ; None
26 26 ; OPTIONAL INPUT PARAMETERS:
... ... @@ -41,7 +41,7 @@ PRO dustem_fit_spin_test_readme,postcript=postcript,mode=mode,help=help
41 41 ; PROCEDURE:
42 42 ; None
43 43 ; EXAMPLES
44   -; dustem_fit_modified_bb_readme
  44 +; dustem_fit_spin_test_readme
45 45 ; MODIFICATION HISTORY:
46 46 ; Written by J.P. Bernard June 9th 2011
47 47 ; see evolution details on the dustem cvs maintained at CESR
... ... @@ -60,7 +60,7 @@ dustem_init,mode='COMPIEGNE_ETAL2010'
60 60 !dustem_show_plot=1
61 61  
62 62 ;=== make a fake one out for given filters
63   -filters=['IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3']
  63 +filters=['DIRBE1','DIRBE2','DIRBE3','DIRBE4','IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3']
64 64 Nfilt=n_elements(filters)
65 65 sed=dustem_initialize_sed(Nfilt)
66 66 sed.filter=filters
... ... @@ -75,6 +75,8 @@ st=dustem_set_data(sed=sed)
75 75 ;amCBEx 25 plaw-ed 1.4500E-03 1.8100E+00 4.0000E-07 2.0000E-04 -2.8000E+00 1.5000E-05 1.5000E-05 2.0000E+00
76 76 ;aSilx 25 plaw-ed 6.7000E-03 3.0000E+00 4.0000E-07 2.0000E-04 -3.4000E+00 2.0000E-05 2.0000E-05 2.0000E+00
77 77  
  78 +goto,freefree_tag
  79 +
78 80 ;=== Set which parameters you want to fit
79 81 pd = [ $
80 82 ; '(*!dustem_params).G0', $ ;G0
... ... @@ -86,23 +88,89 @@ pd = [ $
86 88 '(*!dustem_params).gas.G0', $ ;Gas G0
87 89 '(*!dustem_params).gas.Tgas', $ ;Gas temperature
88 90 '(*!dustem_params).gas.nH', $ ;Gas temperature
  91 +; 'dustem_create_freefree_1', $ ;Gas temperature for free-free emission
  92 +; 'dustem_create_freefree_2', $ ;H-alpha for free-free emission (R)
89 93 'dustem_create_continuum_2'] ;Intensity of NIR continuum
90 94  
91   -p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,100.,1000.,1000.,0.001]
  95 +;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,100.,1000.,1000.,0.001]
  96 +;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,1.e4,0.01,0.1]
  97 +p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.1]
92 98  
93 99 ;=== set initial value of parameters you want to fit
94   -iv=p_truth ;exact solution
95   -iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,90.,500.,600.,0.001] ;shifted
  100 +;iv=p_truth ;exact solution
  101 +;iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,1.e4,0.01,0.001] ;shifted
  102 +iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.001] ;shifted
96 103  
97 104 ulimed=[0 , 0 ,0,0,0,0,0,0,0]
98 105 llimed=[1 , 1 ,1,1,1,1,1,1,1]
99 106 llims =[0. ,0. ,0.,0.,0.,0.,0.,0.,0.]
  107 +goto,out_of_there
  108 +
  109 +freefree_tag:
  110 +;=== Set which parameters you want to fit
  111 +pd = [ $
  112 +; '(*!dustem_params).G0', $ ;G0
  113 + '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
  114 + '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
  115 + '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx
  116 + '(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx
  117 + '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSilx
  118 + '(*!dustem_params).gas.G0', $ ;Gas G0
  119 + '(*!dustem_params).gas.Tgas', $ ;Gas temperature
  120 + '(*!dustem_params).gas.nH', $ ;Gas temperature
  121 + 'dustem_create_freefree_1', $ ;Gas temperature for free-free emission
  122 + 'dustem_create_freefree_2', $ ;H-alpha for free-free emission (R)
  123 + 'dustem_create_continuum_2'] ;Intensity of NIR continuum
  124 +
  125 +;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,100.,1000.,1000.,0.001]
  126 +p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,1.e4,1.e-11,0.1]
  127 +;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.1]
  128 +
  129 +;=== set initial value of parameters you want to fit
  130 +;iv=p_truth ;exact solution
  131 +iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,1.e4,1.e-11,0.001] ;shifted
  132 +;iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.001] ;shifted
  133 +
  134 +ulimed=[0 , 0 ,0,0,0,0,0,0,0,0,0]
  135 +llimed=[1 , 1 ,1,1,1,1,1,1,1,1,1]
  136 +llims =[0. ,0. ,0.,0.,0.,0.,0.,0.,0.,0,0]
  137 +
  138 +synchrotron_tag:
  139 +;=== Set which parameters you want to fit
  140 +pd = [ $
  141 +; '(*!dustem_params).G0', $ ;G0
  142 + '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
  143 + '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
  144 + '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx
  145 + '(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx
  146 + '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSilx
  147 + '(*!dustem_params).gas.G0', $ ;Gas G0
  148 + '(*!dustem_params).gas.Tgas', $ ;Gas temperature
  149 + '(*!dustem_params).gas.nH', $ ;Gas temperature
  150 + 'dustem_create_synchrotron_1', $ ;CR spectral index for synchrotron emission
  151 + 'dustem_create_synchrotron_2', $ ;Amplitude for synchrotron emission (R)
  152 + 'dustem_create_continuum_2'] ;Intensity of NIR continuum
  153 +
  154 +;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,100.,1000.,1000.,0.001]
  155 +p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,3.,0.01,0.1]
  156 +;p_truth=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.1]
  157 +
  158 +;=== set initial value of parameters you want to fit
  159 +;iv=p_truth ;exact solution
  160 +iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,3.,0.01,0.1] ;shifted
  161 +;iv=[7.8000E-04,7.8000E-04,1.6500E-04,1.4500E-03,6.7000E-03,1.e3,1.e5,100.,0.001] ;shifted
  162 +
  163 +ulimed=[0 , 0 ,0,0,0,0,0,0,0,0,0]
  164 +llimed=[1 , 1 ,1,1,1,1,1,1,1,1,1]
  165 +llims =[0. ,0. ,0.,0.,0.,0.,0.,0.,0.,0,0]
  166 +
  167 +out_of_there:
100 168  
101 169 ;== SET THE FITTED PARAMETERS (now done in dustem_init_parinfo)
102 170 dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
103 171  
104   -(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-spin'
105   -(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-spin'
  172 +(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-chrg-spin'
  173 +(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-chrg-spin'
106 174  
107 175 help,*!dustem_fit
108 176  
... ... @@ -114,7 +182,8 @@ sed.error=(sed.spec*0.1);>0.1
114 182  
115 183 ;=== SET THE OBSERVATION STRUCTURE
116 184  
117   -ind=where(sed.wave GE 5.)
  185 +;ind=where(sed.wave GE 5.)
  186 +ind=where(sed.wave GE 0.1)
118 187 st=dustem_set_data(sed=sed[ind])
119 188  
120 189  
... ... @@ -124,25 +193,26 @@ tol=1.e-10
124 193 xtol=1.e-10
125 194 Nitermax=5 ;maximum number of iteration. This is the criterium which will stop the fit procedure
126 195  
  196 +xrange=[1.,2.e4]
  197 +yrange=[1e-4,1.e3]
  198 +
127 199 loadct,13
128 200 ;!y.range=[1e-8,100] ;This is to ajust plot range from outside the routine
129 201 t1=systime(0,/sec)
130   -res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,xrange=[1.,2.e4],/xstyle,yrange=[1e-4,1.e3],/ysty,/ylog,/xlog,xtit='wavelength [mic]',ytit='Brightness []')
  202 +res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,xrange=xrange,/xstyle,yrange=yrange,/ysty,/ylog,/xlog,xtit='wavelength [mic]',ytit='Brightness []')
131 203 t2=systime(0,/sec)
132 204  
133 205 print,(res-p_truth)/p_truth*100.
134 206  
135   -stop
  207 +;stop
136 208  
137 209 ;=== SAVE FIT RESULTS
138 210 ;file_out=getenv('DUSTEM_RES')+'DUSTEM_fit_example.sav'
139   -file_out=!dustem_res+'DUSTEM_bb_fit_example.sav'
  211 +file_out=!dustem_res+'DUSTEM_spinning_fit_example.sav'
140 212 dustem_save_sed_fit,file_out
141 213  
142 214 ;=== Plot best fit
143   -yr=[1e-3,100]
144   -xr=[50,4000]
145   -tit='DUSTEM Modified BB Example'
  215 +tit='DUSTEM spinning Example'
146 216 ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
147 217 xtit=textoidl('\lambda (\mum)')
148 218 ;errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values)
... ... @@ -152,8 +222,8 @@ rchi2=(*!dustem_fit).rchi2
152 222  
153 223 loadct,13
154 224 dustem_sed_plot,(*(*!dustem_fit).current_param_values), $
155   - ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty, $
156   - res=res,errors=errors,chi2=chi2,rchi2=rchi2,function_name='dustem_greybody_mpfit'
  225 + ytit=ytit,xtit=xtit,title=tit,yr=yrange,xr=xrange,/ysty,/xsty, $
  226 + res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog
157 227  
158 228 ;stop
159 229  
... ... @@ -161,13 +231,18 @@ dustem_sed_plot,(*(*!dustem_fit).current_param_values), $
161 231 ;====You can exit IDL here and re-enter
162 232 ;======================================
163 233 window,1
  234 +xrange=[1.,2.e4]
  235 +yrange=[1e-4,1.e3]
164 236  
165 237 ;=== RESTORE FIT RESULTS
166 238 ;=== initialise dustem
167   -dustem_init
168   -file=!dustem_res+'DUSTEM_bb_fit_example.sav'
  239 +dustem_init,mode=use_mode ;CAUTION: MUST use the same mode as earlier !!
  240 +file=!dustem_res+'DUSTEM_spinning_fit_example.sav'
169 241 dustem_restore_sed_fit,file
170 242  
  243 +(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-chrg-spin'
  244 +(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-chrg-spin'
  245 +
171 246 ;=== recover best fit values
172 247 res=*(*!dustem_fit).current_param_values
173 248 chi2=(*!dustem_fit).chi2
... ... @@ -175,20 +250,19 @@ rchi2=(*!dustem_fit).rchi2
175 250 errors=*(*!dustem_fit).current_param_errors
176 251  
177 252 ;=== Plot best fit
178   -yr=[1e-3,100]
179   -xr=[50,4000]
180   -tit='DUSTEM Modified BB Example (Saved)'
  253 +tit='DUSTEM spinning dust Example (Saved)'
181 254 ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
182 255 xtit=textoidl('\lambda (\mum)')
183 256  
184 257 loadct,13
185 258 IF keyword_set(postcript) THEN BEGIN
186 259 set_plot,'PS'
187   - ps_file=!dustem_wrap_soft_dir+'/Docs/Figures/'+'Last_dustem_bb_fit.ps'
  260 + ps_file=!dustem_wrap_soft_dir+'/Docs/Figures/'+'Last_dustem_spinning_fit.ps'
188 261 device,filename=ps_file,/color
189 262 ENDIF
190   -dustem_sed_plot,*(*!dustem_fit).current_param_values, $
191   - ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,function_name='dustem_greybody_mpfit'
  263 +dustem_sed_plot,*(*!dustem_fit).current_param_values,ytit=ytit,xtit=xtit,title=tit,yr=yrange,xr=xrange,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog
  264 +;dustem_sed_plot,*(*!dustem_fit).current_param_values, $
  265 +; ytit=ytit,xtit=xtit,title=tit,yr=yrange,xr=xrange,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog
192 266 IF keyword_set(postcript) THEN BEGIN
193 267 device,/close
194 268 set_plot,'X'
... ...
src/idl/dustem_init.pro
  Diff comments   View file @b5ccb70
... ... @@ -52,6 +52,9 @@ IF keyword_set(help) THEN BEGIN
52 52 goto,the_end
53 53 ENDIF
54 54  
  55 +defsysv,'!run_pol',0 ;0=no polar
  56 +IF keyword_set(pol) THEN !run_pol=1
  57 +
55 58 ;=== color correction system variables
56 59 defsysv,'!dustem_previous_cc',ptr_new() ;store previous color correction values
57 60 defsysv,'!dustem_do_cc',1 ;indicates to dustem if color corrections are needed. Will be changed when fitting
... ... @@ -135,6 +138,8 @@ defsysv, '!dustem_params', ptr_new() ;Contains the values of all Desert Model pa
135 138 ;ENDIF
136 139  
137 140 defsysv, '!dustem_idl_continuum', 0.
  141 +defsysv, '!dustem_idl_freefree', 0.
  142 +defsysv, '!dustem_idl_synchrotron', 0.
138 143  
139 144 defsysv,'!dustem_filters',ptr_new() ;filter information (filled by dustem_filter_init.pro)
140 145 ;; IDL> help,*!dustem_filters,/str
... ... @@ -180,12 +185,26 @@ IF keyword_set(mode) THEN BEGIN
180 185 END
181 186 'DL07':BEGIN
182 187 (*!dustem_inputs).grain='GRAIN_DL07.DAT'
183   - END
184   - 'AJ13':BEGIN
  188 + END
  189 + 'AJ13':BEGIN
185 190 (*!dustem_inputs).grain='GRAIN_J13.DAT'
186   - END
187   -
188   -
  191 + END
  192 + 'G17_MODELA':BEGIN
  193 + (*!dustem_inputs).grain='GRAIN_G17_MODELA.DAT'
  194 + (*!dustem_inputs).align='ALIGN_G17_MODELA.DAT'
  195 + END
  196 + 'G17_MODELB':BEGIN
  197 + (*!dustem_inputs).grain='GRAIN_G17_MODELB.DAT'
  198 + (*!dustem_inputs).align='ALIGN_G17_MODELB.DAT'
  199 + END
  200 + 'G17_MODELC':BEGIN
  201 + (*!dustem_inputs).grain='GRAIN_G17_MODELC.DAT'
  202 + (*!dustem_inputs).align='ALIGN_G17_MODELC.DAT'
  203 + END
  204 + 'G17_MODELD':BEGIN
  205 + (*!dustem_inputs).grain='GRAIN_G17_MODELD.DAT'
  206 + (*!dustem_inputs).align='ALIGN_G17_MODELD.DAT'
  207 + END
189 208 ELSE:BEGIN
190 209 message,'mode '+mode+' unknown'
191 210 END
... ...
src/idl/dustem_init_parinfo.pro
  Diff comments   View file @b5ccb70
... ... @@ -3,7 +3,9 @@ PRO dustem_init_parinfo,pd,iv, $
3 3 up_limited=up_limited,lo_limited=lo_limited,up_limits=up_limits,lo_limits=lo_limits, $
4 4 relstep=relstep,step=step
5 5  
6   -(*!dustem_fit).param_descs=ptr_new(pd) & (*!dustem_fit).param_init_values=ptr_new(iv)
  6 +(*!dustem_fit).param_descs=ptr_new(pd)
  7 +(*!dustem_fit).param_init_values=ptr_new(iv)
  8 +;stop
7 9 dustem_set_func_ind,pd,iv
8 10  
9 11  
... ...
src/idl/dustem_mpfit_data.pro
  Diff comments   View file @b5ccb70
... ... @@ -61,7 +61,7 @@ IF keyword_set(Nitermax) THEN use_Nitermax=Nitermax
61 61 ind_data_tag=dustem_data_tag(count=count_data_tag)
62 62  
63 63 fields=['wav','values','sigma']
64   -for j=0,n_elements(fields)-1 do begin
  64 +FOR j=0,n_elements(fields)-1 DO BEGIN
65 65  
66 66 instruc = fields(j)+' = ([0' ; we add 0 first to allow for ',' after it
67 67  
... ... @@ -96,7 +96,7 @@ for j=0,n_elements(fields)-1 do begin
96 96 ;print,'-> INSTRUC: ',instruc
97 97 toto=execute(instruc) ; Variables wav, values & sigma are created
98 98  
99   -endfor
  99 +ENDFOR
100 100  
101 101 ;stop
102 102 p_min = dustem_mpfitfun(func_name,wav,values,sigma, $
... ...
src/idl/dustem_mpfit_run.pro
  Diff comments   View file @b5ccb70
... ... @@ -52,13 +52,15 @@ p_dim=p_min*(*(*!dustem_fit).param_init_values)
52 52 ;Run dustem with as an interface to mpfitfun
53 53  
54 54 IF !dustem_idl_continuum THEN cont=dustem_create_continuum()
  55 +IF !dustem_idl_freefree THEN freefree=dustem_create_freefree()
  56 +IF !dustem_idl_synchrotron THEN synchrotron=dustem_create_synchrotron()
55 57  
56 58 ;CALL THE DUSTEM_CREATE FUNCTIONS
57 59 ;Actually, res is never used
58 60 ;THIS line below was added by NF, but removed by JPB on Dec 5 2009
59 61 ;dustem_set_params,p_dim ;NF on October 2009 (params have to be updated before calling some functions !)
60 62 ;dustem_create_functions,p_dim/(*!IDO),cont=cont,res=res
61   -dustem_create_functions,p_dim/(*(*!dustem_fit).param_init_values),cont=cont,res=res
  63 +dustem_create_functions,p_dim/(*(*!dustem_fit).param_init_values),cont=cont,freefree=freefree,synchrotron=synchrotron,res=res
62 64  
63 65 ;if ISRF has been modified and IONFRAC has not been modified, then call dustem_create_ionfrac
64 66 ;NF on October 2009 : handling of DUSTEM_CREATE_IONFRAC has been changed
... ... @@ -98,7 +100,7 @@ ind_data_tag=dustem_data_tag(count=count_data_tag)
98 100  
99 101 ;if exist('./noplot') then !dustem_show_plot = 0 else if !dustem_show_plot eq 0 then !dustem_show_plot = 1
100 102  
101   -for i_tag=0,count_data_tag-1 do begin
  103 +FOR i_tag=0,count_data_tag-1 DO BEGIN
102 104  
103 105 ;stop
104 106  
... ... @@ -106,247 +108,245 @@ for i_tag=0,count_data_tag-1 do begin
106 108 ;print,tagnames(ind_data_tag(i_tag))
107 109 CASE tagnames(ind_data_tag(i_tag)) OF
108 110  
109   - 'SED' : BEGIN
110   -
111   - dustem_sed = dustem_compute_sed(p_dim,st,cont=cont)
112   - chi2_sed = total(((*!dustem_data.sed).values-dustem_sed)^2/(*!dustem_data.sed).sigma^2)
113   - n_sed = n_elements((*!dustem_data.sed).values)
114   - rchi2_sed = chi2_sed / n_sed
115   - wrchi2_sed = rchi2_sed * !fit_rchi2_weight.sed
116   - message,"chi2 SED = "+strtrim(chi2_sed,2)+" rchi2 SED = "+strtrim(rchi2_sed,2),/continue;," weighted rchi2 SED=",wrchi2_sed
117   - ;print,'wl ',(*!dustem_data.sed).wav
118   - ;print,"Delta SED per wl: ",(*!dustem_data.sed).values-dustem_sed
119   - ;print,"error per wl: ",(*!dustem_data.sed).sigma
120   - ;print,"chi2 per wl: ",((*!dustem_data.sed).values-dustem_sed)^2/(*!dustem_data.sed).sigma^2
121   - alpha = total(((*!dustem_data.sed).values*dustem_sed)/(*!dustem_data.sed).sigma^2) $
122   - / total(((*!dustem_data.sed).values)^2/(*!dustem_data.sed).sigma^2)
123   - beta = total(((*!dustem_data.sed).values*dustem_sed)/(*!dustem_data.sed).sigma^2) $
124   - / total(dustem_sed^2/(*!dustem_data.sed).sigma^2)
125   -
126   - print,"Best SED chi2 obtained if f_HI = ",alpha, " or if we multiply model by ",beta
127   -
128   - ;dustem_plot_nuinu_em,st,dustem_sed,cont,/print_radiance,p_dim=p_dim
129   -
130   - ; Plot if needed
131   - IF !dustem_show_plot NE 0 THEN BEGIN
132   - window,0
133   - dustem_plot_fit_sed,st,dustem_sed,cont,_extra=extra,res=p_min*(*(*!dustem_fit).param_init_values),chi2=(*!dustem_fit).chi2,rchi2=(*!dustem_fit).rchi2
134   - ;ytit='Total intensity'
135   - ;xtit=textoidl('\lambda (\mum)')
136   - ;win+=1
137   - ;xr=[10,5e3]
138   - ;yr=[5e-33,5.e-28]
139   - ;nodiff = 0
140   - ;dustem_plot_nuinu_em_vg,st,dustem_sed,cont,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1,nodiff=nodiff
141   - ;dustem_plot_nuinu_em_vg,st,dustem_sed,cont,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2,nodiff=nodiff
142   -
143   - ENDIF
144   -
145   - dustem_out = [ dustem_out, dustem_sed ]
146   -
147   - END
148   -
149   - 'EXT' : BEGIN
150   -
151   - dustem_ext = interpol(st.ext.ext_tot,st.ext.wav,(*!dustem_data.ext).wav)
152   - chi2_ext = total(((*!dustem_data.ext).values-dustem_ext)^2/(*!dustem_data.ext).sigma^2)
153   - n_ext = n_elements((*!dustem_data.ext).values)
154   - rchi2_ext = chi2_ext / n_ext
155   - wrchi2_ext = rchi2_ext * !fit_rchi2_weight.ext
156   - print,"chi2 EXT = ",chi2_ext," rchi2 EXT = ",rchi2_ext;," weighted rchi2 EXT=",wrchi2_ext
157   - ;print,"chi2 per wl: ",((*!dustem_data.ext).values-dustem_ext)^2/(*!dustem_data.ext).sigma^2
158   -
159   - ; Plot if needed
160   - IF !dustem_show_plot NE 0 THEN BEGIN
161   -
162   - win+=1
163   - ;dustem_plot_ext,st,yr=[1.e-6,1.e0],/ysty,xr=[0.5,40],/xsty,win=win
164   - xr=[0.3,40]
165   - yr=[1.e-26,1.e-21]
166   - dustem_plot_ext,st,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1
167   - dustem_plot_ext,st,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2
168   -
169   - win+=1
170   - ;dustem_plot_ext,st,/UV,yr=[0,2.5],/ysty,xr=[0,10],/xsty,win=win
171   - xr=[0.2,10]
172   - yr=[0,3e-21]
173   - dustem_plot_ext,st,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1
174   - dustem_plot_ext,st,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2
175   -
176   - ENDIF
  111 + 'SED' : BEGIN
  112 +
  113 + dustem_sed = dustem_compute_sed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
  114 + chi2_sed = total(((*!dustem_data.sed).values-dustem_sed)^2/(*!dustem_data.sed).sigma^2)
  115 + n_sed = n_elements((*!dustem_data.sed).values)
  116 + rchi2_sed = chi2_sed / n_sed
  117 + wrchi2_sed = rchi2_sed * !fit_rchi2_weight.sed
  118 + message,"chi2 SED = "+strtrim(chi2_sed,2)+" rchi2 SED = "+strtrim(rchi2_sed,2),/continue;," weighted rchi2 SED=",wrchi2_sed
  119 + ;print,'wl ',(*!dustem_data.sed).wav
  120 + ;print,"Delta SED per wl: ",(*!dustem_data.sed).values-dustem_sed
  121 + ;print,"error per wl: ",(*!dustem_data.sed).sigma
  122 + ;print,"chi2 per wl: ",((*!dustem_data.sed).values-dustem_sed)^2/(*!dustem_data.sed).sigma^2
  123 + alpha = total(((*!dustem_data.sed).values*dustem_sed)/(*!dustem_data.sed).sigma^2) $
  124 + / total(((*!dustem_data.sed).values)^2/(*!dustem_data.sed).sigma^2)
  125 + beta = total(((*!dustem_data.sed).values*dustem_sed)/(*!dustem_data.sed).sigma^2) $
  126 + / total(dustem_sed^2/(*!dustem_data.sed).sigma^2)
177 127  
178   - dustem_plot_ext,st,/print_Rv
179   -
180   - dustem_out = [ dustem_out, dustem_ext ]
181   -
182   - END
183   -
184   - 'POLEXT': BEGIN
185   -
186   - IF !dustem_show_plot NE 0 and !run_circ then begin
187   - win+=1
188   -
189   - dustem_plot_circ,st,win=win,xr=[0.1,5],aligned=st_model.align.grains.aligned,yr=[-0.03,0.03]
190   - endif
191   -
192   - IF !run_lin then begin
193   -
194   - dustem_polext = interpol(st.polext.ext_tot,st.polext.wav,(*!dustem_data.polext).wav)
195   - chi2_polext = total(((*!dustem_data.polext).values-dustem_polext)^2/(*!dustem_data.polext).sigma^2)
196   - n_polext = n_elements((*!dustem_data.polext).values)
197   - rchi2_polext = chi2_polext / n_polext
198   - wrchi2_polext = rchi2_polext * !fit_rchi2_weight.polext
199   - print,"chi2 POLEXT = ",chi2_polext," rchi2 POLEXT = ",rchi2_polext;," weighted rchi2 POLEXT=",wrchi2_polext
200   - ;print,"chi2 per wl: ",((*!dustem_data.polext).values-dustem_polext)^2/(*!dustem_data.polext).sigma^2
201   -
202   - ; Fit in (pmax,lmax,K)
203   - x=st.polext.wav
204   - logx=alog(x)
205   - logy=alog(st.polext.ext_tot)
206   - ;j=where(st.polext.ext_tot gt 0 and x gt 0.2 and x lt 1)
207   -
208   - ; Whittet+1992 (Table 1)
209   - Uband = 0.36
210   - Bband = 0.43
211   - Vband = 0.55
212   - Rband = 0.63
213   - Iband = 0.78
214   - Jband = 1.21
215   - Hband = 1.64
216   - Kband = 2.04
217   - bands = [Uband,Bband,Vband,Rband,Iband,Jband,Hband]
218   - nbands = n_elements(bands)
219   - jband = indgen(nbands)*0
220   - for k=0,nbands-1 do begin
221   - xx=min(abs(st.polext.wav-bands(k)),j)
222   - jband(k) = j
223   - endfor
224   - ;j=where(st.polext.ext_tot gt 0 and x gt 0.2 and x lt 1)
225   - res = poly_fit(logx(jband),logy(jband),2)
226   -
227   - print
228   - print,"Serkowski Fit with free K over ",string(st.polext(jband).wav,format='(100F5.2)')
229   - print,"-------------------------"
230   - K = -res(2)
231   - lmax = exp(res(1)/(2*K))
232   - pmax = exp(res(0)+K*alog(lmax)^2)
233   - print,"* lmax = ",lmax," micron"
234   - print,"* K = ",K
235   - print,"* pmax = ",pmax," for NH=1e21"
236   - print
237   -
238   - ; Plot if needed
239   - IF !dustem_show_plot NE 0 THEN BEGIN
240   -
241   - ; Polarization fraction in the UV-IR
242   - win+=1
243   -
244   - dustem_plot_polar,st,aligned=st_model.align.grains.aligned,yr=[0,0.15],/ysty,xr=[0.1,4],/xsty,win=win
245   -
246   - ; Polarisation curve (Serkowski)
247   - win+=1
248   - xr=[0.1,50]
249   - yr=[5e-25,5e-23]
250   - dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,cont=cont,/donotclose,multi=1
251   - dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,cont=cont,multi=2
252   -
253   - win+=1
254   - dustem_plot_fpol,win=win,xr=[10,1000],yr=[0,1.1]
255   -
256   -
257   - ENDIF
258   -
259   - ;dustem_plot_polar,st,/print_ratio,cont=cont
260   -
261   - dustem_out = [ dustem_out, dustem_polext ]
  128 + print,"Best SED chi2 obtained if f_HI = ",alpha, " or if we multiply model by ",beta
262 129  
263   - ENDIF
  130 + ;dustem_plot_nuinu_em,st,dustem_sed,cont,/print_radiance,p_dim=p_dim
  131 +
  132 + ; Plot if needed
  133 + IF !dustem_show_plot NE 0 THEN BEGIN
  134 + window,0
  135 + dustem_plot_fit_sed,st,dustem_sed,cont,freefree,synchrotron,_extra=extra,res=p_min*(*(*!dustem_fit).param_init_values),chi2=(*!dustem_fit).chi2,rchi2=(*!dustem_fit).rchi2
  136 + ;ytit='Total intensity'
  137 + ;xtit=textoidl('\lambda (\mum)')
  138 + ;win+=1
  139 + ;xr=[10,5e3]
  140 + ;yr=[5e-33,5.e-28]
  141 + ;nodiff = 0
  142 + ;dustem_plot_nuinu_em_vg,st,dustem_sed,cont,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1,nodiff=nodiff
  143 + ;dustem_plot_nuinu_em_vg,st,dustem_sed,cont,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2,nodiff=nodiff
  144 +
  145 + ENDIF
  146 +
  147 + dustem_out = [ dustem_out, dustem_sed ]
264 148  
265 149 END
266 150  
267   - 'POLFRAC': BEGIN
  151 + 'EXT' : BEGIN
  152 +
  153 + dustem_ext = interpol(st.ext.ext_tot,st.ext.wav,(*!dustem_data.ext).wav)
  154 + chi2_ext = total(((*!dustem_data.ext).values-dustem_ext)^2/(*!dustem_data.ext).sigma^2)
  155 + n_ext = n_elements((*!dustem_data.ext).values)
  156 + rchi2_ext = chi2_ext / n_ext
  157 + wrchi2_ext = rchi2_ext * !fit_rchi2_weight.ext
  158 + print,"chi2 EXT = ",chi2_ext," rchi2 EXT = ",rchi2_ext;," weighted rchi2 EXT=",wrchi2_ext
  159 + ;print,"chi2 per wl: ",((*!dustem_data.ext).values-dustem_ext)^2/(*!dustem_data.ext).sigma^2
268 160  
269   - IF !run_lin then begin
  161 + ; Plot if needed
  162 + IF !dustem_show_plot NE 0 THEN BEGIN
270 163  
271   - dustem_polsed = dustem_compute_polsed(p_dim,st.polsed,cont=cont)
272   - dustem_sed_tmp = dustem_compute_sed(p_dim,st,cont=cont)
273   - dustem_sed = dustem_polsed
274   - ind_filt=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_filt)
275   - for i = 0, count_filt-1 do begin
276   - j=where((*!dustem_data.polsed).filt_names(ind_filt(i)) eq (*!dustem_data.sed).filt_names,count)
277   - if count eq 1 then dustem_sed(i) = dustem_sed_tmp(j(0))
278   - endfor
279   - dustem_polfrac = dustem_polsed/dustem_sed
  164 + win+=1
  165 + ;dustem_plot_ext,st,yr=[1.e-6,1.e0],/ysty,xr=[0.5,40],/xsty,win=win
  166 + xr=[0.3,40]
  167 + yr=[1.e-26,1.e-21]
  168 + dustem_plot_ext,st,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1
  169 + dustem_plot_ext,st,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2
  170 +
  171 + win+=1
  172 + ;dustem_plot_ext,st,/UV,yr=[0,2.5],/ysty,xr=[0,10],/xsty,win=win
  173 + xr=[0.2,10]
  174 + yr=[0,3e-21]
  175 + dustem_plot_ext,st,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1
  176 + dustem_plot_ext,st,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2
  177 +
  178 + ENDIF
  179 +
  180 + dustem_plot_ext,st,/print_Rv
  181 +
  182 + dustem_out = [ dustem_out, dustem_ext ]
  183 +
  184 + END
  185 +
  186 + 'POLEXT': BEGIN
  187 +
  188 + IF !dustem_show_plot NE 0 and !run_circ then begin
  189 + win+=1
  190 +
  191 + dustem_plot_circ,st,win=win,xr=[0.1,5],aligned=st_model.align.grains.aligned,yr=[-0.03,0.03]
  192 + endif
  193 +
  194 + IF !run_lin then begin
  195 +
  196 + dustem_polext = interpol(st.polext.ext_tot,st.polext.wav,(*!dustem_data.polext).wav)
  197 + chi2_polext = total(((*!dustem_data.polext).values-dustem_polext)^2/(*!dustem_data.polext).sigma^2)
  198 + n_polext = n_elements((*!dustem_data.polext).values)
  199 + rchi2_polext = chi2_polext / n_polext
  200 + wrchi2_polext = rchi2_polext * !fit_rchi2_weight.polext
  201 + print,"chi2 POLEXT = ",chi2_polext," rchi2 POLEXT = ",rchi2_polext;," weighted rchi2 POLEXT=",wrchi2_polext
  202 + ;print,"chi2 per wl: ",((*!dustem_data.polext).values-dustem_polext)^2/(*!dustem_data.polext).sigma^2
  203 +
  204 + ; Fit in (pmax,lmax,K)
  205 + x=st.polext.wav
  206 + logx=alog(x)
  207 + logy=alog(st.polext.ext_tot)
  208 + ;j=where(st.polext.ext_tot gt 0 and x gt 0.2 and x lt 1)
  209 +
  210 + ; Whittet+1992 (Table 1)
  211 + Uband = 0.36
  212 + Bband = 0.43
  213 + Vband = 0.55
  214 + Rband = 0.63
  215 + Iband = 0.78
  216 + Jband = 1.21
  217 + Hband = 1.64
  218 + Kband = 2.04
  219 + bands = [Uband,Bband,Vband,Rband,Iband,Jband,Hband]
  220 + nbands = n_elements(bands)
  221 + jband = indgen(nbands)*0
  222 + for k=0,nbands-1 do begin
  223 + xx=min(abs(st.polext.wav-bands(k)),j)
  224 + jband(k) = j
  225 + endfor
  226 + ;j=where(st.polext.ext_tot gt 0 and x gt 0.2 and x lt 1)
  227 + res = poly_fit(logx(jband),logy(jband),2)
  228 +
  229 + print
  230 + print,"Serkowski Fit with free K over ",string(st.polext(jband).wav,format='(100F5.2)')
  231 + print,"-------------------------"
  232 + K = -res(2)
  233 + lmax = exp(res(1)/(2*K))
  234 + pmax = exp(res(0)+K*alog(lmax)^2)
  235 + print,"* lmax = ",lmax," micron"
  236 + print,"* K = ",K
  237 + print,"* pmax = ",pmax," for NH=1e21"
  238 + print
  239 +
  240 + ; Plot if needed
  241 + IF !dustem_show_plot NE 0 THEN BEGIN
  242 +
  243 + ; Polarization fraction in the UV-IR
  244 + win+=1
  245 +
  246 + dustem_plot_polar,st,aligned=st_model.align.grains.aligned,yr=[0,0.15],/ysty,xr=[0.1,4],/xsty,win=win
  247 +
  248 + ; Polarisation curve (Serkowski)
  249 + win+=1
  250 + xr=[0.1,50]
  251 + yr=[5e-25,5e-23]
  252 + dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,/donotclose,multi=1
  253 + dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,multi=2
  254 +
  255 + win+=1
  256 + dustem_plot_fpol,win=win,xr=[10,1000],yr=[0,1.1]
  257 +
  258 +
  259 + ENDIF
280 260  
281   - ; We normalize at 353GHz
282   - ;x=min(((*!dustem_data.polfrac).wav-850)^2,j353)
283   - ;dustem_polfrac = dustem_polfrac/dustem_polfrac(j353)*(*!dustem_data.polfrac).values(j353)
  261 + ;dustem_plot_polar,st,/print_ratio,cont=cont
  262 +
  263 + dustem_out = [ dustem_out, dustem_polext ]
284 264  
285   - chi2_polfrac = total(((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2)
286   - n_polfrac = n_elements((*!dustem_data.polfrac).values)
287   - rchi2_polfrac = chi2_polfrac / n_polfrac
288   - wrchi2_polfrac = rchi2_polfrac * !fit_rchi2_weight.polfrac
289   - print,"chi2 POLFRAC = ",chi2_polfrac," rchi2 POLFRAC = ",rchi2_polfrac;," weighted rchi2 POLFRAC=",wrchi2_polfrac
290   - ;print,"chi2 per wl: ",((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2
  265 + ENDIF
291 266  
292   - dustem_out = [ dustem_out, dustem_polfrac ]
  267 + END
293 268  
294   - ; Plot if needed
295   - IF !dustem_show_plot NE 0 THEN BEGIN
  269 + 'POLFRAC': BEGIN
296 270  
297   - ; Also show prediction for polarized P/I
298   - win+=1
299   - dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/Pfrac,win=win,cont=cont,yr=[0,0.25],/ysty,xr=[10,5d3]
  271 + IF !run_lin then begin
300 272  
301   - ENDIF
  273 + dustem_polsed = dustem_compute_polsed(p_dim,st.polsed,cont=cont,freefree=freefree,synchrotron=synchrotron)
  274 + dustem_sed_tmp = dustem_compute_sed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
  275 + dustem_sed = dustem_polsed
  276 + ind_filt=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_filt)
  277 + for i = 0, count_filt-1 do begin
  278 + j=where((*!dustem_data.polsed).filt_names(ind_filt(i)) eq (*!dustem_data.sed).filt_names,count)
  279 + if count eq 1 then dustem_sed(i) = dustem_sed_tmp(j(0))
  280 + endfor
  281 + dustem_polfrac = dustem_polsed/dustem_sed
302 282  
303   - ENDIF
  283 + ; We normalize at 353GHz
  284 + ;x=min(((*!dustem_data.polfrac).wav-850)^2,j353)
  285 + ;dustem_polfrac = dustem_polfrac/dustem_polfrac(j353)*(*!dustem_data.polfrac).values(j353)
304 286  
305   - END
  287 + chi2_polfrac = total(((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2)
  288 + n_polfrac = n_elements((*!dustem_data.polfrac).values)
  289 + rchi2_polfrac = chi2_polfrac / n_polfrac
  290 + wrchi2_polfrac = rchi2_polfrac * !fit_rchi2_weight.polfrac
  291 + print,"chi2 POLFRAC = ",chi2_polfrac," rchi2 POLFRAC = ",rchi2_polfrac;," weighted rchi2 POLFRAC=",wrchi2_polfrac
  292 + ;print,"chi2 per wl: ",((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2
306 293  
307   - 'POLSED': BEGIN
  294 + dustem_out = [ dustem_out, dustem_polfrac ]
308 295  
309   - IF !run_lin then begin
  296 + ; Plot if needed
  297 + IF !dustem_show_plot NE 0 THEN BEGIN
310 298  
311   - dustem_polsed = dustem_compute_polsed(p_dim,st.polsed,cont=cont)
  299 + ; Also show prediction for polarized P/I
  300 + win+=1
  301 + dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/Pfrac,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,yr=[0,0.25],/ysty,xr=[10,5d3]
312 302  
313   - ; We normalize to 353
314   - ;x=min(((*!dustem_data.polsed).wav-850)^2,j353)
315   - ;dustem_polsed = dustem_polsed/dustem_polsed(j353)*(*!dustem_data.polsed).values(j353)
  303 + ENDIF
316 304  
317   - chi2_polsed = total(((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2)
318   - n_polsed = n_elements((*!dustem_data.polsed).values)
319   - rchi2_polsed = chi2_polsed / n_polsed
320   - wrchi2_polsed = rchi2_polsed * !fit_rchi2_weight.polsed
321   - print,"chi2 POLSED = ",chi2_polsed," rchi2 POLSED = ",rchi2_polsed;," weighted rchi2 POLSED=",wrchi2_polsed
322   - print,"chi2 per wl: ",((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2
323   - ;print,"SNR per wl: ",(*!dustem_data.polsed).values/(*!dustem_data.polsed).sigma
  305 + ENDIF
324 306  
325   - dustem_out = [ dustem_out, dustem_polsed ]
  307 + END
326 308  
327   - ; Plot if needed
328   - IF !dustem_show_plot NE 0 THEN BEGIN
  309 + 'POLSED': BEGIN
  310 +;stop
  311 + IF !run_lin THEN BEGIN
  312 +
  313 + ;dustem_polsed = dustem_compute_polsed(p_dim,st.polsed,cont=cont,freefree=freefree,synchrotron=synchrotron)
  314 + dustem_polsed = dustem_compute_polsed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
  315 + ; We normalize to 353
  316 + ;x=min(((*!dustem_data.polsed).wav-850)^2,j353)
  317 + ;dustem_polsed = dustem_polsed/dustem_polsed(j353)*(*!dustem_data.polsed).values(j353)
  318 + ;try=((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2
  319 + chi2_polsed = total(((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2)
  320 + n_polsed = n_elements((*!dustem_data.polsed).values)
  321 + rchi2_polsed = chi2_polsed / n_polsed
  322 + wrchi2_polsed = rchi2_polsed * !fit_rchi2_weight.polsed
  323 + message,"chi2 POLSED = "+strtrim(chi2_polsed,2)+" rchi2 POLSED = "+strtrim(rchi2_polsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
  324 + print,"chi2 per wl: ",((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2
  325 + ;print,"SNR per wl: ",(*!dustem_data.polsed).values/(*!dustem_data.polsed).sigma
  326 +
  327 + dustem_out = [ dustem_out, dustem_polsed ]
  328 +;stop
  329 + ; Plot if needed
  330 + IF !dustem_show_plot NE 0 THEN BEGIN
329 331  
330   - win+=1
331   - xr=[10,5e3]
332   - yr=[1e-34,1e-28]
333   - dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,xr=xr,yr=yr,/xsty,/donotclose,multi=1
334   - dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,xr=xr,yr=yr,/xsty,multi=2
  332 + win+=1
  333 + xr=[10,5e3]
  334 + yr=[1e-34,1e-28]
  335 + dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,xr=xr,yr=yr,/xsty,/donotclose,multi=1
  336 + dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,xr=xr,yr=yr,/xsty,multi=2
335 337  
336   - ENDIF
  338 + ENDIF
337 339  
338   - ENDIF
  340 + ENDIF
339 341  
340 342 END
341 343  
342   - ELSE : stop
  344 + ELSE : stop
343 345  
344 346 ENDCASE
345 347  
346   -endfor
  348 +ENDFOR
347 349  
348   -; 4 bands fit : IRAS 100 and 857 down to 353 GHz
349   -;fit_sed,"~/tmp/out/SED.RES",['IRAS','HFI'],iband=[0,0,0,1,1,1,1,1,1,0],beta=fixbeta,/noplot
350 350 DOF = n_elements(p_min)
351 351 total_chi2 = 0
352 352 total_rchi2 = 0
... ...
src/idl/dustem_plot_fit_sed.pro
  Diff comments   View file @b5ccb70
1   -PRO dustem_plot_fit_sed,st,sed,cont, $
  1 +PRO dustem_plot_fit_sed,st,sed,cont,freefree,synchrotron, $
2 2 res=res,errors=errors,chi2=chi2,rchi2=rchi2, $
3 3 no_spec_error=no_spec_error, $
4 4 help=help,win=win,ps=ps,_extra=_extra
... ... @@ -83,6 +83,7 @@ IF not keyword_set(line_tot) THEN BEGIN
83 83 ENDIF ELSE BEGIN
84 84 use_line_tot=line_tot
85 85 ENDELSE
  86 +;color and linestyle for continuum
86 87 IF not keyword_set(col_cont) THEN BEGIN
87 88 ;use_col_cont=30
88 89 use_col_cont='dark blue'
... ... @@ -94,11 +95,35 @@ IF not keyword_set(line_cont) THEN BEGIN
94 95 ENDIF ELSE BEGIN
95 96 use_line_cont=line_cont
96 97 ENDELSE
97   -
  98 +;color and linestyle for freefree
  99 +IF not keyword_set(col_freefree) THEN BEGIN
  100 + ;use_col_cont=30
  101 + use_col_freefree='dark red'
  102 +ENDIF ELSE BEGIN
  103 + use_col_freefree=col_freefree
  104 +ENDELSE
  105 +IF not keyword_set(line_freefree) THEN BEGIN
  106 + use_line_freefree=0
  107 +ENDIF ELSE BEGIN
  108 + use_line_freefree=line_freefree
  109 +ENDELSE
  110 +IF not keyword_set(col_synchrotron) THEN BEGIN
  111 + ;use_col_cont=30
  112 + use_col_synchrotron='dark red'
  113 +ENDIF ELSE BEGIN
  114 + use_col_synchrotron=col_synchrotron
  115 +ENDELSE
  116 +IF not keyword_set(line_synchrotron) THEN BEGIN
  117 + use_line_synchrotron=0
  118 +ENDIF ELSE BEGIN
  119 + use_line_synchrotron=line_synchrotron
  120 +ENDELSE
98 121  
99 122 ;spec = st.sed.em_tot * fact + cont
100 123 spec = st.sed.em_tot * fact
101 124 IF !dustem_idl_continuum THEN spec=spec+cont
  125 +IF !dustem_idl_freefree THEN spec=spec+freefree
  126 +IF !dustem_idl_synchrotron THEN spec=spec+synchrotron
102 127  
103 128 ;stop
104 129  
... ... @@ -146,6 +171,7 @@ IF count_spec NE 0 THEN BEGIN
146 171 ENDIF
147 172 ENDIF
148 173 IF count_filt NE 0 THEN BEGIN
  174 + ;stop
149 175 xx=((*!dustem_data.sed).wav)[ind_filt]
150 176 yy=((*!dustem_data.sed).values)[ind_filt]/norm[ind_filt]
151 177 rms=3.*((*!dustem_data.sed).sigma)[ind_filt]/2./norm[ind_filt]
... ... @@ -181,15 +207,30 @@ ENDFOR
181 207 IF !dustem_idl_continuum THEN BEGIN
182 208 cgoplot,st.sed.wav,cont,color=use_col_cont,linestyle=use_line_cont
183 209 ENDIF
  210 +IF !dustem_idl_freefree THEN BEGIN
  211 + cgoplot,st.sed.wav,freefree,color=use_col_freefree,linestyle=use_line_freefree
  212 +ENDIF
  213 +IF !dustem_idl_synchrotron THEN BEGIN
  214 + cgoplot,st.sed.wav,synchrotron,color=use_col_synchrotron,linestyle=use_line_synchrotron
  215 +ENDIF
184 216 cgoplot,st.sed.wav,spec/norm,color=use_col_tot,linestyle=use_line_tot
185 217  
186 218 ;==== print the legend
187 219 frmt0='(A20)'
188 220 frmt1='(1E10.2)'
189 221 frmt2='(F7.2)'
  222 +;WE have to find a way to pass this through keywords ...
  223 +;use_legend_xpos=30.
  224 +;use_legend_ypos=-1
  225 +use_legend_xpos=2.e3
  226 +use_legend_ypos=2.5
  227 +IF keyword_set(legend_xpos) THEN use_legend_xpos=legend_xpos
  228 +IF keyword_set(legend_ypos) THEN use_legend_ypos=legend_ypos
  229 +
190 230 IF keyword_set(res) THEN BEGIN
191 231 offset=0.2
192   - xy_xpos=2.
  232 + ;xy_xpos=30.
  233 + ;ystart=-1
193 234 Npar=n_elements(res)
194 235 FOR i=0L,Npar-1 DO BEGIN
195 236 ; sstr=(str_sep((*!PDO)(i),'(*!dustem_params).'))
... ... @@ -200,14 +241,14 @@ IF keyword_set(res) THEN BEGIN
200 241 IF keyword_set(errors) THEN BEGIN
201 242 str=str+textoidl('\pm')+string(errors(i),format=frmt1)
202 243 ENDIF
203   - xyouts,xy_xpos,10^(1.-(i+1)*offset),str,color=0
  244 + xyouts,use_legend_xpos,10^(use_legend_ypos-(i+1)*offset),str,color=0
204 245 ENDFOR
205 246 ENDIF
206 247 IF keyword_set(chi2) THEN BEGIN
207   - xyouts,xy_xpos,10^(1.-(Npar+1)*offset),string('chi2=',format=frmt0)+string(chi2,format=frmt2),color=0
  248 + xyouts,use_legend_xpos,10^(use_legend_ypos-(Npar+1)*offset),string('chi2=',format=frmt0)+string(chi2,format=frmt2),color=0
208 249 ENDIF
209 250 IF keyword_set(rchi2) THEN BEGIN
210   - xyouts,xy_xpos,10^(1.-(Npar+2)*offset),string('red. chi2=',format=frmt0)+string(rchi2,format=frmt2),color=0
  251 + xyouts,use_legend_xpos,10^(use_legend_ypos-(Npar+2)*offset),string('red. chi2=',format=frmt0)+string(rchi2,format=frmt2),color=0
211 252 ENDIF
212 253  
213 254 IF keyword_set(ps) THEN BEGIN
... ...
src/idl/dustem_plot_polar.pro
  Diff comments   View file @b5ccb70
1 1 PRO dustem_plot_polar,st,ps=ps,_Extra=extra,help=help,UV=UV,SED=SED,Pfrac=Pfrac,print_ratio=print_ratio,win=win,cont=cont,xr=xr,yr=yr,donotclose=donotclose,aligned=aligned,noerrbars=noerrbars,multi=multi,almabands=almabands,nsmooth=nsmooth
2 2  
3   -if not keyword_set(nsmooth) then nsmooth = 1
4   -
5   -nsmooth = 4
6   -almabands = 1.
7 3 ;+
8 4 ; NAME:
9 5 ; dustem_plot_polar
... ... @@ -52,6 +48,11 @@ IF keyword_set(help) THEN BEGIN
52 48 goto,the_end
53 49 ENDIF
54 50  
  51 +if not keyword_set(nsmooth) then nsmooth = 1
  52 +
  53 +nsmooth = 4
  54 +almabands = 1.
  55 +
55 56 if not keyword_set(multi) then multi=0
56 57  
57 58 if multi ne 2 then begin
... ... @@ -64,10 +65,14 @@ if multi ne 2 then begin
64 65 endelse
65 66 endif
66 67  
  68 +;stop
  69 +
67 70 use_col_sed_filt = 70
68 71 use_col_data_filt = 250
69 72  
70   -Ngrains=n_tags(st.sed)-2
  73 +;Ngrains=n_tags(st.sed)-2
  74 +Ngrains=(*!dustem_params).Ngrains
  75 +
71 76 colors=dustem_grains_colors(Ngrains,ls,ps=ps)
72 77  
73 78 xtit=textoidl('\lambda (\mum)')
... ... @@ -82,117 +87,117 @@ xtit=textoidl('\lambda (\mum)')
82 87 ; Are observational data defined ?
83 88 defsysv,'!dustem_data',exists=data_exists
84 89  
85   -if ptr_valid(!dustem_data.polext) then begin
  90 +IF ptr_valid(!dustem_data.polext) THEN BEGIN
86 91 if keyword_set(print_ratio) then begin
87 92  
88   - ; Correction de couleur
89   - cc = 1.11
90   -
91   - Bband = 0.44
92   - Vband = 0.55
93   - Rband = 0.66
94   - Iband = 0.8
95   - Jband = 1.25
96   - Hband = 1.6
97   - Kband = 2.2
98   - ; pV and tauV are In units of 1e-21 cm2/H
99   - pV = interpol(st.polext.ext_tot,st.polext.wav,Vband)
100   - pR = interpol(st.polext.ext_tot,st.polext.wav,Rband)
101   - pI = interpol(st.polext.ext_tot,st.polext.wav,Iband)
102   - pJ = interpol(st.polext.ext_tot,st.polext.wav,Jband)
103   - pH = interpol(st.polext.ext_tot,st.polext.wav,Hband)
104   - pK = interpol(st.polext.ext_tot,st.polext.wav,Kband)
105   - j=where(st.polext.wav gt 5 and st.polext.wav lt 20)
106   - pSil=max(st.polext(j).ext_tot,jmax)
107   - Silband = st.polext(j(jmax)).wav
108   - j=where(st.polext.wav gt 0.4 and st.polext.wav lt 0.75)
109   - pmeanoptical = mean(st.polext(j).ext_tot)
110   -
111   - pV_data = interpol((*!dustem_data.polext).values,(*!dustem_data.polext).wav,Vband)
112   - tauV = interpol(st.ext.ext_tot,st.ext.wav,Vband)
113   - tauV_data = interpol((*!dustem_data.ext).values,(*!dustem_data.ext).wav,Vband)
114   - tauB = interpol(st.ext.ext_tot,st.ext.wav,Bband)
115   - tauB_data = interpol((*!dustem_data.ext).values,(*!dustem_data.ext).wav,Bband)
116   - print
117   -
118   -
119   - ; I353 and P353 are in fact nu_Pnu, in unit of erg/s/H, not erg/s/sr/H
120   - ; 1 W/m2/sr/1d20H = 1e-17 erg/s/sr/H
121   -endif
122   -if ptr_valid(!dustem_data.polsed) then begin
123   -
124   - dustem_polsed = dustem_compute_polsed(p_dim,st.polsed,cont=cont)
125   - x=min(((*!dustem_data.polsed).wav-850)^2,jmin)
126   - P353 = dustem_polsed(jmin)
127   - P353_data = (*!dustem_data.polsed).values(jmin)
128   -
129   - dustem_sed = dustem_compute_sed(p_dim,st,cont=cont)
130   - x=min(((*!dustem_data.sed).wav-850)^2,jmin)
131   - I353 = dustem_sed(jmin)
132   - I353_data = (*!dustem_data.sed).values(jmin)
133   -
134   - ; Data is Inu, not nInu, in MJy/sr/1e20
135   - ; Convert data, which is nu_Inu SED in MJy/sr/1d20H into nu_Inu SED model in erg/s/H
136   - fact = (4*!pi) * 1e-17 / 1e20 * 353e9
137   -
138   - print
139   - print,"Observational constraint"
140   - print,"-------------------------"
141   - print,"I353/NH data ",I353_data," MJy/sr for NH=1e20"
142   - print,"P353/NH data ",P353_data," MJy/sr for NH=1e20"
143   - print,"P353/I353 data ",P353_data/I353_data*100," %"
144   - print
145   - print,"Model results"
146   - print,"-------------"
147   - print,"I353/NH model ",I353," MJy/sr for NH=1e20"
148   - print,"P353/NH model ",P353," MJy/sr for NH=1e20"
149   - print,"P353/I353 model ",P353/I353*100," %"
150   -
151   - if keyword_set(print_ratio) then begin
152   - print
153   - print,"Observational constraint"
154   - print,"-------------------------"
155   -
156   - print,"tau_V data ",tauV_data
157   - print,"(p/tau)_V data ",pV_data/tauV_data*100," %"
158   - print,"p_V/A_V data ",pV_data/(tauV_data*1.086)*100," %"
159   - print,"with p_V data ",pV_data," %"
160   - print,"Rs/v data ",(P353_data/I353_data)/(pv_data/tauv_data)
161   - print,"RP/p data ",(P353_data/1e20) / (pv_data*1e-21)," MJy/sr"
162   - print,"I353/Av data ",(I353_data/1e20) / (1.086*tauv_data*1e-21)," MJy/sr"
163   - print
164   - print,"(p/tau)_V model ",pV/tauV*100," %"
165   - print,"p_V/A_V model ",pV/(tauV*1.086)*100," %"
166   - print,"p_V/E(B-V) model ",pV/((tauB-tauV)*1.086)*100," %"
167   - print,"with p_V model ",pV," %"
168   - print,"Rs/v model ",P353/I353/(pV/tauV)
169   - print,"Rs/EBV model ",P353/I353/(pV/(tauB-tauV))
170   - print,"RP/p model ",(P353/1e20) / (pv*1e-21)," MJy/sr"
171   - print,"- R band = 0.66 um ",(P353/1e20) / (pR*1e-21)," MJy/sr"
172   - print,"- I band = 0.8 um ",(P353/1e20) / (pI*1e-21)," MJy/sr"
173   - print,"- J band = 1.25 um ",(P353/1e20) / (pJ*1e-21)," MJy/sr"
174   - print,"- H band = 1.6 um ",(P353/1e20) / (pH*1e-21)," MJy/sr"
175   - print,"- K band = 2.2 um ",(P353/1e20) / (pK*1e-21)," MJy/sr"
176   - print,"- "+string(Silband,format='(F5.2)')+" um band ",(P353/1e20) / (pSil*1e-21)," MJy/sr"
177   - print,"- mean optical ",(P353/1e20) / (pmeanoptical*1e-21)," MJy/sr"
178   - print,"I353/Av model ",(I353/1e20) / (1.086*tauv*1e-21)," MJy/sr"
179   - print
180   - ;print,"Mixt results"
181   - ;print,"-------------"
182   - ;print,"(pV_model/tauV_data) ",pV/tauV_data*100," %"
183   - ;print,"I353_model/AV_data ",(I353/1e20) / (1.086*tauv_data*1e-21)," MJy/sr"
184   - ;print,"P353model/Idata ",P353/I353_data*100," %"
185   - ;print,"Rs/v P353model/I353data/(pVmodel/tauVdata) ",P353/I353_data/(pV/tauV_data)
186   - ;print,"RP/p P353model/pVdata ",(P353/1e20)/(pV_data*1e-21)," MJy/sr"
187   - ;print
188   - endif
  93 + ; Correction de couleur
  94 + cc = 1.11
  95 +
  96 + Bband = 0.44
  97 + Vband = 0.55
  98 + Rband = 0.66
  99 + Iband = 0.8
  100 + Jband = 1.25
  101 + Hband = 1.6
  102 + Kband = 2.2
  103 + ; pV and tauV are In units of 1e-21 cm2/H
  104 + pV = interpol(st.polext.ext_tot,st.polext.wav,Vband)
  105 + pR = interpol(st.polext.ext_tot,st.polext.wav,Rband)
  106 + pI = interpol(st.polext.ext_tot,st.polext.wav,Iband)
  107 + pJ = interpol(st.polext.ext_tot,st.polext.wav,Jband)
  108 + pH = interpol(st.polext.ext_tot,st.polext.wav,Hband)
  109 + pK = interpol(st.polext.ext_tot,st.polext.wav,Kband)
  110 + j=where(st.polext.wav gt 5 and st.polext.wav lt 20)
  111 + pSil=max(st.polext(j).ext_tot,jmax)
  112 + Silband = st.polext(j(jmax)).wav
  113 + j=where(st.polext.wav gt 0.4 and st.polext.wav lt 0.75)
  114 + pmeanoptical = mean(st.polext(j).ext_tot)
  115 +
  116 + pV_data = interpol((*!dustem_data.polext).values,(*!dustem_data.polext).wav,Vband)
  117 + tauV = interpol(st.ext.ext_tot,st.ext.wav,Vband)
  118 + tauV_data = interpol((*!dustem_data.ext).values,(*!dustem_data.ext).wav,Vband)
  119 + tauB = interpol(st.ext.ext_tot,st.ext.wav,Bband)
  120 + tauB_data = interpol((*!dustem_data.ext).values,(*!dustem_data.ext).wav,Bband)
  121 + print
  122 +
  123 +
  124 + ; I353 and P353 are in fact nu_Pnu, in unit of erg/s/H, not erg/s/sr/H
  125 + ; 1 W/m2/sr/1d20H = 1e-17 erg/s/sr/H
  126 + endif
  127 + IF ptr_valid(!dustem_data.polsed) THEN BEGIN
  128 +
  129 + dustem_polsed = dustem_compute_polsed(p_dim,st.polsed,cont=cont)
  130 + x=min(((*!dustem_data.polsed).wav-850)^2,jmin)
  131 + P353 = dustem_polsed(jmin)
  132 + P353_data = (*!dustem_data.polsed).values(jmin)
  133 +
  134 + dustem_sed = dustem_compute_sed(p_dim,st,cont=cont)
  135 + x=min(((*!dustem_data.sed).wav-850)^2,jmin)
  136 + I353 = dustem_sed(jmin)
  137 + I353_data = (*!dustem_data.sed).values(jmin)
  138 +
  139 + ; Data is Inu, not nInu, in MJy/sr/1e20
  140 + ; Convert data, which is nu_Inu SED in MJy/sr/1d20H into nu_Inu SED model in erg/s/H
  141 + fact = (4*!pi) * 1e-17 / 1e20 * 353e9
  142 +
  143 + print
  144 + print,"Observational constraint"
  145 + print,"-------------------------"
  146 + print,"I353/NH data ",I353_data," MJy/sr for NH=1e20"
  147 + print,"P353/NH data ",P353_data," MJy/sr for NH=1e20"
  148 + print,"P353/I353 data ",P353_data/I353_data*100," %"
  149 + print
  150 + print,"Model results"
  151 + print,"-------------"
  152 + print,"I353/NH model ",I353," MJy/sr for NH=1e20"
  153 + print,"P353/NH model ",P353," MJy/sr for NH=1e20"
  154 + print,"P353/I353 model ",P353/I353*100," %"
  155 +
  156 + IF keyword_set(print_ratio) THEN BEGIN
  157 + print
  158 + print,"Observational constraint"
  159 + print,"-------------------------"
  160 +
  161 + print,"tau_V data ",tauV_data
  162 + print,"(p/tau)_V data ",pV_data/tauV_data*100," %"
  163 + print,"p_V/A_V data ",pV_data/(tauV_data*1.086)*100," %"
  164 + print,"with p_V data ",pV_data," %"
  165 + print,"Rs/v data ",(P353_data/I353_data)/(pv_data/tauv_data)
  166 + print,"RP/p data ",(P353_data/1e20) / (pv_data*1e-21)," MJy/sr"
  167 + print,"I353/Av data ",(I353_data/1e20) / (1.086*tauv_data*1e-21)," MJy/sr"
  168 + print
  169 + print,"(p/tau)_V model ",pV/tauV*100," %"
  170 + print,"p_V/A_V model ",pV/(tauV*1.086)*100," %"
  171 + print,"p_V/E(B-V) model ",pV/((tauB-tauV)*1.086)*100," %"
  172 + print,"with p_V model ",pV," %"
  173 + print,"Rs/v model ",P353/I353/(pV/tauV)
  174 + print,"Rs/EBV model ",P353/I353/(pV/(tauB-tauV))
  175 + print,"RP/p model ",(P353/1e20) / (pv*1e-21)," MJy/sr"
  176 + print,"- R band = 0.66 um ",(P353/1e20) / (pR*1e-21)," MJy/sr"
  177 + print,"- I band = 0.8 um ",(P353/1e20) / (pI*1e-21)," MJy/sr"
  178 + print,"- J band = 1.25 um ",(P353/1e20) / (pJ*1e-21)," MJy/sr"
  179 + print,"- H band = 1.6 um ",(P353/1e20) / (pH*1e-21)," MJy/sr"
  180 + print,"- K band = 2.2 um ",(P353/1e20) / (pK*1e-21)," MJy/sr"
  181 + print,"- "+string(Silband,format='(F5.2)')+" um band ",(P353/1e20) / (pSil*1e-21)," MJy/sr"
  182 + print,"- mean optical ",(P353/1e20) / (pmeanoptical*1e-21)," MJy/sr"
  183 + print,"I353/Av model ",(I353/1e20) / (1.086*tauv*1e-21)," MJy/sr"
  184 + print
  185 + ;print,"Mixt results"
  186 + ;print,"-------------"
  187 + ;print,"(pV_model/tauV_data) ",pV/tauV_data*100," %"
  188 + ;print,"I353_model/AV_data ",(I353/1e20) / (1.086*tauv_data*1e-21)," MJy/sr"
  189 + ;print,"P353model/Idata ",P353/I353_data*100," %"
  190 + ;print,"Rs/v P353model/I353data/(pVmodel/tauVdata) ",P353/I353_data/(pV/tauV_data)
  191 + ;print,"RP/p P353model/pVdata ",(P353/1e20)/(pV_data*1e-21)," MJy/sr"
  192 + ;print
  193 + ENDIF
189 194 endif
190 195  
191 196 endif else if keyword_set(UV) then begin
192 197  
193 198 if not keyword_set(ps) then tit='Dustem polarization cross-section'
194 199 ytit=textoidl('\sigma_{pol}/N_H (cm^2/H)')
195   -if ptr_valid(!dustem_data.polext) then begin
  200 + if ptr_valid(!dustem_data.polext) then begin
196 201 if !dustem_show_plot EQ 2 or multi eq 2 then begin
197 202 ;normpol = polext_tot
198 203 normpol = st.polext.ext_tot
... ... @@ -207,7 +212,7 @@ if ptr_valid(!dustem_data.polext) then begin
207 212  
208 213 if multi eq 0 then plot_oo,st.polext.wav,st.polext.ext_tot/normpol,/nodata,_Extra=extra,xtit=xtit,ytit=ytit,tit=tit,xr=xr,yr=yr,ylog=ylog
209 214 if multi eq 1 then plot_oo,st.polext.wav,st.polext.ext_tot/normpol,/nodata,_Extra=extra,xtit='',ytit=ytit,tit=tit,xr=xr,yr=yr,ylog=ylog,position=[0.17,0.35,0.95,0.95],xtickformat='(A1)'
210   - if multi eq 2 then plot_oo,st.polext.wav,st.polext.ext_tot/normpol,/nodata,_Extra=extra,xtit=xtit,ytit='Normalized',tit='',xr=xr,yr=[-1,2],ylog=ylog,position=[0.17,0.14,0.95,0.35],/noerase,yticks=3,ymino=5,xticklen=0.1
  215 + if multi eq 2 then plot_oo,st.polext.wav,st.polext.ext_tot/normpol,/nodata,_Extra=extra,xtit=xtit,ytit='Normalized',tit='',xr=xr,yr=[-1,2],ylog=ylog,position=[0.17,0.14,0.95,0.35],/noerase,yticks=3,ymino=5,xticklen=0.1
211 216  
212 217  
213 218 FOR i=0L,Ngrains-1 DO if aligned(i) then oplot,st.polext.wav,(st.polext.abs_grain(i)+st.polext.sca_grain(i))/normpol,color=colors(i+1),psym=psym,line=ls(i+1)
... ... @@ -215,11 +220,11 @@ if ptr_valid(!dustem_data.polext) then begin
215 220 if total(aligned gt 0) then oplot,st.polext.wav,st.polext.ext_tot/normpol,psym=psym
216 221  
217 222 if data_exists then if ptr_valid(!dustem_data.polext) then begin
218   - plotsym,0,/FILL
  223 + plotsym,0,/FILL
219 224 oplot, (*!dustem_data.polext).wav,(*!dustem_data.polext).values/dustem_polext,ps=4
220 225 if not keyword_set(noerrbars) then err_bar,(*!dustem_data.polext).wav,(*!dustem_data.polext).values/dustem_polext,yrms=(*!dustem_data.polext).sigma/dustem_polext
221 226 endif
222   - endif
  227 + endif
223 228 ; ; Fit in (pmax,lmax,K)
224 229 ; x=st.polext.wav
225 230 ; logx=alog(x)
... ... @@ -255,10 +260,10 @@ if ptr_valid(!dustem_data.polext) then begin
255 260 ; ;oplot,x(j),exp(res(0)+logx(j)*res(1)-K*logx(j)^2),col=150
256 261  
257 262 endif else if ptr_valid(!dustem_data.polsed) then begin
258   -
259   - ;device, filename = 'plot008.eps', /color, /encapsulated
260   - if not keyword_set(ps) then tit='Polarization intensity in emission'
261   - ytit='Polarization intensity'
  263 +;stop
  264 + ;device, filename = 'plot008.eps', /color, /encapsulated
  265 + if not keyword_set(ps) then tit='Polarization intensity in emission'
  266 + ytit='Polarization intensity'
262 267 ytit=textoidl('\nuP_\nu/N_H (W/m^2/sr/H)')
263 268  
264 269 ; Convert from erg/s/H into W/m2/sr/H
... ... @@ -272,48 +277,50 @@ endif else if ptr_valid(!dustem_data.polsed) then begin
272 277 ylog = 0
273 278 endif else begin
274 279 normpol = st.polsed.(Ngrains+1)*0. + 1.
275   - dustem_polsed = ((*!dustem_data.polsed).values)*0. + 1
  280 + dustem_polsed = ((*!dustem_data.polsed).values)*0. + 1
276 281 ylog=1
277 282 if not keyword_set(yr) then yr=[1e-34,1e-30]
278 283 endelse
279 284  
280   - if multi eq 0 then plot_oo,st.polsed.wav,/nodata,st.polsed.(Ngrains+1)*fact/normpol,xtit=xtit,ytit=ytit,tit=tit,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs
281   - if multi eq 1 then plot_oo,st.polsed.wav,/nodata,st.polsed.(Ngrains+1)*fact/normpol,xtit='',ytit=ytit,tit=tit,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.95,0.95],xtickformat='(A1)'
282   - if multi eq 2 then plot_oo,st.polsed.wav,/nodata,st.polsed.(Ngrains+1)*fact/normpol,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.95,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1
  285 + ;if multi eq 0 then plot_oo,st.polsed.wav,/nodata,st.polsed.(Ngrains+1)*fact/normpol,xtit=xtit,ytit=ytit,tit=tit,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs
  286 + ;if multi eq 1 then plot_oo,st.polsed.wav,/nodata,st.polsed.(Ngrains+1)*fact/normpol,xtit='',ytit=ytit,tit=tit,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.95,0.95],xtickformat='(A1)'
  287 + ;if multi eq 2 then plot_oo,st.polsed.wav,/nodata,st.polsed.(Ngrains+1)*fact/normpol,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.95,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1
  288 + if multi eq 0 then cgplot,st.polsed.wav,/nodata,st.polsed.(Ngrains+1)*fact/normpol,xtit=xtit,ytit=ytit,tit=tit,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs
  289 + if multi eq 1 then cgplot,st.polsed.wav,/nodata,st.polsed.(Ngrains+1)*fact/normpol,xtit='',ytit=ytit,tit=tit,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.95,0.95],xtickformat='(A1)'
  290 + if multi eq 2 then cgplot,st.polsed.wav,/nodata,st.polsed.(Ngrains+1)*fact/normpol,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.95,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1
283 291  
284   - for i = 0, Ngrains-1 do begin
285   - if aligned(i) then oplot,st.polsed.wav,st.polsed.(i+1)*fact/normpol,color=colors(i+1),psym=psym,line=ls(i+1)
  292 + for i = 0, Ngrains-1 do begin
  293 + if aligned(i) then cgoplot,st.polsed.wav,st.polsed.(i+1)*fact/normpol,color=colors(i+1),psym=psym,line=ls(i+1)
286 294 endfor
287   - if (total(aligned) gt 0) then oplot,st.polsed.wav,st.polsed.(Ngrains+1)*fact/normpol
  295 + ;stop
  296 + if (total(aligned) gt 0) then cgoplot,st.polsed.wav,st.polsed.(Ngrains+1)*fact/normpol
288 297  
289   - if data_exists then if ptr_valid(!dustem_data.polsed) then begin
  298 + if data_exists then if ptr_valid(!dustem_data.polsed) then begin
290 299  
291 300 ind_filt=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_filt)
292 301 ;=== Plot the data
293 302 IF count_filt NE 0 THEN BEGIN
294 303  
295   - if !dustem_show_plot EQ 2 or multi eq 2 then begin
296   - endif else begin
  304 + if !dustem_show_plot EQ 2 or multi eq 2 then begin
  305 + endif else begin
297 306 ; Convert data which Pnu in MJy/sr/1d20 into nu*Pnu in W/m2/sr/1d20
298 307 fact = (3.e10/(1.e-4*(*!dustem_data.polsed).wav(ind_filt)))/1.d40
299   - endelse
300   -
301   - dustem_polsed = dustem_compute_polsed(p_dim,st.polsed,cont=cont)*fact
302   - if !dustem_show_plot EQ 2 or multi eq 2 then begin
  308 + endelse
  309 + dustem_polsed = dustem_compute_polsed(p_dim,st,cont=cont)*fact
  310 + if !dustem_show_plot EQ 2 or multi eq 2 then begin
303 311 normpol = dustem_polsed
304   - endif else begin
  312 + endif else begin
305 313 normpol = dustem_polsed * 0. + 1
306   - endelse
307   -
308   - plotsym,0,/FILL
309   - oplot,((*!dustem_data.polsed).wav)(ind_filt),((*!dustem_data.polsed).values)(ind_filt)*fact/normpol(ind_filt),psym=8,_extra=extra,color=use_col_sed_filt,symsize=0.7
310   - if not keyword_set(noerrbars) then err_bar,((*!dustem_data.polsed).wav)(ind_filt),((*!dustem_data.polsed).values)(ind_filt)*fact/normpol(ind_filt),yrms=3.*((*!dustem_data.polsed).sigma)(ind_filt)/2.*fact/normpol(ind_filt),color=use_col_sed_filt
  314 + endelse
  315 + plotsym,0,/FILL
  316 + cgoplot,((*!dustem_data.polsed).wav)(ind_filt),((*!dustem_data.polsed).values)(ind_filt)*fact/normpol(ind_filt),psym=8,_extra=extra,color=use_col_sed_filt,symsize=0.7
  317 + if not keyword_set(noerrbars) then err_bar,((*!dustem_data.polsed).wav)(ind_filt),((*!dustem_data.polsed).values)(ind_filt)*fact/normpol(ind_filt),yrms=3.*((*!dustem_data.polsed).sigma)(ind_filt)/2.*fact/normpol(ind_filt),color=use_col_sed_filt
311 318  
312 319 ; PLot color corrected POLSED
313   - if keyword_set(ps) then plotsym,8,thick=8 else plotsym,8
314   - if multi eq 0 then oplot,((*!dustem_data.polsed).wav)(ind_filt),dustem_polsed(ind_filt)/normpol(ind_filt),psym=8,syms=2,_extra=extra,col=use_col_data_filt
  320 + if keyword_set(ps) then plotsym,8,thick=8 else plotsym,8
  321 + if multi eq 0 then cgoplot,((*!dustem_data.polsed).wav)(ind_filt),dustem_polsed(ind_filt)/normpol(ind_filt),psym=8,syms=2,_extra=extra,col=use_col_data_filt
315 322  
316   - endif
  323 + endif
317 324 endif
318 325  
319 326 endif else if keyword_set(Pfrac) then begin
... ...
src/idl/dustem_read_align.pro
  Diff comments   View file @b5ccb70
... ... @@ -16,6 +16,7 @@ WHILE not eof(unit) DO BEGIN
16 16 ENDWHILE
17 17 close,unit
18 18 free_lun,unit
  19 +
19 20 ;== now read the file
20 21 openr,unit,file,/get_lun
21 22 ncurrent=0L
... ... @@ -38,58 +39,108 @@ anisG0 = float(xx(0))
38 39 ncurrent+=+1
39 40 nalig=nlines-ncurrent
40 41  
  42 +;stop
  43 +
41 44 one_st={aligned:0.,law:'',athresh:0.,plev:0.,pstiff:0.,TdsTg:0.,a0:0.,rvcut:0.}
42   -st=replicate(one_st,st_grains.Ngrains)
  45 +;Ngrains=st_grains.Ngrains
  46 +Ngrains=n_elements(st_grains)
43 47  
44   -st.aligned = stregex(st_grains.grains.type_keywords, 'pol', /bool)
  48 +st=replicate(one_st,Ngrains)
  49 +
  50 +;st.aligned = stregex(st_grains.grains.type_keywords, 'pol', /bool)
45 51 ;count=0
46   -FOR i=0L,st_grains.Ngrains-1 DO BEGIN
47   - IF st(i).aligned then begin
48   - ; if count EQ 0. then begin
49   - readf,unit,str,format='(A100)'
50   - ;endif
  52 +
  53 +FOR i=0L,Ngrains-1 DO BEGIN
  54 + ;stop
  55 + is_aligned = stregex(st_grains(i).type_keywords, 'pol', /bool)
  56 +
  57 + ;IF st(i).aligned then begin
  58 + IF is_aligned THEN BEGIN
  59 + ;stop
  60 + readf,unit,str,format='(A100)'
  61 + strv=str_sep(strcompress(strtrim(str,2)),' ')
  62 + Nstrv=n_elements(strv)
  63 +
  64 + ii=0L
  65 + st(i).law=strv(ii) & ii=ii+1
  66 +
  67 + ; IDG
  68 + IF stregex(st(i).law, 'idg', /bool) THEN BEGIN
  69 + st[i].TdsTg = strv[ii] & ii=ii+1
  70 + st[i].a0 = strv[ii] & ii=ii+1
  71 + st[i].rvcut = strv[ii] & ii=ii+1
  72 + ENDIF
  73 +
  74 + ; RAT
  75 + IF stregex(st(i).law, 'rat', /bool) THEN BEGIN
  76 + ENDIF
  77 +
  78 + ; PARAMETRIC
  79 + IF stregex(st(i).law, 'par', /bool) THEN BEGIN
  80 + ; Grain radius Threshold for alignment, given in microns and converted in cm
  81 + st[i].athresh = strv(ii) & ii=ii+1
  82 + st[i].pstiff = strv(ii) & ii=ii+1
  83 + st[i].plev = strv(ii) & ii=ii+1
  84 + ENDIF
  85 + st[i].aligned = stregex(st_grains(i).type_keywords, 'pol', /bool)
  86 + ;stop
  87 +
  88 + IF !run_univ eq 1 THEN BEGIN
  89 + st = replicate(st(i),Ngrains)
  90 + ;st.aligned = stregex(st_grains.grains.type_keywords, 'pol', /bool)
  91 + st.aligned = stregex(st_grains(i).type_keywords, 'pol', /bool)
  92 + break
  93 + ENDIF
  94 + ENDIF
  95 +ENDFOR
  96 +
  97 +; FOR i=0L,st_grains.Ngrains-1 DO BEGIN
  98 +; IF st(i).aligned then begin
  99 +; ; if count EQ 0. then begin
  100 +; readf,unit,str,format='(A100)'
  101 +; ;endif
51 102  
52   - strv=str_sep(strcompress(strtrim(str,2)),' ')
53   - Nstrv=n_elements(strv)
54   -
55   - ii=0L
56   - st(i).law=strv(ii) & ii=ii+1
57   -
58   - ; IDG
59   - IF stregex(st(i).law, 'idg', /bool) THEN begin
60   - ii=ii+1
61   - st(i).TdsTg = strv(ii) & ii=ii+1
62   - st(i).a0 = strv(ii) & ii=ii+1
63   - st(i).rvcut = strv(ii) & ii=ii+1
64   - endif
65   -
66   - ; RAT
67   - IF stregex(st(i).law, 'rat', /bool) THEN begin
68   - endif
69   -
70   - ; PARAMETRIC
71   - IF stregex(st(i).law, 'par', /bool) THEN begin
72   - ; Grain radius Threshold for
73   - ; alignment, given in microns and
74   - ; converted in cm
75   - st(i).athresh = strv(ii) & ii=ii+1
76   - st(i).pstiff = strv(ii) & ii=ii+1
77   - st(i).plev = strv(ii) & ii=ii+1
78   - ; count=1
79   - endif
80   -
81   - IF !run_univ eq 1 then begin
82   - st = replicate(st(i),st_grains.Ngrains)
83   - st.aligned = stregex(st_grains.grains.type_keywords, 'pol', /bool)
84   - ;stop
85   - break
86   - endif
  103 +; strv=str_sep(strcompress(strtrim(str,2)),' ')
  104 +; Nstrv=n_elements(strv)
  105 +
  106 +; ii=0L
  107 +; st(i).law=strv(ii) & ii=ii+1
  108 +
  109 +; ; IDG
  110 +; IF stregex(st(i).law, 'idg', /bool) THEN begin
  111 +; ii=ii+1
  112 +; st(i).TdsTg = strv(ii) & ii=ii+1
  113 +; st(i).a0 = strv(ii) & ii=ii+1
  114 +; st(i).rvcut = strv(ii) & ii=ii+1
  115 +; endif
  116 +
  117 +; ; RAT
  118 +; IF stregex(st(i).law, 'rat', /bool) THEN begin
  119 +; endif
  120 +
  121 +; ; PARAMETRIC
  122 +; IF stregex(st(i).law, 'par', /bool) THEN begin
  123 +; ; Grain radius Threshold for
  124 +; ; alignment, given in microns and
  125 +; ; converted in cm
  126 +; st(i).athresh = strv(ii) & ii=ii+1
  127 +; st(i).pstiff = strv(ii) & ii=ii+1
  128 +; st(i).plev = strv(ii) & ii=ii+1
  129 +; ; count=1
  130 +; endif
  131 +
  132 +; IF !run_univ eq 1 then begin
  133 +; st = replicate(st(i),st_grains.Ngrains)
  134 +; st.aligned = stregex(st_grains.grains.type_keywords, 'pol', /bool)
  135 +; ;stop
  136 +; break
  137 +; endif
87 138  
88 139  
89 140  
90   - ENDIF
  141 +; ENDIF
91 142  
92   -ENDFOR
  143 +; ENDFOR
93 144  
94 145 close,unit
95 146 free_lun,unit
... ...
src/idl/dustem_read_all_web3p8.pro
  Diff comments   View file @b5ccb70
... ... @@ -147,7 +147,8 @@ IF !run_pol THEN BEGIN
147 147 IF stregex(st_grains(i).type_keywords, 'pol', /bool) THEN BEGIN
148 148 IF not stregex(st_align.keywords, 'rrf', /bool) THEN BEGIN
149 149 FOR i_axis = 1, 2 DO BEGIN
150   - Q_file=dir_in_qabs+'Q'+strtrim(i_axis,2)+'_'+st_grains.grains(i).grain_type+'.DAT'
  150 + ;Q_file=dir_in_qabs+'Q'+strtrim(i_axis,2)+'_'+st_grains.grains(i).grain_type+'.DAT'
  151 + Q_file=dir_in_qabs+'Q'+strtrim(i_axis,2)+'_'+st_grains(i).grain_type+'.DAT'
151 152 st=dustem_read_qabs_lv(Q_file,silent=silent)
152 153 st_qpol(i_axis-1,i)=ptr_new(st)
153 154 ENDFOR
... ... @@ -155,14 +156,16 @@ IF !run_pol THEN BEGIN
155 156  
156 157 IF (st_align.anisG0 > 0) THEN BEGIN
157 158 FOR i_axis = 1, 2 DO BEGIN
158   - Q_file=dir_in_qabs+'QH'+strtrim(i_axis,2)+'_'+st_grains.grains(i).grain_type+'.DAT'
  159 + ;Q_file=dir_in_qabs+'QH'+strtrim(i_axis,2)+'_'+st_grains.grains(i).grain_type+'.DAT'
  160 + Q_file=dir_in_qabs+'QH'+strtrim(i_axis,2)+'_'+st_grains(i).grain_type+'.DAT'
159 161 st=dustem_read_qabs_lv(Q_file,silent=silent)
160 162 st_qH(i_axis-1,i)=ptr_new(st)
161 163 ENDFOR
162 164 ENDIF
163 165  
164 166 IF stregex(st_align.keywords, 'circ', /bool) THEN BEGIN
165   - Q_file=dir_in_qabs+'Qc_'+st_grains.grains(i).grain_type+'.DAT'
  167 + ;Q_file=dir_in_qabs+'Qc_'+st_grains.grains(i).grain_type+'.DAT'
  168 + Q_file=dir_in_qabs+'Qc_'+st_grains(i).grain_type+'.DAT'
166 169 st=dustem_read_qabs_lv(Q_file,silent=silent)
167 170 st_qcirc(i)=ptr_new(st)
168 171 ENDIF
... ... @@ -192,8 +195,16 @@ if !run_tls then begin
192 195 ENDIF
193 196  
194 197 ;=== This will become !dustem_params
195   -st={Ngrains:Ngrains,G0:G0,Keywords:key_str,grains:st_grains, $
196   - isrf:st_isrf,qabs:st_qabs,calor:st_calor,lambda:st_lambda,size:st_size,mix:st_mix,chrg:st_chrg,gas:st_gas,spin:st_spin,pol:st_pol,tls:st_tls}
  198 +;st={Ngrains:Ngrains,G0:G0,Keywords:key_str,grains:st_grains, $
  199 +; isrf:st_isrf,qabs:st_qabs,calor:st_calor,lambda:st_lambda,size:st_size,mix:st_mix,chrg:st_chrg,gas:st_gas,spin:st_spin,pol:st_pol,tls:st_tls}
  200 +;VG Add NH
  201 +if !run_pol then begin
  202 + st={Ngrains:Ngrains,G0:G0,Keywords:key_str,grains:st_grains,isrf:st_isrf,qabs:st_qabs,calor:st_calor,lambda:st_lambda,size:st_size,mix:st_mix,chrg:st_chrg,gas:st_gas,spin:st_spin,align:st_align,qpol:st_qpol,qcirc:st_qcirc,qH:st_qH,tls:st_tls}
  203 +endif else begin
  204 + st={Ngrains:Ngrains,G0:G0,Keywords:key_str,grains:st_grains,isrf:st_isrf,qabs:st_qabs,calor:st_calor,lambda:st_lambda,size:st_size,mix:st_mix,chrg:st_chrg,gas:st_gas,spin:st_spin,tls:st_tls}
  205 +endelse
  206 +
  207 +;stop
197 208  
198 209 the_end:
199 210  
... ...
src/idl/dustem_read_filters.pro
  Diff comments   View file @b5ccb70
... ... @@ -73,6 +73,7 @@ dir_wmap=filter_dir+'WMAP'+'/'
73 73 dir_spire=filter_dir+'SPIRE'+'/'
74 74 dir_pacs=filter_dir+'PACS'+'/'
75 75 dir_akari=filter_dir+'AKARI'+'/'
  76 +dir_spass=filter_dir+'SPASS'+'/'
76 77  
77 78 ;Read Herschel PACS filter Transmission
78 79 ;wfilt_pacs=[75,110,170.] ;microns
... ... @@ -767,6 +768,25 @@ gismo_struct={Name:'GISMO',Nbands:Nband,filter_names:filt_names,central_waveleng
767 768 use_frequencies:use_nu,use_transmissions:use_T,use_wavelengths:use_w,use_wmin:dblarr(Nband),use_wmax:dblarr(Nband),cc_method:'dustem_cc_gismo'}
768 769  
769 770  
  771 +;Read S-PASS filters
  772 +Nband=1
  773 +filt_names=['SPASS1']
  774 +wfilt_spass=dustem_filter2wav(filt_names)
  775 +wav_ptr=ptrarr(Nband)
  776 +nu_ptr=ptrarr(Nband) & T_ptr=ptrarr(Nband) & wav_ptr=ptrarr(Nband)
  777 +use_nu=ptrarr(Nband) & use_T=ptrarr(Nband) & use_w=ptrarr(Nband)
  778 +
  779 +st=read_xcat(dir_spass+'SPASS_TRANSMISSION.xcat',/silent)
  780 +
  781 +freq=st.frequency*1.e9 ;Hz
  782 +Trans=st.transmission
  783 +
  784 +T_ptr(0)=ptr_new(Trans/max(Trans)) & nu_ptr(0)=ptr_new(freq) & wav_ptr(0)=ptr_new(cmic/freq)
  785 +
  786 +spass_struct={Name:'SPASS',Nbands:Nband,filter_names:filt_names,central_wavelengths:wfilt_spass,central_frequencies:cmic/wfilt_spass, $
  787 + filter_wavelengths:wav_ptr,filter_frequencies:nu_ptr,filter_transmissions:T_ptr, $
  788 + use_frequencies:use_nu,use_transmissions:use_T,use_wavelengths:use_w,use_wmin:dblarr(Nband),use_wmax:dblarr(Nband),cc_method:'dustem_cc_spass'}
  789 +
770 790 ;== Make the filter strcuture
771 791 dm_filter_struct={IRAC:irac_struct,MIPS:mips_struct,MSX:msx_struct,IRAS:iras_struct,DIRBE:dirbe_struct, $
772 792 SPM:spm_struct, $
... ... @@ -777,7 +797,7 @@ dm_filter_struct={IRAC:irac_struct,MIPS:mips_struct,MSX:msx_struct,IRAS:iras_str
777 797 HFI:hfi_struct,LFI:lfi_struct, $
778 798 WMAP:wmap_struct, $
779 799 SPIRE:spire_struct,PACS:pacs_struct, $
780   - AKARI:akari_struct, BOLOCAM:bolocam_struct, WISE:wise_struct, LABOCA:laboca_struct, GISMO:gismo_struct}
  800 + AKARI:akari_struct, BOLOCAM:bolocam_struct, WISE:wise_struct, LABOCA:laboca_struct, GISMO:gismo_struct,SPASS:spass_struct}
781 801  
782 802 ;== If requested, plot transmissions
783 803  
... ...
src/idl/dustem_read_grain.pro
  Diff comments   View file @b5ccb70
... ... @@ -81,7 +81,8 @@ IF keyword_set(help) THEN BEGIN
81 81 goto,the_end
82 82 ENDIF
83 83  
84   -defsysv, '!run_pol', 0. ; Global flag to know if polarisation is being run
  84 +;!run_pol now defined in dustem_init.pro
  85 +;defsysv, '!run_pol', 0. ; Global flag to know if polarisation is being run
85 86 defsysv,'!run_tls',0. ; Global flag to know if TLS is being run
86 87  
87 88 CASE !dustem_which OF
... ...
src/idl/dustem_restore_system_variables.pro 0 → 100644
  Diff comments   View file @b5ccb70
... ... @@ -0,0 +1,85 @@
  1 +PRO dustem_restore_system_variables,file,help=help
  2 +
  3 +;+
  4 +; NAME:
  5 +; dustem_restore_system_variables
  6 +; PURPOSE:
  7 +; Restores Dustem system variables a file
  8 +; CATEGORY:
  9 +; Dustem
  10 +; CALLING SEQUENCE:
  11 +; dustem_restore_system_variables,file[,/help]
  12 +; INPUTS:
  13 +; file = File name
  14 +; OPTIONAL INPUT PARAMETERS:
  15 +; None
  16 +; OUTPUTS:
  17 +; None
  18 +; OPTIONAL OUTPUT PARAMETERS:
  19 +; None
  20 +; ACCEPTED KEY-WORDS:
  21 +; help = If set, print this help
  22 +; COMMON BLOCKS:
  23 +; None
  24 +; SIDE EFFECTS:
  25 +; File is restored.
  26 +; The following system variables are resstored:
  27 +; !dustem_fit
  28 +; !dustem_data
  29 +; !dustem_filters
  30 +; !run_ionfrac
  31 +; !dustem_verbose
  32 +; RESTRICTIONS:
  33 +; The dustem idl wrapper must be installed
  34 +; PROCEDURE:
  35 +; None
  36 +; EXAMPLES
  37 +;
  38 +; MODIFICATION HISTORY:
  39 +; Written by J.-Ph. Bernard
  40 +; see evolution details on the dustem cvs maintained at CESR
  41 +; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
  42 +;-
  43 +
  44 +IF keyword_set(help) THEN BEGIN
  45 + doc_library,'dustem_restore_system_variables'
  46 + goto,the_end
  47 +ENDIF
  48 +
  49 +restore,file
  50 +
  51 +defsysv,'!dustem_dat',saved_dustem_dat
  52 +defsysv,'!dustem_data',saved_dustem_data
  53 +defsysv,'!DUSTEM_DO_CC',saved_DUSTEM_DO_CC
  54 +defsysv,'!DUSTEM_F90_EXEC',saved_DUSTEM_F90_EXEC
  55 +defsysv,'!dustem_filters',saved_dustem_filters
  56 +defsysv,'!dustem_fit',saved_dustem_fit
  57 +defsysv,'!DUSTEM_F_HI',saved_DUSTEM_F_HI
  58 +defsysv,'!DUSTEM_IDL_CONTINUUM',saved_DUSTEM_IDL_CONTINUUM
  59 +defsysv,'!DUSTEM_IDL_FREEFREE',saved_DUSTEM_IDL_FREEFREE
  60 +defsysv,'!DUSTEM_IDL_SYNCHROTRON',saved_DUSTEM_IDL_SYNCHROTRON
  61 +defsysv,'!DUSTEM_INPUTS',saved_DUSTEM_INPUTS
  62 +defsysv,'!DUSTEM_INSTRUMENT_DESCRIPTION',saved_DUSTEM_INSTRUMENT_DESCRIPTION
  63 +defsysv,'!DUSTEM_NEVER_DO_CC',saved_DUSTEM_NEVER_DO_CC
  64 +defsysv,'!DUSTEM_PARAMS',saved_DUSTEM_PARAMS
  65 +defsysv,'!DUSTEM_PARINFO',saved_DUSTEM_PARINFO
  66 +defsysv,'!DUSTEM_PREVIOUS_CC',saved_DUSTEM_PREVIOUS_CC
  67 +defsysv,'!DUSTEM_RES',saved_DUSTEM_RES
  68 +defsysv,'!DUSTEM_SHOW_PLOT',saved_DUSTEM_SHOW_PLOT
  69 +defsysv,'!DUSTEM_SOFT_DIR',saved_DUSTEM_SOFT_DIR
  70 +defsysv,'!DUSTEM_VERBOSE',saved_DUSTEM_VERBOSE
  71 +defsysv,'!DUSTEM_WHICH',saved_DUSTEM_WHICH
  72 +defsysv,'!DUSTEM_WRAP_SOFT_DIR',saved_DUSTEM_WRAP_SOFT_DIR
  73 +defsysv,'!FIT_RCHI2_WEIGHT',saved_FIT_RCHI2_WEIGHT
  74 +defsysv,'!RUN_ANIS',saved_RUN_ANIS
  75 +defsysv,'!RUN_CIRC',saved_RUN_CIRC
  76 +defsysv,'!RUN_IONFRAC',saved_RUN_IONFRAC
  77 +defsysv,'!RUN_LIN',saved_RUN_LIN
  78 +defsysv,'!RUN_POL',saved_RUN_POL
  79 +defsysv,'!RUN_RRF',saved_RUN_RRF
  80 +defsysv,'!RUN_TLS',saved_RUN_TLS
  81 +defsysv,'!RUN_UNIV',saved_RUN_UNIV
  82 +
  83 +the_end:
  84 +
  85 +END
... ...
src/idl/dustem_save_sed_fit.pro
  Diff comments   View file @b5ccb70
... ... @@ -76,7 +76,8 @@ errors=*(*!dustem_fit).current_param_errors
76 76  
77 77 save,param_desc,init_value,fixed_param_desc,fixed_init_value, $
78 78 dustem_fit_chi2,dustem_fit_rchi2,dfp,func_ind, $
79   - lambda_sed,run_ionfrac, $
  79 + ;lambda_sed, $
  80 + run_ionfrac, $
80 81 dustem_data,dustem_filters,dustem_verbose,errors, $
81 82 ; dustem_fit_sed, $
82 83 file=file
... ...
src/idl/dustem_save_system_variables.pro 0 → 100644
  Diff comments   View file @b5ccb70
... ... @@ -0,0 +1,149 @@
  1 +PRO dustem_save_system_variables,file,help=help
  2 +
  3 +;+
  4 +; NAME:
  5 +; dustem_save_system_variables
  6 +; PURPOSE:
  7 +; Save Dustem system variables a file
  8 +; CATEGORY:
  9 +; Dustem
  10 +; CALLING SEQUENCE:
  11 +; dustem_save_system_variables,file[,/help]
  12 +; INPUTS:
  13 +; file = File name
  14 +; OPTIONAL INPUT PARAMETERS:
  15 +; None
  16 +; OUTPUTS:
  17 +; None
  18 +; OPTIONAL OUTPUT PARAMETERS:
  19 +; None
  20 +; ACCEPTED KEY-WORDS:
  21 +; help = If set, print this help
  22 +; COMMON BLOCKS:
  23 +; None
  24 +; SIDE EFFECTS:
  25 +; File is stored.
  26 +; The following system variables are stored:
  27 +; !dustem_fit
  28 +; !dustem_data
  29 +; !dustem_filters
  30 +; !run_ionfrac
  31 +; !dustem_verbose
  32 +; RESTRICTIONS:
  33 +; The dustem idl wrapper must be installed
  34 +; PROCEDURE:
  35 +; None
  36 +; EXAMPLES
  37 +;
  38 +; MODIFICATION HISTORY:
  39 +; Written by J.-Ph. Bernard
  40 +; see evolution details on the dustem cvs maintained at CESR
  41 +; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
  42 +;-
  43 +
  44 +; !DUSTEM_DAT = '/tmp/ramdisk8MB/'
  45 +; !DUSTEM_DATA = -> <Anonymous> Array[1]
  46 +; !DUSTEM_DO_CC = 1
  47 +; !DUSTEM_F90_EXEC = '$HOME/Soft_Librairies/dustem4.3_wk/src/dustem'
  48 +; !DUSTEM_FILTERS = <PtrHeapVar570>
  49 +; !DUSTEM_FIT = <PtrHeapVar1>
  50 +; !DUSTEM_F_HI = 1.00000
  51 +; !DUSTEM_IDL_CONTINUUM = 0.00000
  52 +; !DUSTEM_IDL_FREEFREE = 0.00000
  53 +; !DUSTEM_IDL_SYNCHROTRON = 0.00000
  54 +; !DUSTEM_INPUTS = <PtrHeapVar2>
  55 +; !DUSTEM_INSTRUMENT_DESCRIPTION = -> <Anonymous> Array[120]
  56 +; !DUSTEM_NEVER_DO_CC = 0
  57 +; !DUSTEM_PARAMS = <PtrHeapVar15>
  58 +; !DUSTEM_PARINFO = <PtrHeapVar943>
  59 +; !DUSTEM_PREVIOUS_CC = <PtrHeapVar3346>
  60 +; !DUSTEM_RES = '/tmp/ramdisk8MB/'
  61 +; !DUSTEM_SHOW_PLOT = 1
  62 +; !DUSTEM_SOFT_DIR = '$HOME/Soft_Librairies/dustem4.3_wk/'
  63 +; !DUSTEM_VERBOSE = 1
  64 +; !DUSTEM_WHICH = 'WEB3p8'
  65 +; !DUSTEM_WRAP_SOFT_DIR = '$HOME/Soft_Librairies/dustem-wrapper_idl/'
  66 +; !FIT_RCHI2_WEIGHT = -> <Anonymous> Array[1]
  67 +; !RUN_ANIS = 0.00000
  68 +; !RUN_CIRC = 0.00000
  69 +; !RUN_IONFRAC = 0.00000
  70 +; !RUN_LIN = 1.00000
  71 +; !RUN_POL = 1
  72 +; !RUN_RRF = 0.00000
  73 +; !RUN_TLS = 0.00000
  74 +; !RUN_UNIV = 0.00000
  75 +
  76 +IF keyword_set(help) THEN BEGIN
  77 + doc_library,'dustem_save_system_variables'
  78 + goto,the_end
  79 +ENDIF
  80 +
  81 +saved_dustem_dat=!dustem_dat
  82 +saved_dustem_data=!dustem_data
  83 +saved_dustem_do_cc=!DUSTEM_DO_CC
  84 +saved_dustem_f90_exec=!DUSTEM_F90_EXEC
  85 +saved_dustem_filters=!dustem_filters
  86 +saved_dustem_fit=!dustem_fit
  87 +saved_dustem_f_hi=!DUSTEM_F_HI
  88 +saved_dustem_IDL_CONTINUUM=!DUSTEM_IDL_CONTINUUM
  89 +saved_DUSTEM_IDL_FREEFREE=!DUSTEM_IDL_FREEFREE
  90 +saved_DUSTEM_IDL_SYNCHROTRON=!DUSTEM_IDL_SYNCHROTRON
  91 +saved_DUSTEM_INPUTS=!DUSTEM_INPUTS
  92 +saved_DUSTEM_INSTRUMENT_DESCRIPTION=!DUSTEM_INSTRUMENT_DESCRIPTION
  93 +saved_DUSTEM_NEVER_DO_CC=!DUSTEM_NEVER_DO_CC
  94 +saved_DUSTEM_PARAMS=!DUSTEM_PARAMS
  95 +saved_DUSTEM_PARINFO=!DUSTEM_PARINFO
  96 +saved_DUSTEM_PREVIOUS_CC=!DUSTEM_PREVIOUS_CC
  97 +saved_DUSTEM_RES=!DUSTEM_RES
  98 +saved_DUSTEM_SHOW_PLOT=!DUSTEM_SHOW_PLOT
  99 +saved_DUSTEM_SOFT_DIR=!DUSTEM_SOFT_DIR
  100 +saved_DUSTEM_VERBOSE=!DUSTEM_VERBOSE
  101 +saved_DUSTEM_WHICH=!DUSTEM_WHICH
  102 +saved_DUSTEM_WRAP_SOFT_DIR=!DUSTEM_WRAP_SOFT_DIR
  103 +saved_FIT_RCHI2_WEIGHT=!FIT_RCHI2_WEIGHT
  104 +saved_RUN_ANIS=!RUN_ANIS
  105 +saved_RUN_CIRC=!RUN_CIRC
  106 +saved_RUN_IONFRAC=!RUN_IONFRAC
  107 +saved_RUN_LIN=!RUN_LIN
  108 +saved_RUN_POL=!RUN_POL
  109 +saved_RUN_RRF=!RUN_RRF
  110 +saved_RUN_TLS=!RUN_TLS
  111 +saved_RUN_UNIV=!RUN_UNIV
  112 +
  113 +save,saved_dustem_dat,$
  114 + saved_dustem_data, $
  115 + saved_dustem_fit, $
  116 + saved_dustem_do_cc, $
  117 + saved_dustem_f90_exec, $
  118 + saved_dustem_filters, $
  119 + saved_dustem_fit, $
  120 + saved_dustem_f_hi, $
  121 + saved_dustem_IDL_CONTINUUM, $
  122 + saved_DUSTEM_IDL_FREEFREE, $
  123 + saved_DUSTEM_IDL_SYNCHROTRON, $
  124 + saved_DUSTEM_INPUTS, $
  125 + saved_DUSTEM_INSTRUMENT_DESCRIPTION, $
  126 + saved_DUSTEM_NEVER_DO_CC, $
  127 + saved_DUSTEM_PARAMS, $
  128 + saved_DUSTEM_PARINFO, $
  129 + saved_DUSTEM_PREVIOUS_CC, $
  130 + saved_DUSTEM_RES, $
  131 + saved_DUSTEM_SHOW_PLOT, $
  132 + saved_DUSTEM_SOFT_DIR, $
  133 + saved_DUSTEM_VERBOSE, $
  134 + saved_DUSTEM_WHICH, $
  135 + saved_DUSTEM_WRAP_SOFT_DIR, $
  136 + saved_FIT_RCHI2_WEIGHT, $
  137 + saved_RUN_ANIS, $
  138 + saved_RUN_CIRC, $
  139 + saved_RUN_IONFRAC, $
  140 + saved_RUN_LIN, $
  141 + saved_RUN_POL, $
  142 + saved_RUN_RRF, $
  143 + saved_RUN_TLS, $
  144 + saved_RUN_UNIV, $
  145 + file=file
  146 +
  147 +the_end:
  148 +
  149 +END
... ...
src/idl/dustem_sed_plot.pro
  Diff comments   View file @b5ccb70
1   -PRO dustem_sed_plot,p_min,_EXTRA=extra,function_name=function_name
  1 +PRO dustem_sed_plot,p_min,_EXTRA=extra,function_name=function_name,pol=pol
2 2  
3 3 IF not keyword_set(function_name) THEN BEGIN
4 4 ;Run dustem with as an interface to mpfitfun
5 5 ;stop
6 6 ;== JPB: st used to come out from here
7 7 ;even when not set. Had to replavce by out_st ...
8   - dustem_sed=dustem_compute_sed(p_min,st=st,cont=cont,out_st=out_st)
  8 + dustem_sed=dustem_compute_sed(p_min,st=st,cont=cont,freefree=freefree,synchrotron=synchrotron,out_st=out_st)
9 9 st=out_st
10   - dustem_plot_fit_sed,st,dustem_sed,cont,_extra=extra
  10 + dustem_plot_fit_sed,st,dustem_sed,cont,freefree,synchrotron,_extra=extra
  11 + IF keyword_set(pol) THEN BEGIN
  12 + stop
  13 + ;dustem_polsed=dustem_compute_polsed(p_min,st=pst,cont=cont,freefree=freefree,synchrotron=synchrotron,out_st=out_st)
  14 + dustem_plot_polar,st,ps=ps,_Extra=extra,help=help,UV=UV,SED=SED,Pfrac=Pfrac,print_ratio=print_ratio,win=win,cont=cont,xr=xr,yr=yr,donotclose=donotclose,aligned=aligned,noerrbars=noerrbars,multi=multi,almabands=almabands,nsmooth=nsmooth
  15 + ENDIF
11 16 ENDIF ELSE BEGIN
12 17 CASE function_name OF
13 18 'dustem_greybody_mpfit':BEGIN
... ...
src/idl/dustem_set_func_ind.pro
  Diff comments   View file @b5ccb70
... ... @@ -48,7 +48,7 @@ IDO = ptr_new(iv) ;init_desc ordered (IDO)
48 48 ;==== DEAL WITH THE ISRF<->IONFRAC DEPENDENCY
49 49 fct = strarr(nparam)
50 50 FOR i=0L,nparam-1 DO BEGIN
51   - fct(i) = strmid((*PDO)(i),0,strlen((*PDO)(i))-2.)
  51 + fct[i] = strmid((*PDO)(i),0,strlen((*PDO)(i))-2.)
52 52 ENDFOR
53 53 wisrf = where(fct eq 'dustem_create_isrf' or fct eq 'dustem_create_isrf2',cisrf) ;DUSTEM_CREATE_ISRF is being called as a param ?
54 54 IF cisrf NE 0 THEN BEGIN
... ... @@ -66,6 +66,15 @@ IF cisrf NE 0 THEN BEGIN
66 66 ENDIF
67 67 ENDIF
68 68  
  69 +;=== added by JPB on June 1st 2021 to set !dustem_idl_continuum when needed
  70 +ind=where(fct EQ 'dustem_create_continuum',count)
  71 +IF count NE 0 THEN !dustem_idl_continuum=1
  72 +ind=where(fct EQ 'dustem_create_freefree',count)
  73 +IF count NE 0 THEN !dustem_idl_freefree=1
  74 +ind=where(fct EQ 'dustem_create_synchrotron',count)
  75 +IF count NE 0 THEN !dustem_idl_synchrotron=1
  76 +;stop
  77 +
69 78 ;=== SET !FUNC_IND VALUES
70 79 ; !FUNC_IND = 0 for the DUSTEM_PARAMS
71 80 ; !FUNC_IND = 1 for the first DUSTEM_CREATE function, !FUNC_IND = 2
... ... @@ -79,12 +88,12 @@ FOR i=0L,nparam-1 DO BEGIN
79 88 ;made possible
80 89 ; IF (strmid((*PDO)(i),0,3) NE '(*!') THEN BEGIN
81 90 ;This is the new definition of a function
82   - IF (strmid((*PDO)(i),0,7) EQ 'dustem_') THEN BEGIN
83   - fct = strmid((*PDO)(i),14,strlen((*PDO)(i))-16)
  91 + IF (strmid((*PDO)[i],0,7) EQ 'dustem_') THEN BEGIN
  92 + fct = strmid((*PDO)[i],14,strlen((*PDO)[i])-16)
84 93 IF (fct EQ fct_prev) THEN BEGIN
85   - (*(*!dustem_fit).param_func)(i)=((*!dustem_fit).param_func)(i-1)
  94 + (*(*!dustem_fit).param_func)[i]=(*(*!dustem_fit).param_func)[i-1]
86 95 ENDIF ELSE BEGIN
87   - (*(*!dustem_fit).param_func)(i)=ind
  96 + (*(*!dustem_fit).param_func)[i]=ind
88 97 ind=ind+1
89 98 ENDELSE
90 99 fct_prev=fct
... ...
src/idl/dustem_set_params.pro
  Diff comments   View file @b5ccb70
... ... @@ -23,7 +23,7 @@ FOR i=0L,Nparams-1 DO BEGIN
23 23 length=strlen((*(*!dustem_fit).param_descs)(i))
24 24 strr=strmid((*(*!dustem_fit).param_descs)(i),14,length-16)
25 25 ; stop
26   - IF strr NE 'continuum' THEN BEGIN
  26 + IF strr NE 'continuum' AND strr NE 'freefree' AND strr NE 'synchrotron' THEN BEGIN ;This is not very clean programming ....
27 27 IF strr eq 'isrf2' then strr = 'isrf'
28 28 IF !dustem_which EQ 'VERSTRAETE' THEN BEGIN
29 29 ddir=!dustem_dat+'les_DAT/'
... ...
src/idl/dustem_write_align.pro
  Diff comments   View file @b5ccb70
... ... @@ -24,10 +24,10 @@ printf,unit,st_align.keywords,format=&#39;(A-40)&#39;
24 24 printf,unit,st_align.anisG0,format='(E10.2)'
25 25  
26 26 if !run_univ eq 1 then begin
27   - printf,unit,st_align.grains(0).law,st_align.grains(0).athresh,st_align.grains(0).pstiff,st_align.grains(0).plev,format='(A-10,2X,10(E18.10,2X))'
  27 + printf,unit,st_align.grains(0).law,st_align.grains(0).athresh,st_align.grains(0).pstiff,st_align.grains(0).plev,format='(A-10,2X,10(E18.10,2X))'
28 28 endif else begin
29 29 for i=0,n_elements(st_align.grains)-1 do begin
30   - if st_align.grains(i).aligned then if st_align.grains(i).law ne '' then printf,unit,st_align.grains(i).law,st_align.grains(i).athresh,st_align.grains(i).pstiff,st_align.grains(i).plev,format='(A-10,2X,10(E18.10,2X))'
  30 + if st_align.grains(i).aligned then if st_align.grains(i).law ne '' then printf,unit,st_align.grains(i).law,st_align.grains(i).athresh,st_align.grains(i).pstiff,st_align.grains(i).plev,format='(A-10,2X,10(E18.10,2X))'
31 31 endfor
32 32 endelse
33 33  
... ...
src/idl/dustem_write_all_web3p8.pro
  Diff comments   View file @b5ccb70
... ... @@ -41,8 +41,31 @@ dustem_write_qabs_lv,dir_out_qabs,st.qabs
41 41 dustem_write_calor_lv,dir_out_capa,st.calor
42 42  
43 43 ;== POL
44   -dustem_write_pol,dir_out_dat,st.pol
45   -FOR i_axis = 1, 3 DO dustem_write_qabspol,dir_out_qabs,st,i_axis
  44 +;dustem_write_pol,dir_out_dat,st.pol
  45 +;FOR i_axis = 1, 3 DO dustem_write_qabspol,dir_out_qabs,st,i_axis
  46 +
  47 +;== POL
  48 +if !run_pol then begin
  49 + ;stop
  50 + file_out=dir_out_dat+'ALIGN.DAT'
  51 + ;stop
  52 + ; Pol file
  53 + dustem_write_align,file_out,st.align
  54 +
  55 + ; Linear polarization
  56 + file_out=dir_out_qabs
  57 + IF not stregex(st.align.keywords, 'rrf', /bool) THEN for i_axis = 1, 2 do dustem_write_qpol,file_out,st,i_axis
  58 +
  59 + ; Circular polarization
  60 + if !run_circ then dustem_write_qcirc,file_out,st
  61 +
  62 + ; Anisotropic heating
  63 + IF (st.align.anisG0 > 0) THEN for i_axis = 1, 2 do dustem_write_qh,file_out,st,i_axis
  64 +
  65 + ; RRF Files
  66 + IF stregex(st.align.keywords, 'rrf', /bool) THEN for i_axis = 1, 2 do dustem_write_qpol_rrf,file_out,st,i_axis
  67 +
  68 +endif
46 69  
47 70 ;== TLS
48 71 file_out=dir_out_dat
... ...
src/idl/dustem_write_grain_web3p8.pro
  Diff comments   View file @b5ccb70
... ... @@ -37,11 +37,25 @@ printf,unit,(*!dustem_params).G0
37 37 ;printf,unit,(*!dustem_params).ngrains
38 38 ;===== JPB: updating for longer keywords
39 39 ;frmt='(A12,I4,A12,11E14.6)'
40   -frmt='(A12,I4,A20,11E14.6)'
  40 +;frmt='(A12,I4,A20,11E14.6)'
  41 +frmt='(A20,I4,A20,11E14.6)'
41 42  
42 43 FOR i=0L,(*!dustem_params).ngrains-1 DO BEGIN
43   - printf,unit,st(i).grain_type,st(i).nsize,st(i).type_keywords,st(i).mdust_o_mh,st(i).rho,st(i).amin,st(i).amax, $
44   - st(i).alpha_o_a0, st(i).at,st(i).ac,st(i).gamma,st(i).au,st(i).zeta,st(i).eta,format=frmt
  44 + printf,unit,st(i).grain_type, $
  45 + st(i).nsize, $
  46 + st(i).type_keywords, $
  47 + st(i).mdust_o_mh, $
  48 + st(i).rho, $
  49 + st(i).amin, $
  50 + st(i).amax, $
  51 + st(i).alpha_o_a0, $
  52 + st(i).at, $
  53 + st(i).ac, $
  54 + st(i).gamma, $
  55 + st(i).au, $
  56 + st(i).zeta, $
  57 + st(i).eta, $
  58 + format=frmt
45 59 IF !dustem_verbose EQ 1 THEN BEGIN
46 60 ; stop
47 61 print,st(i).grain_type,st(i).nsize,st(i).type_keywords,st(i).mdust_o_mh,st(i).rho,st(i).amin,st(i).amax, $
... ... @@ -52,5 +66,7 @@ ENDFOR
52 66 close,unit
53 67 free_lun,unit
54 68  
  69 +;stop
  70 +
55 71 END
56 72  
... ...
src/idl/dustem_write_qpol.pro
  Diff comments   View file @b5ccb70
... ... @@ -3,15 +3,21 @@ PRO dustem_write_qpol,dir,st,i_axis
3 3 ;Caution: dir is the directory only
4 4 ;Filenames are constructed inside the routine
5 5  
  6 +;stop
  7 +
6 8 Estring='.DAT'
7   -Nfiles=st.grain.ngrains
  9 +;Nfiles=st.grain.ngrains
  10 +Nfiles=st.ngrains
8 11 frmt="(100E15.5)"
9 12  
10 13 FOR i=0L,Nfiles-1 DO BEGIN
11 14  
  15 +;stop
  16 +
12 17 IF st.align.grains(i).aligned THEN BEGIN
13 18  
14   - filename=dir+'Q'+strtrim(i_axis,2)+'_'+st.grain.grains(i).grain_type+Estring
  19 + ;filename=dir+'Q'+strtrim(i_axis,2)+'_'+st.grain.grains(i).grain_type+Estring
  20 + filename=dir+'Q'+strtrim(i_axis,2)+'_'+st.grains(i).grain_type+Estring
15 21 Nsizes=n_elements((*st.qabs(i)).sizes)
16 22 OPENW,unit,filename,/get_lun
17 23 printf,unit,Nsizes
... ... @@ -20,20 +26,20 @@ FOR i=0L,Nfiles-1 DO BEGIN
20 26 NQabs=n_tags((*st.qpol(i_axis-1,i)).qabs)
21 27 vars=fltarr(NQabs,Nlines)
22 28  
23   - ; Separator
24   - printf,unit,"# QABS"
  29 + ;Separator
  30 + printf,unit,"# QABS"
25 31 FOR k=0L,NQabs-1 DO vars(k,*)=((*st.qpol(i_axis-1,i)).qabs).(k)
26 32 FOR j=0L,Nlines-1 DO printf,unit,vars(*,j),format=frmt
27 33  
28   - ; Separator
29   - printf,unit,"# QSCA"
  34 + ; Separator
  35 + printf,unit,"# QSCA"
30 36 FOR k=0L,NQabs-1 DO vars(k,*)=((*st.qpol(i_axis-1,i)).qsca).(k)
31 37 FOR j=0L,Nlines-1 DO printf,unit,vars(*,j),format=frmt
32 38  
33 39 close,unit
34 40 free_lun,unit
35 41  
36   - endif
  42 + ENDIF
37 43  
38 44 ENDFOR
39 45  
... ...
src/idl/gaunt_ff.pro 0 → 100644
  Diff comments   View file @b5ccb70
... ... @@ -0,0 +1,11 @@
  1 +FUNCTION gaunt_ff,nu,T
  2 +
  3 +;nu in Hz
  4 +nu40=nu/1.e9/40 ;in units of 40 GHz
  5 +T4=T/1.e4
  6 +
  7 +mean_gff=4.4*(T4)^0.21*(nu40)^(-0.14)
  8 +
  9 +RETURN,mean_gff
  10 +
  11 +END
0 12 \ No newline at end of file
... ...
src/idl/intensity_free_free.pro 0 → 100644
  Diff comments   View file @b5ccb70
... ... @@ -0,0 +1,57 @@
  1 +FUNCTION intensity_free_free,nu,T,em,I_halpha_R=I_halpha_R,nHe=nHe,ergs=ergs,mjy=mjy
  2 +
  3 +;nu=frequency in Hz
  4 +;em is emission measure in cm^-6*pc
  5 +;T is the gas temperature
  6 +;example:
  7 +;nu=range_gen(1000,[1.,1000.])*1.e9 & T=5.e4 & em=1.
  8 +;I_ff=intensity_free_free(nu,T,em,I_halpha=I_halpha_R)
  9 +;cgplot,nu/1.e9,I_ff/I_halpha_R,/xlog,/ylog,yr=[1.e-2,1.e4],/ysty,xr=[1,1e3],/xsty
  10 +;print,I_halpha_R
  11 +; 0.0750229
  12 +
  13 +nu_GHz=nu/1.e9
  14 +nu_10=nu_GHz/10. ;in units of 10 GHz
  15 +c=2.998e8 ;m/s
  16 +k=1.38e-23 ;J/K
  17 +
  18 +;compute I Halpha for caseB recombination in Rayleigh
  19 +I_halpha=intensity_h_alpha(em,T) ;this is in erg/cm2/s/sr
  20 +R=5.661e-18 ;1R in ergs/s/cm^2/arcsec^2
  21 +sr=4.254517e10 ;1sr in arcsec^2
  22 +RR=R*sr
  23 +I_halpha_R=I_halpha/RR ;in Rayleigh
  24 +
  25 +Z=1.
  26 +
  27 +;stop
  28 +
  29 +T4=T/1.e4
  30 +;ionization potentials of He
  31 +chi_he=24.5874 ;eV
  32 +chi_Heii=54.41778 ;eV
  33 +use_nHe=0.
  34 +IF keyword_set(nhe) THEN use_nHe=nHe
  35 +
  36 +
  37 +;not sure what the double 0.08 really means ....
  38 +T_muK=I_halpha_R*14.*T4^0.317*10^(0.029/T4)*(1.+0.08*chi_Heii/chi_he*use_nHe/0.08)*Z^2*gaunt_ff(nu,T)/nu_10^2
  39 +;T_muK=I_halpha_R*14.*T4^0.317*10.^(0.029/T4)*Z^2*gaunt_ff(nu,T)/nu_10^2
  40 +;T_muK is the free-free brightness temperature (I=2*h*nu^2*k*Tb/c^2) in muK
  41 +Ifreefree=T_muK
  42 +
  43 +IF keyword_set(ergs) THEN BEGIN
  44 + Ifreefree=2.*nu^2*k*T_muK*1.e6/c^2 ;This is ergs/s/m^2/Hz/sr
  45 +ENDIF
  46 +
  47 +;convert to MJy/sr if needed
  48 +IF keyword_set(mjy) THEN BEGIN
  49 + lambda=c/nu*1.e6 ;mic
  50 + T_mK=T_muK/1.e3
  51 + convert_mk_mjy, c/nu, T_mK, I_Mjy, /RJ
  52 + Ifreefree=I_Mjy
  53 +ENDIF
  54 +
  55 +RETURN,Ifreefree
  56 +
  57 +END
0 58 \ No newline at end of file
... ...
src/idl/intensity_h_alpha.pro 0 → 100644
  Diff comments   View file @b5ccb70
... ... @@ -0,0 +1,21 @@
  1 +FUNCTION intensity_h_alpha,em,T,case_ab=case_ab
  2 +
  3 +;em is emission measure in cm^-6*pc
  4 +;T is gas the temperature in K
  5 +;intensity is in erg/cm2/s/sr
  6 +
  7 +use_case='B'
  8 +IF keyword_set(case_ab) THEN use_case=case_ab
  9 +
  10 +T4=T/1.e4
  11 +
  12 +CASE use_case OF
  13 + ;This is caseB eq 9 of Walls-Gabaud 1998
  14 + 'B':intensity=9.41*1.e-8*T4^(-1.017)*10.^(-0.029/T4)*em
  15 + ;This is caseA eq 10 of Walls-Gabaud 198
  16 + 'A':intensity=6.36*1.e-8*T4^(-1.134)*10.^(-0.038/T4)*em
  17 +ENDCASE
  18 +
  19 +RETURN,intensity
  20 +
  21 +END
0 22 \ No newline at end of file
... ...
src/idl/test_intensity_free_free.pro 0 → 100644
  Diff comments   View file @b5ccb70
... ... @@ -0,0 +1,31 @@
  1 +PRO test_intensity_free_free
  2 +
  3 +;This should reproduce the plot of Fig 2 in Walls-Gabaud PASA 1998, 15 111
  4 +nu=range_gen(1000,[1.,1000.])*1.e9 & T=5.e3 & em=1.
  5 +I_ff=intensity_free_free(nu,T,em,I_halpha=I_halpha_R)
  6 +print,I_halpha_R
  7 +; 0.0750229
  8 +
  9 +cgplot,alog10(nu/1.e9),alog10(I_ff/I_halpha_R),yr=[-2.,4],/ysty,xr=[0,3],/xsty,ytit='log(Tb [muK]/IHalpha [R])',xtit='log(nu [GHz])'
  10 +
  11 +T=1.e4
  12 +I_ff=intensity_free_free(nu,T,em,I_halpha=I_halpha_R)
  13 +cgoplot,alog10(nu/1.e9),alog10(I_ff/I_halpha_R),linestyle=2
  14 +
  15 +T=2.e4
  16 +I_ff=intensity_free_free(nu,T,em,I_halpha=I_halpha_R)
  17 +cgoplot,alog10(nu/1.e9),alog10(I_ff/I_halpha_R),linestyle=3
  18 +
  19 +;==== test converion to MJy/sr or ergs
  20 +nu=range_gen(1000,[1.,1000.])*1.e9 & T=5.e3 & em=1.
  21 +I_ff_erg=intensity_free_free(nu,T,em,I_halpha=I_halpha_R,/erg)
  22 +print,I_halpha_R
  23 +
  24 +cgplot,alog10(nu/1.e9),alog10(I_ff_erg/I_halpha_R),/ysty,xr=[0,3],/xsty,ytit='log(Tb [erg/s/m2/sr/Hz]/IHalpha [R])',xtit='log(nu [GHz])'
  25 +I_ff_mmjjyy=I_ff_erg*1.e20 ;in MJy/sr
  26 +
  27 +I_ff_mjy=intensity_free_free(nu,T,em,I_halpha=I_halpha_R,/mjy)
  28 +cgplot,alog10(nu/1.e9),alog10(I_ff_mjy/I_halpha_R),/ysty,xr=[0,3],/xsty,ytit='log(Tb [MJy/sr]/IHalpha [R])',xtit='log(nu [GHz])'
  29 +cgoplot,alog10(nu/1.e9),alog10(I_ff_mmjjyy/I_halpha_R),color='red'
  30 +
  31 +END
0 32 \ No newline at end of file
... ...