Commit 452c334ec7988ef353fb6f1d7dc14a61c8bd1927

Authored by Ilyes Choubani
1 parent 0e608856
Exists in master

Implementation Of Plugin Formalism

Showing 1341 changed files with 1000 additions and 368 deletions   Show diff stats

Too many changes.

To preserve performance only 100 of 1341 files displayed.

src/gdl_startup_sample 100644 → 100755
src/idl/HFI_UC_CC/get_hfibolo_list.pro 100644 → 100755
src/idl/HFI_UC_CC/hfi_color_correction_dustem.pro 100644 → 100755
src/idl/HFI_UC_CC/hfi_read_bandpass.pro 100644 → 100755
src/idl/Other/read_draine.pro 100644 → 100755
src/idl/Other/read_draine_qabs.pro 100644 → 100755
src/idl/all_dustem_dependencies.xcat 100644 → 100755
src/idl/bg.pro 100644 → 100755
src/idl/bg_mpfit.pro 100644 → 100755
src/idl/dustem_cc.pro 100644 → 100755
src/idl/dustem_cc_mips_vs_handbook.pro 100644 → 100755
src/idl/dustem_collect_wav.pro 100644 → 100755
src/idl/dustem_compute_gb_sed.pro 100644 → 100755
src/idl/dustem_compute_gb_sed_fast.pro 100644 → 100755
src/idl/dustem_compute_polsed.pro 100644 → 100755
1   -FUNCTION dustem_compute_polsed,p_dim,stp,_extra=extra,out_st=out_st
  1 +FUNCTION dustem_compute_polsed,p_dim,stp,_extra=extra,out_st=out_st,dustem_qsed,dustem_used,Q_sed,U_sed
2 2  
3 3 ;+
4 4 ; NAME:
... ... @@ -42,20 +42,75 @@ IF keyword_set(help) THEN BEGIN
42 42 ENDIF
43 43  
44 44 IF not keyword_set(stp) THEN BEGIN
45   - ;st=dustem_run(p_dim)
46   - dustem_create_functions,p_dim/(*(*!dustem_fit).param_init_values),res=res
47 45 stp=dustem_run(p_dim)
  46 + dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),res=res,stp
48 47 ENDIF
49   -If keyword_set(result) THEN result=st
50   -fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/stp.polsed.wav)*1.e20/1.e7
  48 +
  49 +If keyword_set(result) THEN result=stp
  50 +fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/(stp.polsed).wav)*1.e20/1.e7
51 51 spec = stp.polsed.em_tot * fact
52 52 ;stop
  53 +
  54 +;ADDING PLUGIN TO SPECTRUM----------------
  55 +for i=0L,n_tags(*!dustem_scope)-1 do begin
  56 +if total(strsplit((*!dustem_scope).(i),'+',/extract) eq 'ADD_POLSED') then begin
  57 +if tag_exist(*!dustem_scope,'SYNCHROTRON') then begin
  58 +;little segment for the polarization fraction of synchrotron
  59 +indx=where((*(*!dustem_fit).param_descs) EQ 'dustem_create_stokes_3')
  60 +frac_synch=p_dim(indx) ; not initial value for the fake polsed. You need to find a solution
  61 +spec+=(*(*!dustem_plugin).(i))*frac_synch(0)
  62 +endif else begin
  63 +frac_synch=0.7 ;default value
  64 +spec+=(*(*!dustem_plugin).(i))*frac_synch
  65 +endelse
  66 +endif
  67 +endfor
  68 +;-----------------------------------------
  69 +
  70 +
53 71 if not isarray(spec) THEN stop
54 72 ;stop
55 73 ;COMPUTE THE MODEL SED
56 74  
57 75 dustem_polsed = (*!dustem_data.polsed).values * 0.
58 76  
  77 +
  78 +;MODIFYING VIA PLUGIN-------------- THIS BLOCK NEEDS TO BE UNDER THE DUSTEM_POLSED LINE
  79 +for i=0L,n_tags(*!dustem_scope)-1 do begin
  80 +
  81 +if total(strsplit((*!dustem_scope).(i),'+',/extract) eq 'MODIFY') then begin
  82 +
  83 +;making sure it's the stokes plugin that gets used here
  84 +
  85 +iidx=where(tag_names(*!dustem_scope) eq 'SYNCHROTRON',countiidx)
  86 +
  87 +IF countiidx ne 0 THEN BEGIN
  88 +
  89 +Q_dust=(*(*!dustem_plugin).(i))(0:n_elements((*(*!dustem_plugin).(i)))/4-1)
  90 +U_dust=(*(*!dustem_plugin).(i))(n_elements((*(*!dustem_plugin).(i)))/4:n_elements((*(*!dustem_plugin).(i)))/2-1)
  91 +Q_synch=(*(*!dustem_plugin).(i))(n_elements((*(*!dustem_plugin).(i)))/2:3*n_elements((*(*!dustem_plugin).(i)))/4-1)
  92 +U_synch=(*(*!dustem_plugin).(i))(3*n_elements((*(*!dustem_plugin).(i)))/4:*)
  93 +
  94 +ENDIF ELSE BEGIN
  95 +Q_dust=(*(*!dustem_plugin).(i))(0:n_elements((*(*!dustem_plugin).(i)))/2-1)
  96 +U_dust=(*(*!dustem_plugin).(i))(n_elements((*(*!dustem_plugin).(i)))/2:*)
  97 +Q_synch=Q_dust*0
  98 +U_synch=U_dust*0
  99 +ENDELSE
  100 +
  101 +;THE ABOVE VARIABLES SHOULD BE ACCESSED NATURALLY IN THE MPFIT_RUN.PRO
  102 +Q_sed=Q_dust*fact+Q_synch
  103 +U_sed=U_dust*fact+U_synch
  104 +
  105 +;data the dustem_wrap will use
  106 +dustem_qsed = dustem_polsed
  107 +dustem_used = dustem_polsed
  108 +
  109 +endif
  110 +endfor
  111 +;-------------------------------------------
  112 +
  113 +
59 114 ;stop
60 115 IF !dustem_do_cc NE 0 AND !dustem_never_do_cc EQ 0 THEN BEGIN
61 116 message,'DOING color correction calculations',/info
... ... @@ -64,13 +119,19 @@ ENDIF ELSE BEGIN
64 119 ENDELSE
65 120  
66 121 ind_polsed=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_polsed)
67   -
68   -IF count_polsed NE 0 THEN BEGIN
  122 +
  123 +IF count_polsed NE 0 THEN BEGIN ; INCLUDE COLOR CORRECTIONS OF (Q,U) for dust & synch
69 124 filter_names=((*!dustem_data.polsed).filt_names)(ind_polsed)
70   - ;FOR ii=0L,n_elements(filter_names)-1 DO BEGIN
71   - spolsed=dustem_cc(stp.sed.wav,spec,filter_names,cc=cc)
72   - ;ENDFOR
73   - dustem_polsed(ind_polsed)=spolsed
  125 +
  126 + spolsed=dustem_cc(stp.sed.wav,spec,filter_names,cc=cc) & dustem_polsed(ind_polsed)=spolsed
  127 +
  128 +
  129 + ;PLUGIN 'STOKES' W/ scope='MODIFY' PART------------------------------------
  130 +
  131 + if isa(dustem_qsed) then sqsed=dustem_cc(stp.sed.wav,Q_sed,filter_names,cc=cc) & dustem_qsed(ind_polsed)=sqsed
  132 + if isa(dustem_used) then sused=dustem_cc(stp.sed.wav,U_sed,filter_names,cc=cc) & dustem_used(ind_polsed)=sused
  133 + ;--------------------------------------------------------------------------
  134 +
74 135 ENDIF
75 136  
76 137 ;For spectrum data points, interpolate in log-log
... ... @@ -79,9 +140,15 @@ ENDIF
79 140 ind_spec=where((*!dustem_data.polsed).filt_names EQ 'SPECTRUM',count_spec)
80 141 IF count_spec NE 0 THEN BEGIN
81 142 dustem_polsed(ind_spec)=interpol(spec,stp.wav,(((*!dustem_data.polsed).wav)(ind_spec)))
  143 +
  144 + ;STOKES PARAMS TREATMENT
  145 + if isa(dustem_qsed) then dustem_qsed(ind_spec)=interpol(Q_sed,stp.wav,(((*!dustem_data.polsed).wav)(ind_spec)))
  146 + if isa(dustem_used) then dustem_used(ind_spec)=interpol(U_sed,stp.wav,(((*!dustem_data.polsed).wav)(ind_spec)))
  147 +
82 148 ENDIF
83 149 out_st=stp
84 150  
  151 +
85 152 ;clean pointers
86 153 heap_gc
87 154  
... ...
src/idl/dustem_compute_sed.pro 100644 → 100755
1   -FUNCTION dustem_compute_sed,p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron,_extra=extra,out_st=out_st,wave0=wave0
  1 +FUNCTION dustem_compute_sed,p_dim,st,_extra=extra,out_st=out_st,wave0=wave0
2 2  
3 3 ;+
4 4 ; NAME:
... ... @@ -42,16 +42,12 @@ IF keyword_set(help) THEN BEGIN
42 42 goto,the_end
43 43 ENDIF
44 44  
45   -;stop
46 45  
47 46 IF not keyword_set(st) THEN BEGIN
48   - ;st=dustem_run(p_dim)
49   - IF !dustem_idl_continuum THEN cont=dustem_create_continuum()
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
  47 +
53 48 st=dustem_run(p_dim)
54   - ;stop
  49 + dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),res=res,st
  50 +
55 51 ENDIF
56 52  
57 53 ;If keyword_set(result) THEN result=st
... ... @@ -60,10 +56,17 @@ if not keyword_set(wave0) then wave0 = 0.
60 56 ; Convert into MJy/sr/1d20
61 57 fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/st.sed.wav)*1.e20/1.e7
62 58 spec = st.sed.em_tot * fact
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
  59 +
  60 +;ADDING PLUGIN TO SPECTRUM----------------
  61 +;if n_tags(!dustem_data.sed) gt 1 then begin ; executing only when !dustem_data.sed is set meaning that this allows for the creation of fake seds
  62 +for i=0L,n_tags(*!dustem_scope)-1 do begin
  63 +if total(strsplit((*!dustem_scope).(i),'+',/extract) eq 'ADD_SED') then spec+=(*(*!dustem_plugin).(i))
  64 +endfor
  65 +;endif
  66 +;------------------------------------------
  67 +
66 68 if not isarray(spec) THEN stop
  69 +
67 70 ;COMPUTE THE MODEL SED
68 71 dustem_sed = (*!dustem_data.sed).values * 0.
69 72  
... ... @@ -81,7 +84,7 @@ IF count_sed NE 0 THEN BEGIN
81 84 ;ENDFOR
82 85 dustem_sed[ind_sed]=ssed
83 86 ENDIF ELSE BEGIN
84   - print,"**** WARNING **** NO FILTER defined in DATA for COLOUR CORRECTION"
  87 + print,"**** WARNING **** NO FILTER defined in DATA for COLOUR CORRECTION"
85 88 ENDELSE
86 89  
87 90 ;For spectrum data points, interpolate in log-log
... ...
src/idl/dustem_create_continuum.pro 100644 → 100755
1   -FUNCTION dustem_create_continuum,key=key,val=val,help=help
  1 +FUNCTION dustem_create_continuum,key=key,val=val,scope,help=help
2 2  
3 3 ;+
4 4 ; NAME:
... ... @@ -34,6 +34,7 @@ FUNCTION dustem_create_continuum,key=key,val=val,help=help
34 34 ; Written by JPB
35 35 ;-
36 36  
  37 +
37 38 IF keyword_set(help) THEN BEGIN
38 39 doc_library,'dustem_create_continuum'
39 40 output=0.
... ... @@ -46,11 +47,16 @@ ENDIF
46 47 temp=1000.
47 48 ampl=1.d-2
48 49  
  50 +
  51 +
49 52 IF keyword_set(key) THEN BEGIN
50 53 a=where(key EQ 1,count1)
51 54 b=where(key EQ 2,count2)
  55 +
52 56 IF count1 NE 0 then temp=(val(a))(0)
53 57 IF count2 NE 0 then ampl=(val(b))(0)
  58 +
  59 +
54 60 ENDIF
55 61  
56 62 IF !dustem_which EQ 'DESERT' THEN BEGIN
... ... @@ -63,6 +69,19 @@ ENDELSE
63 69 norm = max(dustem_planck_function(temp,lambir))
64 70 output = ampl*dustem_planck_function(temp,lambir)/norm
65 71  
  72 +
  73 +; if keyword_set(scopes) then begin
  74 +; k=execute(scopes)
  75 +; k=execute('((*!dustem_plugins).nir_continuum)=l') ; maybe some changes will have to be done at the dustem_init level
  76 +; print, 'l is:'
  77 +; print, l
  78 +; help, l
  79 +; endif
  80 +; if not keyword_set (scopes) && ((*!dustem_plugins).nir_continuum) eq '' then ((*!dustem_plugins).nir_continuum)=['ADD_SED'] ; DEFAULT VALUE FOR THIS PLUGIN
  81 +
  82 +scope=((*!dustem_scope).continuum) ; FUNCTION RETURNS SCOPE
  83 +;if keyword_set(scope) then
  84 +
66 85 the_end:
67 86 RETURN,output
68 87  
... ...
src/idl/dustem_create_freefree.pro 100644 → 100755
1   -FUNCTION dustem_create_freefree,key=key,val=val,help=help
  1 +FUNCTION dustem_create_freefree,key=key,val=val,scope,help=help
2 2  
3 3 ;+
4 4 ; NAME:
... ... @@ -101,6 +101,10 @@ CASE use_method OF
101 101 END
102 102 ENDCASE
103 103  
  104 +
  105 +scope=((*!dustem_scope).free_free)
  106 +;if keyword_set(scope) then
  107 +
104 108 the_end:
105 109 RETURN,output
106 110  
... ...
src/idl/dustem_create_functions.pro deleted
... ... @@ -1,49 +0,0 @@
1   -PRO dustem_create_functions,p_min,cont=cont,freefree=freefree,synchrotron=synchrotron,res=res
2   -
3   -p_dim = p_min * (*(*!dustem_fit).param_init_values)
4   -
5   -f=1
6   -FOR i=0L,n_elements(p_min)-1 DO BEGIN
7   - ;stop
8   - parameter_type=dustem_parameter_description2type((*(*!dustem_fit).param_descs)[i],string_name=string_name)
9   - ;IF ((*(*!dustem_fit).param_func)(i) NE 0.) THEN BEGIN ;test if this is a parameter for a pluggin
10   - IF parameter_type EQ 'PLUGIN' THEN BEGIN ;test if this is a parameter for a plugin
11   - ftn = strmid((*(*!dustem_fit).param_descs)(i), 0, strlen((*(*!dustem_fit).param_descs)(i))-2) ;caution, this would not allow more than 9 parameters
12   - tmp = (where(*(*!dustem_fit).param_func eq f, count))
13   - PDO_tmp = (*(*!dustem_fit).param_descs)(tmp)
14   - p_dim_tmp = p_dim(tmp)
15   - index = fltarr(count) & value = fltarr(count)
16   - FOR k=0, count-1 DO BEGIN
17   - index[k] = strmid(PDO_tmp[k],0,1,/REV )
18   - value[k] = p_dim_tmp[k]
19   - ENDFOR
20   - IF ftn EQ 'dustem_create_ionfrac' THEN BEGIN
21   - key_ionfrac=index
22   - val_ionfrac=value
23   - ENDIF
24   - IF ftn NE 'dustem_create_continuum' and ftn NE 'dustem_create_freefree' and ftn NE 'dustem_create_synchrotron' THEN BEGIN
25   - str='res='+ftn+'(key=index,val=value)'
26   - toto=execute(str)
27   - IF !dustem_verbose NE 0 THEN message,str,/info
28   - ENDIF ELSE BEGIN
29   - IF ftn EQ 'dustem_create_continuum' THEN BEGIN
30   - str='cont='+ftn+'(key=index,val=value)'
31   - toto=execute(str)
32   - IF !dustem_verbose NE 0 THEN message,str,/info
33   - ENDIF
34   - IF ftn EQ 'dustem_create_freefree' THEN BEGIN
35   - str='freefree='+ftn+'(key=index,val=value)'
36   - toto=execute(str)
37   - IF !dustem_verbose NE 0 THEN message,str,/info
38   - ENDIF
39   - IF ftn EQ 'dustem_create_synchrotron' THEN BEGIN
40   - str='synchrotron='+ftn+'(key=index,val=value)'
41   - toto=execute(str)
42   - IF !dustem_verbose NE 0 THEN message,str,/info
43   - ENDIF
44   - ENDELSE
45   - f=f+1 & i=i+count-1
46   - ENDIF
47   -ENDFOR
48   -
49   -END
src/idl/dustem_create_ionfrac.pro 100644 → 100755
src/idl/dustem_create_synchrotron.pro 100644 → 100755
1   -FUNCTION dustem_create_synchrotron,key=key,val=val,help=help
  1 +FUNCTION dustem_create_synchrotron,key=key,val=val,scope,help=help
2 2  
3 3 ;+
4 4 ; NAME:
5 5 ; dustem_create_synchrotron
6 6 ; PURPOSE:
7   -; DUSTEM pluggin to compute synchrotron emission
  7 +; DUSTEM plugin to compute synchrotron emission
8 8 ; CATEGORY:
9 9 ; DUSTEM Wrapper
10 10 ; CALLING SEQUENCE:
... ... @@ -42,31 +42,38 @@ ENDIF
42 42  
43 43 ;default values of input parameters
44 44 s=3 ;power spectrum of CR E distribution
45   -A=1. ;Synchrotron radiation amplitude at 10 mm
  45 +A=1 ;Synchrotron radiation amplitude at 10 mm
  46 +
46 47  
47 48 IF keyword_set(key) THEN BEGIN
48 49 a=where(key EQ 1,count1)
49 50 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 51  
  52 + IF count1 NE 0 then s=(val(a))(0) ;else s=3
  53 + IF count2 NE 0 then A=(val(b))(0) ;else A=1
  54 +ENDIF
  55 +
54 56 IF !dustem_which EQ 'DESERT' THEN BEGIN
55 57 lambir=((*!dustem_params).gemissiv.lambir)
56 58 ENDIF ELSE BEGIN
57 59 lambir=((*!dustem_params).lambda.lambda)
58 60 ENDELSE
59 61  
  62 +
60 63 cmic=3.e14
61 64 nu=cmic/lambir
62 65 lambir_ref=10000.
63 66  
  67 +
  68 +
64 69 ;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
  70 +output=nu^(-(s+3.)/2.)
  71 +norm=interpol(output,lambir,lambir_ref)
  72 +output=A*output/norm
  73 +
68 74  
69   -output=synch
  75 +scope=((*!dustem_scope).synchrotron)
  76 +;if keyword_set(scope) then
70 77  
71 78 the_end:
72 79 RETURN,output
... ...
src/idl/dustem_data_tag.pro 100644 → 100755
src/idl/dustem_define_isrf_stellar_modifyer_variable.pro 100644 → 100755
src/idl/dustem_extsed_plot.pro 100644 → 100755
src/idl/dustem_filter2fluxconv.pro 100644 → 100755
src/idl/dustem_filter2freq.pro 100644 → 100755
src/idl/dustem_filter2instru.pro 100644 → 100755
src/idl/dustem_filter2reso.pro 100644 → 100755
src/idl/dustem_filter2sig2.pro 100644 → 100755
src/idl/dustem_filter2sigabs.pro 100644 → 100755
src/idl/dustem_filter2unit.pro 100644 → 100755
src/idl/dustem_filter2wav.pro 100644 → 100755
src/idl/dustem_filter_init.pro 100644 → 100755
... ... @@ -85,7 +85,7 @@ IF not keyword_set(dir) THEN BEGIN
85 85 'DESERT':dir_in=!dustem_soft_dir+'DESERT_POST_ISO_MORE/'
86 86 'COMPIEGNE': dir_in=!dustem_soft_dir+'MC_DAT/'
87 87 'VERSTRAETE': dir_in=!dustem_soft_dir+'d_3.5/'
88   - 'WEB3p8': dir_in=!dustem_soft_dir
  88 + 'WEB3p8': dir_in=!dustem_soft_dir ; I CHANGED THIS and put it back to its original form
89 89 ELSE: dir_in=!dustem_soft_dir+'/Data/les_DAT/'
90 90 ENDCASE
91 91 ENDIF ELSE BEGIN
... ...
src/idl/dustem_fit.pro 100644 → 100755
src/idl/dustem_fit_J13.pro 100644 → 100755
src/idl/dustem_fit_modified_bb_readme.pro 100644 → 100755
src/idl/dustem_fit_redshift_readme.pro 100644 → 100755
src/idl/dustem_fit_sed.pro 100644 → 100755
src/idl/dustem_fit_sed_extsed_readme.pro 100644 → 100755
... ... @@ -92,29 +92,41 @@ st=dustem_set_data(ext=extsed,sed=spec,rchi2_weight=rchi2_weight,f_HI=f_HI)
92 92  
93 93 ;=== Set which parameters you want to fit
94 94 pd = [ $
  95 + 'dustem_create_continuum_2',$
  96 + 'dustem_create_powerlaw_1',$
95 97 '(*!dustem_params).gas.G0', $ ;G0
96 98 '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
97 99 '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
98 100 '(*!dustem_params).grains(2).mdust_o_mh' $ ;amCBEx
99 101 ]
100 102  
101   -p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04]
  103 +p_truth=[0.002,0.5,1.,7.8000E-04,7.8000E-04,1.6500E-04]
102 104  
103 105 ;=== set initial value of parameters you want to fit
104 106 ;iv=p_truth ;exact solution
105   -iv=p_truth+[0.,3.e-5,1.e-5,-2.e-5] ;shifted from solution
  107 +iv=p_truth+[0.01,0.07,0.,3.e-5,1.e-5,-2.e-5] ;shifted from solution
  108 +
  109 +ulimed=[0, 0 ,0 , 0 ,0 ,0 ]
  110 +llimed=[1 ,1 ,1 , 1 ,1 ,1 ]
  111 +llims =[0. ,0., 0. ,0. ,0. ,0.]
  112 +
  113 +
  114 +;
  115 +dustem_init_plugins, pd ;will have to be uncommented for this procedure to work.
  116 +
  117 +
106 118  
107   -ulimed=[0 , 0 ,0 ,0 ]
108   -llimed=[1 , 1 ,1 ,1 ]
109   -llims =[0. ,0. ,0. ,0.]
110 119  
111 120 ;== SET THE FITTED PARAMETERS (now done in dustem_init_parinfo)
112 121 dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
113 122  
  123 +
  124 +
  125 +
114 126 ;=== RUN fit
115   -tol=1.e-10
116   -xtol=1.e-10
117   -Nitermax=2 ;maximum number of iteration. This is the criterium which will stop the fit procedure
  127 +tol=1.e-5
  128 +xtol=1.e-5
  129 +Nitermax=6 ;maximum number of iteration. This is the criterium which will stop the fit procedure
118 130 IF keyword_set(itermax) THEN Nitermax=itermax
119 131  
120 132 xrange=[1.,2.e4]
... ...
src/idl/dustem_fit_sed_polsed_readme.pro 100644 → 100755
... ... @@ -68,29 +68,63 @@ sed.filter=filters
68 68 sed.wave=dustem_filter2wav(filters)
69 69 sed.instru=dustem_filter2instru(filters)
70 70 polsed=sed
71   -st=dustem_set_data(ext=ext,polsed=polsed,polext=polext,sed=sed,polfrac=polfrac,rchi2_weight=rchi2_weight,f_HI=f_HI)
  71 +;polext=sed
  72 +
  73 +;stokes fake parameters
  74 +qsed=sed
  75 +used=sed
  76 +
  77 +; ;stokes extinction fake parameters
  78 +; qext=sed
  79 +; uext=sed
  80 +
  81 +st=dustem_set_data(ext=ext,polsed=polsed,polext=polext,sed=sed,polfrac=polfrac,rchi2_weight=rchi2_weight,f_HI=f_HI,qsed=qsed,used=used,qext=qext,uext=uext)
72 82  
73 83 ;=== Set which parameters you want to fit
74 84 pd = [ $
75 85 '(*!dustem_params).G0', $ ;G0
76 86 '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
77 87 '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
78   - '(*!dustem_params).grains(2).mdust_o_mh' $ ;amCBEx
  88 + '(*!dustem_params).grains(2).mdust_o_mh',$ ;amCBEx
  89 + 'dustem_create_synchrotron_2',$
  90 + 'dustem_create_stokes_1',$ ;Polarization angle psi_ref_dust
  91 + 'dustem_create_stokes_2',$;Polarization angle psi_ref_synch
  92 + 'dustem_create_stokes_3'$;Polarization fraction of synchrotron spectrum
  93 + ;'dustem_create_stext_1'$;Polarization angle in extinction PSI_ext
79 94 ]
80 95  
81   -p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04]
  96 +p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04,5e-3,70,20,0.2]
  97 +;p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04,25]
  98 +
  99 +
  100 +
  101 +
  102 +
  103 +dustem_init_plugins, pd
82 104  
83 105 ;=== set initial value of parameters you want to fit
84 106 ;iv=p_truth ;exact solution
85   -iv=p_truth+[0.,3.e-5,1.e-5,-2.e-5] ;shifted from solution
  107 +iv=p_truth+[0.,3.e-5,1.e-5,-2.e-5,5e-2,5,5,0.2];shifted from solution
  108 +
  109 +;iv=p_truth+[0.,3.e-5,1.e-5,-2.e-5,5];shifted from solution
  110 +
  111 +
  112 +ulimed=[0 ,0 ,0 ,0 ,0 ,0 , 0 ,0] ;Polarization angles exist in the range [0-180] degress
  113 +
  114 +;ulimed=[0 ,0 ,0 ,0 ,0] ;Polarization angles exist in the range [0-180] degress
86 115  
87   -ulimed=[0 , 0 ,0 ,0 ]
88   -llimed=[1 , 1 ,1 ,1 ]
89   -llims =[0. ,0. ,0. ,0.]
  116 +
  117 +llimed=[1 , 1 ,1 ,1 ,1 ,1 ,1 ,1]
  118 +;llimed=[1 , 1 ,1 ,1 ,1]
  119 +
  120 +llims =[0. ,0. ,0. ,0. , 0. ,0. ,0. ,0.]
  121 +;llims =[0. ,0. ,0. ,0. ,0.]
90 122  
91 123 ;== SET THE FITTED PARAMETERS (now done in dustem_init_parinfo)
92 124 dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
93 125  
  126 +;dustem_init_fixed_params,fpd,fiv
  127 +
94 128 ;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-chrg-spin'
95 129 ;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-chrg-spin'
96 130 ;(*!dustem_params).grains[4].TYPE_KEYWORDS='pol-plaw-ed'
... ... @@ -102,19 +136,50 @@ help,!run_pol
102 136 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
103 137 sed.error=(sed.spec*0.1);>0.1
104 138  
105   -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
  139 +polsed.spec=dustem_compute_polsed(p_truth,sstp,_extra=extra,dustem_qsed,dustem_used) ;sst is not set on purpose, for dustem_compute_sed to compute the sed using fortran code
106 140 polsed.error=(polsed.spec*0.1)
107 141  
  142 +
  143 +;creating fake polext
  144 +
  145 +;polext.spec=dustem_compute_polext(p_truth,ssttp,dustem_qext,dustem_uext)
  146 +;polext.error=(polext.spec*0.1)
  147 +;
  148 +
  149 +
  150 +;SET qsed and used data
  151 +qsed.spec=dustem_qsed
  152 +qsed.error=qsed.spec*0.1
  153 +
  154 +
  155 +used.spec=dustem_used
  156 +used.error=used.spec*0.1
  157 +
  158 +
  159 +; ;SET qext and uext data
  160 +; qext.spec=dustem_qext
  161 +; qext.error=qext.spec*0.1
  162 +;
  163 +;
  164 +; uext.spec=dustem_uext
  165 +; uext.error=uext.spec*0.1
  166 +;
  167 +
108 168 ;=== SET THE OBSERVATION STRUCTURE
109 169  
110 170 ind=where(sed.wave GE 0.1)
111 171  
112   -st=dustem_set_data(polsed=polsed[ind],sed=sed[ind])
  172 +st=dustem_set_data(polsed=polsed[ind],sed=sed[ind],qsed=qsed[ind],used=used[ind])
  173 +;st=dustem_set_data(sed=sed[ind],polext=polext[ind],qext=qext[ind],uext=uext[ind])
113 174  
114 175 ;=== RUN fit
115   -tol=1.e-10
116   -xtol=1.e-10
117   -Nitermax=5 ;maximum number of iteration. This is the criterium which will stop the fit procedure
  176 +tol=1.e-12
  177 +xtol=1.e-12
  178 +Nitermax=0 ; just an initialization
  179 +for i=0L,n_tags(!dustem_data)-1 do begin
  180 +if ptr_valid(!dustem_data.(i)) then Nitermax+=1
  181 +endfor
  182 +;Nitermax=4 ;maximum number of iteration. This is the criterion which will stop the fit procedure
118 183  
119 184 xrange=[1.,2.e4]
120 185 yrange=[1e-4,1.e3]
... ... @@ -180,3 +245,4 @@ print,'dustem_mpfit_sed executed in ',t2-t1,' sec'
180 245 the_end:
181 246  
182 247 END
  248 +
... ...
src/idl/dustem_fit_sed_readme.pro 100644 → 100755
... ... @@ -129,6 +129,7 @@ CASE use_model OF
129 129 pd = [ $
130 130 ;'(*!dustem_params).gas.G0', $ ;G0
131 131 'dustem_create_continuum_2', $ ;intensity of NIR continuum
  132 +
132 133 ;'!dustem_isrf_star_add[0].distance', $ ;distance of first stellar contribution to ISRF
133 134 '!dustem_isrf_star_add[0].amplitude', $ ;amplitude of first stellar contribution to ISRF
134 135 '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
... ... @@ -138,17 +139,18 @@ CASE use_model OF
138 139 '(*!dustem_params).grains(4).mdust_o_mh' $
139 140 ]
140 141 ;initial parameter values for parameters to be fitted
141   - iv = [ $
  142 + iv = [ $
142 143 ;1.0, $
143   - 0.002, $ ;intensity of NIR continuum
144   - ;0.003, $ ;distance of first stellar contribution to ISRF
145   - 1.e-3, $ ;amplitude of first stellar contribution to ISRF
146   - 7.8e-4, $ ;mass fraction of PAH0
147   - 7.8e-4, $ ;mass fraction of PAH1
  144 + 0.002, $ ;intensity of NIR continuum
  145 + ;0.003, $ ;distance of first stellar contribution to ISRF
  146 + 1.e-3, $ ;amplitude of first stellar contribution to ISRF
  147 + 7.8e-4, $ ;mass fraction of PAH0
  148 + 7.8e-4, $ ;mass fraction of PAH1
148 149 1.65e-4, $ ;mass fraction of amCBEx
149 150 1.45e-3, $ ;mass fraction of amCBEx
150 151 7.8e-3 $ ;mass fraction of
151 152 ]
  153 +
152 154 Npar=n_elements(pd)
153 155 ulimed=replicate(0,Npar)
154 156 llimed=replicate(1,Npar)
... ... @@ -222,9 +224,11 @@ CASE use_model OF
222 224 ENDCASE
223 225  
224 226 dustem_init,model=use_model,pol=polarization
  227 +
225 228 !dustem_verbose=1
226 229 !dustem_show_plot=1
227 230  
  231 +
228 232 ;=== Read sample SED
229 233 ;=== Composite SED from Compiegne et al 2010, gathered by C. Bot
230 234 dir=!dustem_wrap_soft_dir+'/Data/SEDs/'
... ... @@ -237,6 +241,7 @@ ind=where(spec.instru EQ 'FIRAS',count)
237 241 IF count NE 0 THEN spec(ind).error=0.2*spec(ind).spec
238 242  
239 243  
  244 +
240 245 ;=== SET THE OBSERVATION STRUCTURE
241 246 st=dustem_set_data(sed=spec)
242 247  
... ... @@ -245,6 +250,11 @@ dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo
245 250  
246 251 dustem_init_fixed_params,fpd,fiv
247 252  
  253 +dustem_init_plugins, pd ;initializing the scopes of the plugins to be read later
  254 +
  255 +
  256 +
  257 +
248 258 ;=== RUN fit
249 259 tol=1.e-14
250 260 use_Nitermax=2 ;maximum number of iteration. This is the criterium which will stop the fit procedure
... ...
src/idl/dustem_fit_sed_spinning_readme.pro 100644 → 100755
src/idl/dustem_fit_spin_test_readme.pro 100644 → 100755
src/idl/dustem_gb_plot_fit_sed.pro 100644 → 100755
src/idl/dustem_get_filter_name.pro 100644 → 100755
src/idl/dustem_grains_colors.pro 100644 → 100755
src/idl/dustem_greybody_mpfit.pro 100644 → 100755
src/idl/dustem_init.pro 100644 → 100755
... ... @@ -97,9 +97,14 @@ if keyword_set(pol) then begin
97 97 ext: ptr_new(), $
98 98 polext: ptr_new(), $
99 99 polsed: ptr_new(), $
100   - polfrac: ptr_new() $
  100 + polfrac: ptr_new(), $ ; after this line you can put the new variables to fit
  101 + qsed: ptr_new(), $
  102 + used: ptr_new(), $
  103 + qext: ptr_new(), $
  104 + uext: ptr_new() $
101 105 }
102   - rchi2_weight = {sed: 0.,ext: 0.,polext: 0.,polsed: 0.,polfrac:0.}
  106 +; rchi2_weight = {sed: 0.,ext: 0.,polext: 0.,polsed: 0.,polfrac:0.} ;old line
  107 + rchi2_weight = {sed: 0.,ext: 0.,polext: 0.,polsed: 0.,polfrac:0.,qsed:0.,used:0.,qext:0.,uext:0.}
103 108  
104 109 endif else begin
105 110 defsysv, '!dustem_data', { $ ;Data to fit
... ... @@ -152,10 +157,23 @@ defsysv, '!dustem_params', ptr_new() ;Contains the values of all Desert Model pa
152 157 ; 2:call to dustem_create_ionfrac with the SIZE_xxx.DAT files
153 158 ;ENDIF
154 159  
  160 +
  161 +
  162 +;WE DO NOT NEED THESE VARIABLES ANYMORE BUT HAVEN'T TESTED YET
155 163 defsysv, '!dustem_idl_continuum', 0.
156 164 defsysv, '!dustem_idl_freefree', 0.
157 165 defsysv, '!dustem_idl_synchrotron', 0.
158 166  
  167 +;
  168 +;
  169 +; ;a new plugin will have to be added to this list
  170 +; stt={synchrotron:ptr_new(),free_free:ptr_new(), nir_continuum:ptr_new(), grey_body:ptr_new()}
  171 +; defsysv, '!dustem_plugins', ptr_new(stt)
  172 +
  173 +
  174 +
  175 +
  176 +
159 177 defsysv,'!dustem_filters',ptr_new() ;filter information (filled by dustem_filter_init.pro)
160 178 ;; IDL> help,*!dustem_filters,/str
161 179 ;; ** Structure <1b95608>, 17 tags, length=8232, data length=8118, refs=1:
... ... @@ -237,7 +255,9 @@ IF not keyword_set(dir) THEN BEGIN
237 255 ; 'VERSTRAETE': dir_in=!dustem_wrap_soft_dir+'/Data/d_3.5/'
238 256 'VERSTRAETE': dir_in=!dustem_wrap_soft_dir+'/Data/LV_DAT/'
239 257 'WEB3p8': dir_in=!dustem_soft_dir
240   - ELSE: dir_in=!dustem_wrap_soft_dir+'/Data/les_DAT/'
  258 + ;'WEB4p3': dir_in=!dustem_soft_dir ; I ADDED THIS
  259 + ;'WEB4p3': dir_in='/Users/ilyeschoubani/Soft_Librairies/DUSTEM/Data/les_DAT/' ;dir_in=!dustem_soft_dir+ '/Data/' ; I added this
  260 + ELSE: dir_in=!dustem_wrap_soft_dir+'/Data/les_DAT/';!dustem_wrap_soft_dir+'/Data/les_DAT/';I CHANGED THIS!!!!
241 261 ENDCASE
242 262 ENDIF ELSE BEGIN
243 263 dir_in=dir
... ... @@ -270,7 +290,7 @@ CASE !dustem_which OF
270 290 spawn,'rm -rf '+!dustem_res+'/les_RES'
271 291 ENDIF
272 292 END
273   - 'WEB3p8':BEGIN
  293 + 'WEB3p8':BEGIN ; I CHANGED THIS
274 294 IF file_test(!dustem_dat,/dir) THEN BEGIN
275 295 spawn,'rm -rf '+!dustem_dat+'/data'
276 296 spawn,'rm -rf '+!dustem_dat+'/hcap'
... ... @@ -295,7 +315,7 @@ CASE !dustem_which OF
295 315 spawn,'mkdir '+!dustem_res+'/les_RES'
296 316 ENDIF
297 317 END
298   - 'WEB3p8':BEGIN
  318 + 'WEB3p8':BEGIN ; IN CHANGED THIS AS WELL
299 319 IF file_test(!dustem_dat,/dir) THEN BEGIN
300 320 spawn,'mkdir '+!dustem_dat+'/data'
301 321 spawn,'mkdir '+!dustem_dat+'/hcap'
... ... @@ -326,7 +346,7 @@ ENDIF
326 346 dustem_filter_init,dir=dir
327 347  
328 348 ;=== test plugins
329   -plugins_desc_file=getenv('DUSTEM_SOFT_DIR')+'plugins_description.xcat'
  349 +plugins_desc_file=getenv('DUSTEM_WRAP_SOFT_DIR')+'/plugins_description.xcat' ;added a slash and changed DUSTEM_SOFT_DIR to DUSTEM_WRAP_SOFT_DIR
330 350 plugins_st=read_xcat(plugins_desc_file,/silent)
331 351 Nplugins=n_elements(plugins_st)
332 352 FOR i=0L,Nplugins-1 DO BEGIN
... ...
src/idl/dustem_init_fixed_params.pro 100644 → 100755
... ... @@ -24,7 +24,8 @@ dustem_set_func_ind,fpd,fiv
24 24 ind=where(fpd NE '',count)
25 25 IF count NE 0 THEN BEGIN
26 26 dustem_set_params,[fiv]
27   - dustem_create_functions,[fiv]/(*(*!dustem_fit).param_init_values)
  27 + ;dustem_create_functions,[fiv]/(*(*!dustem_fit).param_init_values)
  28 + dustem_activate_plugins,[fiv]/(*(*!dustem_fit).param_init_values)
28 29 dustem_set_params,[fiv]
29 30 ENDIF
30 31  
... ...
src/idl/dustem_init_parinfo.pro 100644 → 100755
... ... @@ -4,6 +4,8 @@ PRO dustem_init_parinfo,pd,iv, $
4 4 relstep=relstep,step=step
5 5  
6 6 (*!dustem_fit).param_descs=ptr_new(pd)
  7 +;A chance to include a plugin dependency here since pd is not really used in this routine
  8 +
7 9 (*!dustem_fit).param_init_values=ptr_new(iv)
8 10 ;stop
9 11 dustem_set_func_ind,pd,iv
... ...
src/idl/dustem_init_st.pro 100644 → 100755
src/idl/dustem_initialize_sed.pro 100644 → 100755
src/idl/dustem_instru2filters.pro 100644 → 100755
src/idl/dustem_mpfit_data.pro 100644 → 100755
src/idl/dustem_mpfit_run.pro 100644 → 100755
... ... @@ -47,45 +47,52 @@ IF keyword_set(help) THEN BEGIN
47 47 ENDIF
48 48  
49 49 p_dim=p_min*(*(*!dustem_fit).param_init_values)
50   -;Run dustem with as an interface to mpfitfun
51 50  
52   -;==== JPB: not sure why this is required since plugins will be called in dustem_create_functions below ...
53   -;IF !dustem_idl_continuum THEN cont=dustem_create_continuum()
54   -;IF !dustem_idl_freefree THEN freefree=dustem_create_freefree()
55   -;IF !dustem_idl_synchrotron THEN synchrotron=dustem_create_synchrotron()
56   -
57   -;CALL THE DUSTEM_CREATE FUNCTIONS
58   -dustem_create_functions,p_dim/(*(*!dustem_fit).param_init_values),cont=cont,freefree=freefree,synchrotron=synchrotron,res=res
59   -
60   -;if ISRF has been modified and IONFRAC has not been modified, then call dustem_create_ionfrac
61   -;NF on October 2009 : handling of DUSTEM_CREATE_IONFRAC has been changed
62   -;this call only if IONFRACPAH master-keyword in GRAIN.DAT
63   -;Removed for now (JPB). Use of ionfrac should be fully revisited.
64 51 IF !run_ionfrac gt 0. THEN toto = dustem_create_ionfrac()
  52 +
65 53 ;RUN THE MODEL AND GET THE SPECTRUM
66 54 st = dustem_run(p_dim)
67 55 st_model=dustem_read_all(!dustem_soft_dir)
68 56  
  57 +dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),res=res,st
  58 +
  59 +
69 60 chi2_sed = 0
70 61 chi2_ext = 0
71 62 chi2_polext = 0
72 63 chi2_polsed = 0
73   -chi2_polfrac = 0
  64 +chi2_polfrac = 0
  65 +chi2_qsed =0
  66 +chi2_used =0
  67 +chi2_qext =0
  68 +chi2_uext =0
74 69 rchi2_sed = 0
75 70 rchi2_ext = 0
76 71 rchi2_polext = 0
77 72 rchi2_polsed = 0
78 73 rchi2_polfrac = 0
  74 +rchi2_qsed =0
  75 +rchi2_used =0
  76 +rchi2_qext =0
  77 +rchi2_uext =0
79 78 wrchi2_sed = 0
80 79 wrchi2_ext = 0
81 80 wrchi2_polext = 0
82 81 wrchi2_polsed = 0
83 82 wrchi2_polfrac = 0
  83 +wrchi2_qsed =0
  84 +wrich2_used =0
  85 +wrchi2_qext =0
  86 +wrich2_uext =0
84 87 n_sed = 0
85 88 n_ext = 0
86 89 n_polext = 0
87 90 n_polsed = 0
88 91 n_polfrac = 0
  92 +n_qsed =0
  93 +n_used =0
  94 +n_qext =0
  95 +n_uext =0
89 96  
90 97 win = 0
91 98  
... ... @@ -104,8 +111,10 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
104 111 CASE tagnames(ind_data_tag(i_tag)) OF
105 112  
106 113 'SED' : BEGIN
107   -
108   - dustem_sed = dustem_compute_sed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
  114 +
  115 + ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED
  116 +
  117 + dustem_sed = dustem_compute_sed(p_dim,st);dustem_compute_sed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
109 118 chi2_sed = total(((*!dustem_data.sed).values-dustem_sed)^2/(*!dustem_data.sed).sigma^2)
110 119 n_sed = n_elements((*!dustem_data.sed).values)
111 120 rchi2_sed = chi2_sed / n_sed
... ... @@ -126,8 +135,17 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
126 135  
127 136 ; Plot if needed
128 137 IF !dustem_show_plot NE 0 THEN BEGIN
129   - window,0
130   - 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
  138 + ;title='DUSTEM WRAP('+strtrim(string(win),2)+')'
  139 + window,win,title='DUSTEM WRAP (SED)'
  140 +
  141 +
  142 + ;cgwindow,WMULTI=[2,2,1],/CURRENT
  143 +
  144 +
  145 + ;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
  146 + dustem_plot_fit_sed,st,dustem_sed,_extra=extra,res=p_min*(*(*!dustem_fit).param_init_values),chi2=(*!dustem_fit).chi2,rchi2=(*!dustem_fit).rchi2
  147 +
  148 +
131 149 ;stop
132 150 ENDIF
133 151  
... ... @@ -136,8 +154,17 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
136 154 END
137 155  
138 156 'EXT' : BEGIN
139   -
  157 +
  158 + ;ADDING PLUGIN TO SPECTRUM----------------
  159 + ;if n_tags(!dustem_data.ext) gt 1 then begin
  160 + for i=0L,n_tags(*!dustem_scope)-1 do begin
  161 + if total(strsplit((*!dustem_scope).(i),'+',/extract) eq 'ADD_EXT') then st.ext.ext_tot+=(*(*!dustem_plugin).(i))
  162 + endfor
  163 + ;endif
  164 + ;------------------------------------------
  165 +
140 166 dustem_ext = interpol(st.ext.ext_tot,st.ext.wav,(*!dustem_data.ext).wav)
  167 +
141 168 chi2_ext = total(((*!dustem_data.ext).values-dustem_ext)^2/(*!dustem_data.ext).sigma^2)
142 169 n_ext = n_elements((*!dustem_data.ext).values)
143 170 rchi2_ext = chi2_ext / n_ext
... ... @@ -154,7 +181,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
154 181 yr=[1.e-26,1.e-21]
155 182 dustem_plot_ext,st,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,/donotclose,multi=1
156 183 dustem_plot_ext,st,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,mergePAH=mergePAH,multi=2
157   -
  184 +
158 185 win+=1
159 186 ;dustem_plot_ext,st,/UV,yr=[0,2.5],/ysty,xr=[0,10],/xsty,win=win
160 187 xr=[0.2,10]
... ... @@ -179,10 +206,10 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
179 206 endif
180 207  
181 208 IF !run_lin then begin
182   -
183   - dustem_polext = interpol(st.polext.ext_tot,st.polext.wav,(*!dustem_data.polext).wav)
  209 +
  210 + dustem_polext = dustem_compute_polext(p_dim,st,dustem_qext,dustem_uext,Q_ext,U_ext)
184 211 chi2_polext = total(((*!dustem_data.polext).values-dustem_polext)^2/(*!dustem_data.polext).sigma^2)
185   - n_polext = n_elements((*!dustem_data.polext).values)
  212 + n_polext = n_elements((*!dustem_data.polext).values)
186 213 rchi2_polext = chi2_polext / n_polext
187 214 wrchi2_polext = rchi2_polext * !fit_rchi2_weight.polext
188 215 print,"chi2 POLEXT = ",chi2_polext," rchi2 POLEXT = ",rchi2_polext;," weighted rchi2 POLEXT=",wrchi2_polext
... ... @@ -227,21 +254,22 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
227 254 ; Plot if needed
228 255 IF !dustem_show_plot NE 0 THEN BEGIN
229 256  
230   - ; Polarization fraction in the UV-IR
231   - win+=1
232   -
233   - dustem_plot_polar,st,aligned=st_model.align.grains.aligned,yr=[0,0.15],/ysty,xr=[0.1,4],/xsty,win=win
234   -
  257 +; ; Polarization fraction in the UV-IR
  258 +; win+=1
  259 +;
  260 +; dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,yr=[0,0.15],/ysty,xr=[0.1,50],/xsty,win=win;multi=0,ps=filename
  261 +; ;
235 262 ; Polarisation curve (Serkowski)
236   - win+=1
237   - xr=[0.1,50]
238   - yr=[5e-25,5e-23]
239   - 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
240   - 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
  263 + win+=1
  264 + ;xr=[0.3,1e4]
  265 + xr=[1e-4,3.5]
  266 + yr=[1e-30,5e-24]
  267 + dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,/donotclose,multi=1
  268 + dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/UV,yr=yr,/ysty,xr=xr,/xsty,ps=filename,win=win,multi=2
241 269  
242   - win+=1
243   - dustem_plot_fpol,win=win,xr=[10,1000],yr=[0,1.1]
244   -
  270 +; win+=1
  271 +; dustem_plot_fpol,win=win,xr=[10,1000],yr=[0,1.1]
  272 +;
245 273  
246 274 ENDIF
247 275  
... ... @@ -253,50 +281,60 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
253 281  
254 282 END
255 283  
256   - 'POLFRAC': BEGIN
257   -
258   - IF !run_lin then begin
259   -
260   - dustem_polsed = dustem_compute_polsed(p_dim,st.polsed,cont=cont,freefree=freefree,synchrotron=synchrotron)
261   - dustem_sed_tmp = dustem_compute_sed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
262   - dustem_sed = dustem_polsed
263   - ind_filt=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_filt)
264   - for i = 0, count_filt-1 do begin
265   - j=where((*!dustem_data.polsed).filt_names(ind_filt(i)) eq (*!dustem_data.sed).filt_names,count)
266   - if count eq 1 then dustem_sed(i) = dustem_sed_tmp(j(0))
267   - endfor
268   - dustem_polfrac = dustem_polsed/dustem_sed
269   -
270   - ; We normalize at 353GHz
271   - ;x=min(((*!dustem_data.polfrac).wav-850)^2,j353)
272   - ;dustem_polfrac = dustem_polfrac/dustem_polfrac(j353)*(*!dustem_data.polfrac).values(j353)
273   -
274   - chi2_polfrac = total(((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2)
275   - n_polfrac = n_elements((*!dustem_data.polfrac).values)
276   - rchi2_polfrac = chi2_polfrac / n_polfrac
277   - wrchi2_polfrac = rchi2_polfrac * !fit_rchi2_weight.polfrac
278   - print,"chi2 POLFRAC = ",chi2_polfrac," rchi2 POLFRAC = ",rchi2_polfrac;," weighted rchi2 POLFRAC=",wrchi2_polfrac
279   - ;print,"chi2 per wl: ",((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2
280   -
281   - dustem_out = [ dustem_out, dustem_polfrac ]
282   -
283   - ; Plot if needed
284   - IF !dustem_show_plot NE 0 THEN BEGIN
285   -
286   - ; Also show prediction for polarized P/I
287   - win+=1
288   - 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]
289   -
290   - ENDIF
291   -
292   - ENDIF
293   -
294   - END
295   -
  284 +;USE OF POLFRAC DEPRECATED---------------------------------------
  285 +
  286 +; 'POLFRAC': BEGIN
  287 +;
  288 +; ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_POLSED AND DUSTEM_COMPUTE_SED
  289 +;
  290 +; IF !run_lin then begin
  291 +;
  292 +; dustem_polsed = dustem_compute_polsed(p_dim,st);dustem_compute_polsed(p_dim,st.polsed,cont=cont,freefree=freefree,synchrotron=synchrotron)
  293 +; dustem_sed_tmp = dustem_compute_sed(p_dim,st);dustem_compute_sed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
  294 +; dustem_sed = dustem_polsed
  295 +; ind_filt=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_filt)
  296 +; for i = 0, count_filt-1 do begin
  297 +; j=where((*!dustem_data.polsed).filt_names(ind_filt(i)) eq (*!dustem_data.sed).filt_names,count)
  298 +; if count eq 1 then dustem_sed(i) = dustem_sed_tmp(j(0))
  299 +; endfor
  300 +; dustem_polfrac = dustem_polsed/dustem_sed
  301 +;
  302 +; ; We normalize at 353GHz
  303 +; ;x=min(((*!dustem_data.polfrac).wav-850)^2,j353)
  304 +; ;dustem_polfrac = dustem_polfrac/dustem_polfrac(j353)*(*!dustem_data.polfrac).values(j353)
  305 +;
  306 +; chi2_polfrac = total(((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2)
  307 +; n_polfrac = n_elements((*!dustem_data.polfrac).values)
  308 +; rchi2_polfrac = chi2_polfrac / n_polfrac
  309 +; wrchi2_polfrac = rchi2_polfrac * !fit_rchi2_weight.polfrac
  310 +; print,"chi2 POLFRAC = ",chi2_polfrac," rchi2 POLFRAC = ",rchi2_polfrac;," weighted rchi2 POLFRAC=",wrchi2_polfrac
  311 +; ;print,"chi2 per wl: ",((*!dustem_data.polfrac).values-dustem_polfrac)^2/(*!dustem_data.polfrac).sigma^2
  312 +;
  313 +; dustem_out = [ dustem_out, dustem_polfrac ]
  314 +;
  315 +; ; Plot if needed
  316 +; IF !dustem_show_plot NE 0 THEN BEGIN
  317 +;
  318 +; ; Also show prediction for polarized P/I
  319 +; win+=1
  320 +; dustem_plot_polar,st,aligned=st_model.align.grains.aligned,/Pfrac,win=win,yr=[0,0.25],/ysty,xr=[10,5d3];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]
  321 +;
  322 +; ENDIF
  323 +;
  324 +; ENDIF
  325 +;
  326 +; END
  327 +
  328 +
  329 +;-----------------------------------------------------------------
296 330 'POLSED': BEGIN
  331 +
  332 + ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED
  333 +
297 334 IF !run_lin THEN BEGIN
298 335  
299   - dustem_polsed = dustem_compute_polsed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
  336 + dustem_polsed = dustem_compute_polsed(p_dim,st,dustem_qsed,dustem_used,Q_sed,U_sed);dustem_compute_polsed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
  337 +
300 338 chi2_polsed = total(((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2)
301 339 n_polsed = n_elements((*!dustem_data.polsed).values)
302 340 rchi2_polsed = chi2_polsed / n_polsed
... ... @@ -304,20 +342,229 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
304 342 message,"chi2 POLSED = "+strtrim(chi2_polsed,2)+" rchi2 POLSED = "+strtrim(rchi2_polsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
305 343 print,"chi2 per wl: ",((*!dustem_data.polsed).values-dustem_polsed)^2/(*!dustem_data.polsed).sigma^2
306 344 ;print,"SNR per wl: ",(*!dustem_data.polsed).values/(*!dustem_data.polsed).sigma
307   -
308   - dustem_out = [ dustem_out, dustem_polsed ]
  345 +
  346 + dustem_out = [ dustem_out, dustem_polsed ] ; classic command
  347 + ;-----------------------------------------------
  348 +
309 349 ; Plot if needed
310 350 IF !dustem_show_plot NE 0 THEN BEGIN
311 351 win+=1
312   - xr=[10,5e3]
  352 +;
  353 + xr=[10,8e3]
313 354 yr=[1e-34,1e-28]
314   - 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
315   - 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
  355 +
  356 +
  357 + dustem_plot_polar,st,p_dim,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
  358 + dustem_plot_polar,st,p_dim,aligned=st_model.align.grains.aligned,/SED,ps=filename,win=win,cont=cont,freefree=freefree,synchrotron=synchrotron,xr=xr,yr=yr,/xsty,multi=2
  359 +
  360 +
  361 +
316 362 ENDIF
317 363 ENDIF
318 364 END
319   -
320   - ELSE : stop
  365 +
  366 + 'QSED': BEGIN
  367 +
  368 + ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED
  369 + chi2_qsed =total(((*!dustem_data.qsed).values-dustem_qsed)^2/(*!dustem_data.qsed).sigma^2)
  370 + n_qsed = n_elements((*!dustem_data.qsed).values)
  371 + rchi2_qsed = chi2_qsed / n_qsed
  372 + wrchi2_qsed = rchi2_qsed * !fit_rchi2_weight.qsed
  373 + dustem_out = [ dustem_out, dustem_qsed ] ; classic command
  374 + message,"chi2 QSED = "+strtrim(chi2_qsed,2)+" rchi2 QSED = "+strtrim(rchi2_qsed,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
  375 + print,"chi2 per wl: ",((*!dustem_data.qsed).values-dustem_qsed)^2/(*!dustem_data.qsed).sigma^2
  376 +
  377 +
  378 +
  379 + if !dustem_show_plot NE 0 then begin
  380 + fact = (3.e10/(1.e-4*(st.polsed.wav)))/1.d40
  381 + fact1 = (3.e10/(1.e-4*((*!dustem_data.qsed).wav)))/1.d40 ; for data plotting
  382 +
  383 +
  384 + win+=1
  385 + window ,win ,title='DUSTEM WRAP (QSED)'
  386 + titstq='Stokes Q Polarization Intensity'
  387 + ytitstq=textoidl('\nuQ_\nu/N_H (W/m^2/sr/H)')
  388 + xtit=textoidl('\lambda (\mum)')
  389 + xr=[10,8e3]
  390 + normpol=Q_sed*fact
  391 + normpol1=dustem_qsed*fact1
  392 + ;normpol1=(*!dustem_data.qsed).values*fact1
  393 +
  394 + yr=[1e-34,1e-28]
  395 + indq=where(Q_sed lt 0, countq)
  396 + ylog=1
  397 + if countq ne 0 then begin
  398 + ylog=0
  399 + yr=[-5e-31,5e-31]
  400 + endif
  401 +
  402 +
  403 + cgplot,st.polsed.wav,Q_sed*fact,xtit='',ytit=ytitstq,tit=titstq,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)'
  404 + cgoplot,st.polsed.wav,Q_sed*fact,color='black'
  405 + cgoplot,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1,color='Tomato',psym=16,symsize=1,thick=2
  406 + err_bar,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1,yrms=3.*((*!dustem_data.qsed).sigma)/2.*fact1,color=cgColor('Tomato')
  407 + cgoplot,(*!dustem_data.qsed).wav,dustem_qsed*fact1,psym=6,color='red',symsize=2
  408 +
  409 + cgplot,st.polsed.wav,Q_sed*fact/normpol,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1
  410 + cgoplot,st.polsed.wav,Q_sed*fact/normpol,color='black'
  411 + cgoplot,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1/normpol1,color='Tomato',psym=16,symsize=1,thick=2
  412 + err_bar,(*!dustem_data.qsed).wav,(*!dustem_data.qsed).values*fact1/normpol1,yrms=3.*((*!dustem_data.qsed).sigma)/2.*fact1/normpol1,color=cgColor('Tomato')
  413 +
  414 + endif
  415 +
  416 + END
  417 +
  418 + 'USED': BEGIN
  419 +
  420 + ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED
  421 + chi2_used =total(((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2)
  422 + n_used = n_elements((*!dustem_data.used).values)
  423 + rchi2_used = chi2_used / n_used
  424 + wrchi2_used = rchi2_used * !fit_rchi2_weight.used
  425 + dustem_out = [ dustem_out, dustem_used] ; classic command
  426 + message,"chi2 USED = "+strtrim(chi2_used,2)+" rchi2 USED = "+strtrim(rchi2_used,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
  427 + print,"chi2 per wl: ",((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2
  428 +
  429 +
  430 + if !dustem_show_plot ne 0 then begin
  431 + fact = (3.e10/(1.e-4*(st.polsed.wav)))/1.d40
  432 + fact1 = (3.e10/(1.e-4*((*!dustem_data.used).wav)))/1.d40 ; for data plotting
  433 +
  434 +
  435 + win+=1
  436 + window ,win ,title='DUSTEM WRAP (USED)'
  437 + titstu='Stokes U Polarization Intensity'
  438 + ytitstu=textoidl('\nuU_\nu/N_H (W/m^2/sr/H)')
  439 + xtit=textoidl('\lambda (\mum)')
  440 + xr=[10,8e3]
  441 + normpol2=U_sed*fact
  442 + normpol3=dustem_used*fact1
  443 + ;yr=[min(normpol2)-min(normpol2)/10,max(normpol2)+max(normpol2)*10]
  444 +
  445 + yr=[1e-34,1e-28]
  446 + indu=where(U_sed lt 0, countu)
  447 + ylog=1
  448 + if countu ne 0 then begin
  449 + ylog=0
  450 + yr=[-5e-31,5e-31]
  451 + endif
  452 +
  453 + cgplot,st.polsed.wav,U_sed*fact,xtit='',ytit=ytitstu,tit=titstu,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black'
  454 + cgoplot,st.polsed.wav,U_sed*fact,color='black'
  455 + cgoplot,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1,color='Teal',psym=16,symsize=1,thick=2
  456 + err_bar,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1,yrms=3.*((*!dustem_data.used).sigma)/2.*fact1,color=cgColor('Teal')
  457 + cgoplot,(*!dustem_data.used).wav,dustem_used*fact1,psym=6,color='red',symsize=2
  458 +
  459 +
  460 +
  461 + cgplot,st.polsed.wav,U_sed*fact/normpol2,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1
  462 + cgoplot,st.polsed.wav,U_sed*fact/normpol2,color='black'
  463 + cgoplot,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1/normpol3,color='Teal',psym=16,symsize=1,thick=2
  464 + err_bar,(*!dustem_data.used).wav,(*!dustem_data.used).values*fact1/normpol3,yrms=3.*((*!dustem_data.used).sigma)/2.*fact1/normpol3,color=cgColor('Teal')
  465 +
  466 + ;cgwindow, 'plot',
  467 +
  468 + endif
  469 + END
  470 +
  471 +
  472 +
  473 +
  474 + 'QEXT': BEGIN ; Q - EXTINCTION CROSS SECTION
  475 +
  476 + chi2_qext =total(((*!dustem_data.qext).values-dustem_qext)^2/(*!dustem_data.qext).sigma^2)
  477 + n_qext = n_elements((*!dustem_data.qext).values)
  478 + rchi2_qext = chi2_qext / n_qext
  479 + wrchi2_qext = rchi2_qext * !fit_rchi2_weight.qext
  480 + dustem_out = [ dustem_out, dustem_qext] ; classic command
  481 + message,"chi2 QEXT = "+strtrim(chi2_qext,2)+" rchi2 QEXT = "+strtrim(rchi2_qext,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
  482 + print,"chi2 per wl: ",((*!dustem_data.qext).values-dustem_qext)^2/(*!dustem_data.qext).sigma^2
  483 + if !dustem_show_plot ne 0 then begin
  484 +
  485 + win+=1
  486 + window ,win ,title='DUSTEM WRAP (QEXT)'
  487 + titqex='Q-polarized Cross-section'
  488 + ytitqex=textoidl('\sigma_{Q}/N_H (cm^2/H)')
  489 + xtit=textoidl('\lambda^{-1} (\mum^{-1})')
  490 + ;xr=[0.3,1e4]
  491 + xr=[1e-4,3.5]
  492 +
  493 + yr=[1e-34,5e-24]
  494 + indqex=where(Q_ext lt 0, countqex)
  495 + ylog=1
  496 + if countqex ne 0 then begin
  497 + ylog=0
  498 + yr=[-9e-27,5e-27]
  499 + endif
  500 + normpol = st.polext.ext_tot * 0. + 1.d21
  501 + cgplot,1/st.polext.wav,Q_ext,xtit='',ytit=ytitqex,tit=titqex,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black'
  502 + cgoplot,1/st.polext.wav,Q_ext/normpol,color='black'
  503 + cgoplot,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/normpol,color='Indian Red',psym=16,symsize=1,thick=2
  504 + err_bar,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/normpol,yrms=3.*((*!dustem_data.qext).sigma)/2/normpol,color=cgColor('Indian Red')
  505 + cgoplot,1/(*!dustem_data.qext).wav,dustem_qext/normpol,psym=6,color='red',symsize=2
  506 +
  507 +
  508 +
  509 + cgplot,1/st.polsed.wav,Q_ext/Q_ext,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1
  510 + cgoplot,1/st.polsed.wav,Q_ext/Q_ext,color='black'
  511 + cgoplot,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/dustem_qext,color='Indian Red',psym=16,symsize=1,thick=2
  512 + err_bar,1/(*!dustem_data.qext).wav,(*!dustem_data.qext).values/dustem_qext,yrms=3.*((*!dustem_data.qext).sigma)/2.,color=cgColor('Indian Red')
  513 +
  514 + endif
  515 +
  516 + END
  517 +;
  518 + 'UEXT': BEGIN
  519 +
  520 +
  521 + chi2_uext =total(((*!dustem_data.uext).values-dustem_uext)^2/(*!dustem_data.uext).sigma^2)
  522 + n_uext = n_elements((*!dustem_data.uext).values)
  523 + rchi2_uext = chi2_uext / n_uext
  524 + wrchi2_uext = rchi2_uext * !fit_rchi2_weight.uext
  525 + dustem_out = [ dustem_out, dustem_uext] ; classic command
  526 + message,"chi2 UEXT = "+strtrim(chi2_uext,2)+" rchi2 UEXT = "+strtrim(rchi2_uext,2),/continue;," weighted rchi2 POLSED=",wrchi2_polsed
  527 + print,"chi2 per wl: ",((*!dustem_data.uext).values-dustem_uext)^2/(*!dustem_data.uext).sigma^2
  528 + if !dustem_show_plot ne 0 then begin
  529 +
  530 + win+=1
  531 + window ,win ,title='DUSTEM WRAP (UEXT)'
  532 + tituex='U-polarized Cross-section'
  533 + ytituex=textoidl('\sigma_{U}/N_H (cm^2/H)')
  534 + xtit=textoidl('\lambda^{-1} (\mum^{-1})')
  535 + ;xr=[0.3,1e4]
  536 + xr=[1e-4,3.5]
  537 +
  538 + yr=[1e-34,5e-24]
  539 + induex=where(U_ext lt 0, countuex)
  540 + ylog=1
  541 + if countuex ne 0 then begin
  542 + ylog=0
  543 + yr=[-9e-27,5e-27]
  544 + endif
  545 + normpol = st.polext.ext_tot * 0. + 1.d21
  546 + cgplot,1/st.polext.wav,U_ext,xtit='',ytit=ytituex,tit=tituex,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)',color='black'
  547 + cgoplot,1/st.polext.wav,U_ext/normpol,color='black'
  548 + cgoplot,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/normpol,color='Steel Blue',psym=16,symsize=1,thick=2
  549 + err_bar,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/normpol,yrms=3.*((*!dustem_data.uext).sigma)/2/normpol,color=cgColor('Steel Blue')
  550 + cgoplot,1/(*!dustem_data.uext).wav,dustem_uext/normpol,psym=6,color='red',symsize=2
  551 +
  552 +
  553 +
  554 + cgplot,1/st.polsed.wav,U_ext/U_ext,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1
  555 + cgoplot,1/st.polsed.wav,U_ext/U_ext,color='black'
  556 + cgoplot,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/dustem_uext,color='Steel Blue',psym=16,symsize=1,thick=2
  557 + err_bar,1/(*!dustem_data.uext).wav,(*!dustem_data.uext).values/dustem_uext,yrms=3.*((*!dustem_data.uext).sigma)/2/normpol,color=cgColor('Steel Blue')
  558 +
  559 + endif
  560 +
  561 +
  562 +
  563 +
  564 + END
  565 +
  566 +
  567 + ELSE : stop
321 568  
322 569 ENDCASE
323 570  
... ... @@ -344,23 +591,49 @@ if !fit_rchi2_weight.polext ne 0 then begin
344 591 total_rchi2 += rchi2_polext
345 592 endif
346 593 if !fit_rchi2_weight.polsed ne 0 then begin
347   - total_n += n_polsed
348   - total_chi2 += chi2_polsed
349   - total_rchi2 += rchi2_polsed
  594 + total_n += n_polsed
  595 + total_chi2 += chi2_polsed
  596 + total_rchi2 += rchi2_polsed
350 597 endif
351 598 if !fit_rchi2_weight.polfrac ne 0 then begin
352 599 total_n += n_polfrac
353 600 total_chi2 += chi2_polfrac
354 601 total_rchi2 += rchi2_polfrac
355 602 endif
  603 +if !fit_rchi2_weight.qsed ne 0 then begin
  604 + total_n += n_qsed
  605 + total_chi2 += chi2_qsed
  606 + total_rchi2 += rchi2_qsed
  607 +endif
  608 +if !fit_rchi2_weight.used ne 0 then begin
  609 + total_n += n_used
  610 + total_chi2 += chi2_used
  611 + total_rchi2 += rchi2_used
356 612 endif
  613 +if !fit_rchi2_weight.qext ne 0 then begin
  614 + total_n += n_qext
  615 + total_chi2 += chi2_qext
  616 + total_rchi2 += rchi2_qext
  617 +endif
  618 +if !fit_rchi2_weight.uext ne 0 then begin
  619 + total_n += n_uext
  620 + total_chi2 += chi2_uext
  621 + total_rchi2 += rchi2_uext
  622 +endif
  623 +
  624 +endif
  625 +
357 626  
358 627 total_wchi2 = chi2_ext * !fit_rchi2_weight.ext $
359 628 + chi2_sed * !fit_rchi2_weight.sed
360 629 if !run_pol then begin
361 630 total_wchi2 += chi2_polext * !fit_rchi2_weight.polext $
362 631 + chi2_polsed * !fit_rchi2_weight.polsed $
363   - + chi2_polfrac * !fit_rchi2_weight.polfrac
  632 + + chi2_polfrac * !fit_rchi2_weight.polfrac $
  633 + + chi2_qsed * !fit_rchi2_weight.qsed $
  634 + + chi2_used * !fit_rchi2_weight.used $
  635 + + chi2_qext * !fit_rchi2_weight.qext $
  636 + + chi2_uext * !fit_rchi2_weight.uext
364 637 endif
365 638  
366 639 message,'======================================================================',/continue
... ... @@ -373,7 +646,7 @@ if !run_pol then print,&quot;Weights&quot;,!fit_rchi2_weight.ext,!fit_rchi2_weight.sed,!f
373 646 else print,"Weights",!fit_rchi2_weight.ext,!fit_rchi2_weight.sed
374 647 print,"Mean weighted reduced chi2 = ", total_wchi2/(total_n-DOF)
375 648 print,"Unweighted reduced Total chi2 = ", total_chi2/(total_n-DOF)
376   -print,"Unweighted mean reduced chi2 = ", total_rchi2/4*total_n/(total_n-DOF)
  649 +print,"Unweighted mean reduced chi2 = ", total_rchi2/4*total_n/(total_n-DOF);total_rchi2/4*total_n/(total_n-DOF) original
377 650 message,'======================================================================',/continue
378 651  
379 652 ;print,"Total reduced chi2 = ", (wrchi2_sed+wrchi2_ext+wrchi2_polext+wrchi2_polsed+wrchi2_polfrac)/(n_sed+n_ext+n_polext+n_polsed+n_polfrac-DOF)
... ...
src/idl/dustem_mpfit_sed.pro 100644 → 100755
src/idl/dustem_nir_continuum.pro 100644 → 100755
... ... @@ -55,6 +55,7 @@ ENDIF
55 55  
56 56 scope=['ADD_SED']
57 57  
  58 +
58 59 lambir=((*!dustem_params).lambda.lambda)
59 60  
60 61 ;=== normalize spectrum to the requested amplitude.
... ...
src/idl/dustem_parameter_description2type.pro 100644 → 100755
src/idl/dustem_planck_function.pro 100644 → 100755
src/idl/dustem_plot_data_sed.pro 100644 → 100755
src/idl/dustem_plot_ext.pro 100644 → 100755
... ... @@ -172,8 +172,19 @@ ENDIF ELSE BEGIN
172 172 cgoplot,xrange,(st.ext.abs_grain(i)+st.ext.sca_grain(i))/norm,color=colors(i+1),psym=psym,line=ls(i+1)
173 173 ENDELSE
174 174 ENDFOR
175   - ;oplot,xrange,(abs_tot+sca_tot)/norm,psym=psym
176   - cgoplot,xrange,(abs_tot+sca_tot)/norm,psym=psym
  175 + ;cgoplot,xrange,(abs_tot+sca_tot)/norm,psym=psym
  176 +
  177 + ;cgoplot,xrange,(abs_tot+sca_tot+powerlaw)/norm,psym=psym
  178 +
  179 + cgoplot,xrange,(st.ext.ext_tot)/norm,psym=psym
  180 +
  181 +
  182 + ; LINES OF CODE TO PLOT THE NEW ADDED PLUGINS TO THE EXTINCTION WILL HAVE TO BE ADDED HERE
  183 +
  184 + cgoplot,xrange,(*(*!dustem_plugin).powerlaw)/norm,psym=psym,color='cadetblue'
  185 +
  186 +
  187 +
177 188  
178 189 IF !dustem_show_plot EQ 2 or multi EQ 2 THEN BEGIN
179 190 norm = interpol(st.ext.ext_tot,st.ext.wav,(*!dustem_data.ext).wav)
... ...
src/idl/dustem_plot_ext_notgood.pro 100644 → 100755
src/idl/dustem_plot_extinction.pro 100644 → 100755
src/idl/dustem_plot_fit_extsed.pro 100644 → 100755
src/idl/dustem_plot_fit_sed.pro 100644 → 100755
1   -PRO dustem_plot_fit_sed,st,sed,cont,freefree,synchrotron, $
2   - res=res,errors=errors,chi2=chi2,rchi2=rchi2, $
3   - no_spec_error=no_spec_error, $
4   - help=help,win=win,ps=ps, $
5   - legend_xpos=legend_xpos,legend_ypos=legend_ypos,legend_offset=legend_offset, $
6   - _extra=_extra
  1 + PRO dustem_plot_fit_sed,st,sed, $
  2 + res=res,errors=errors,chi2=chi2,rchi2=rchi2, $
  3 + no_spec_error=no_spec_error, $
  4 + help=help,win=win,ps=ps, $
  5 + legend_xpos=legend_xpos,legend_ypos=legend_ypos,legend_offset=legend_offset, $
  6 + _extra=_extra
7 7  
8 8 ;+
9 9 ; NAME:
... ... @@ -55,12 +55,13 @@ IF keyword_set(help) THEN BEGIN
55 55 ENDIF
56 56  
57 57 IF keyword_set(ps) THEN BEGIN
58   - set_plot, 'PS'
59   - device, filename=ps, /color, /encapsulated
60   -ENDIF ;ELSE BEGIN
61   -; set_plot,'X'
62   -; IF keyword_set(win) then window,win
63   -;ENDELSE
  58 + set_plot, 'PS'
  59 + device, filename=ps, /color, /encapsulated
  60 +ENDIF ELSE BEGIN
  61 +
  62 + set_plot,'X'
  63 + IF keyword_set(win) then window,win
  64 +ENDELSE
64 65  
65 66 fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/st.sed.wav)*1.e20/1.e7
66 67  
... ... @@ -86,47 +87,17 @@ IF not keyword_set(line_tot) THEN BEGIN
86 87 ENDIF ELSE BEGIN
87 88 use_line_tot=line_tot
88 89 ENDELSE
89   -;color and linestyle for continuum
90   -IF not keyword_set(col_cont) THEN BEGIN
91   - ;use_col_cont=30
92   - use_col_cont='Navy'
93   -ENDIF ELSE BEGIN
94   - use_col_count=col_cont
95   -ENDELSE
96   -IF not keyword_set(line_cont) THEN BEGIN
97   - use_line_cont=0
98   -ENDIF ELSE BEGIN
99   - use_line_cont=line_cont
100   -ENDELSE
101   -;color and linestyle for freefree
102   -IF not keyword_set(col_freefree) THEN BEGIN
103   - ;use_col_cont=30
104   - use_col_freefree='Dark Red'
105   -ENDIF ELSE BEGIN
106   - use_col_freefree=col_freefree
107   -ENDELSE
108   -IF not keyword_set(line_freefree) THEN BEGIN
109   - use_line_freefree=0
110   -ENDIF ELSE BEGIN
111   - use_line_freefree=line_freefree
112   -ENDELSE
113   -IF not keyword_set(col_synchrotron) THEN BEGIN
114   - ;use_col_cont=30
115   - use_col_synchrotron='Dark Red'
116   -ENDIF ELSE BEGIN
117   - use_col_synchrotron=col_synchrotron
118   -ENDELSE
119   -IF not keyword_set(line_synchrotron) THEN BEGIN
120   - use_line_synchrotron=0
121   -ENDIF ELSE BEGIN
122   - use_line_synchrotron=line_synchrotron
123   -ENDELSE
124 90  
125   -;spec = st.sed.em_tot * fact + cont
  91 +
126 92 spec = st.sed.em_tot * fact
127   -IF !dustem_idl_continuum THEN spec=spec+cont
128   -IF !dustem_idl_freefree THEN spec=spec+freefree
129   -IF !dustem_idl_synchrotron THEN spec=spec+synchrotron
  93 +
  94 +;ADDING PLUGIN(S) TO SPECTRUM----------------
  95 +;if n_tags(!dustem_data.sed) gt 1 then begin
  96 +for i=0L,n_tags(*!dustem_scope)-1 do begin
  97 +if total(strsplit((*!dustem_scope).(i),'+',/extract) eq 'ADD_SED') then spec+=(*(*!dustem_plugin).(i))
  98 +endfor
  99 +;endif
  100 +;------------------------------------------
130 101  
131 102 ;stop
132 103  
... ... @@ -141,8 +112,18 @@ use_lines=replicate(0,Ngrains)
141 112  
142 113 norm = sed * 0. + 1
143 114  
  115 +xr=[1,10000]
  116 +yr=[1e-5,1e3]
  117 +
144 118 ;====== PLOT THE SED
145   -cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/norm,/nodata,_extra=_extra
  119 +cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/norm,/nodata,xtit='',ytit=textoidl('Brightness/N_H (MJy/sr/H)'),tit='Spectral Energy Distribution',ylog=1,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.93,0.93],xtickformat='(A1)'
  120 +
  121 +
  122 +;cgplot,st.polsed.wav,Q_sed*fact,xtit='',ytit=ytitstq,tit=titstq,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.95,0.95],xtickformat='(A1)'
  123 +
  124 +;this is where you will add the normalized sed
  125 +
  126 +
146 127 ind_filt=where((*!dustem_data.sed).filt_names NE 'SPECTRUM',count_filt)
147 128 ind_spec=where((*!dustem_data.sed).filt_names EQ 'SPECTRUM',count_spec)
148 129 ;=== Plot the data
... ... @@ -151,7 +132,7 @@ IF count_spec NE 0 THEN BEGIN
151 132 xx=((*!dustem_data.sed).wav)[ind_spec]
152 133 yy=((*!dustem_data.sed).values)[ind_spec]/norm[ind_spec]
153 134 rms=3.*((*!dustem_data.sed).sigma)[ind_spec]/2./norm[ind_spec]
154   - cgoplot,xx,yy,psym=8,_extra=extra,syms=0.5,color=use_col_sed_spec
  135 + cgoplot,xx,yy,psym=8,syms=0.5,color=use_col_sed_spec
155 136 IF not keyword_set(no_spec_error) THEN BEGIN
156 137 cgerrplot,xx,yy-rms,yy+rms,color=use_col_sed_spec
157 138 ;err_bar,((*!dustem_data.sed).wav)(ind_spec),((*!dustem_data.sed).values)(ind_spec)/norm(ind_spec),yrms=3.*((*!dustem_data.sed).sigma)(ind_spec)/2./norm(ind_spec)
... ... @@ -163,9 +144,9 @@ IF count_filt NE 0 THEN BEGIN
163 144 yy=((*!dustem_data.sed).values)[ind_filt]/norm[ind_filt]
164 145 rms=3.*((*!dustem_data.sed).sigma)[ind_filt]/2./norm[ind_filt]
165 146 plotsym,0,/fill
166   - cgoplot,xx,yy,psym=8,_extra=extra,color=use_col_data_filt
  147 + cgoplot,xx,yy,psym=8,color='Dodger Blue';use_col_data_filt
167 148 ;err_bar,((*!dustem_data.sed).wav)(ind_filt),((*!dustem_data.sed).values)(ind_filt)/norm(ind_filt),yrms=3.*((*!dustem_data.sed).sigma)(ind_filt)/2./norm(ind_filt),color=use_col_data_filt
168   - cgerrplot,xx,yy-rms,yy+rms,color=use_col_data_filt
  149 + cgerrplot,xx,yy-rms,yy+rms,color='Dodger Blue';use_col_data_filt
169 150 ENDIF
170 151 ;=== Plot the computed SED
171 152 IF count_filt NE 0 THEN BEGIN
... ... @@ -190,26 +171,32 @@ ENDELSE
190 171 FOR i=0L,Ngrains-1 DO BEGIN
191 172 cgoplot,st.sed.wav,st.sed.(i+1)*fact/norm,color=use_cols[i],linestyle=use_lines[i]
192 173 ENDFOR
193   -;NIR continuum
194   -IF !dustem_idl_continuum THEN BEGIN
195   - cgoplot,st.sed.wav,cont,color=use_col_cont,linestyle=use_line_cont
  174 +
  175 +
  176 +;PLOTTING OF THE PLUGIN(S)---------------
  177 +IF tag_exist(*!dustem_scope,'CONTINUUM') THEN BEGIN
  178 + cgoplot,st.sed.wav,(*(*!dustem_plugin).continuum),color='Navy',linestyle=3
196 179 ENDIF
197   -IF !dustem_idl_freefree THEN BEGIN
198   - cgoplot,st.sed.wav,freefree,color=use_col_freefree,linestyle=use_line_freefree
  180 +IF tag_exist(*!dustem_scope,'FREE_FREE') THEN BEGIN
  181 + cgoplot,st.sed.wav,(*(*!dustem_plugin).freefree),color='Dark Red',linestyle=3
199 182 ENDIF
200   -IF !dustem_idl_synchrotron THEN BEGIN
201   - cgoplot,st.sed.wav,synchrotron,color=use_col_synchrotron,linestyle=use_line_synchrotron
  183 +IF tag_exist(*!dustem_scope,'SYNCHROTRON') THEN BEGIN
  184 + cgoplot,st.sed.wav,(*(*!dustem_plugin).synchrotron),color='Crimson',linestyle=3
202 185 ENDIF
  186 +
203 187 cgoplot,st.sed.wav,spec/norm,color=use_col_tot,linestyle=use_line_tot
  188 +;plot the normalized data as well.
  189 +;----------------------------------------
204 190  
205 191 ;==== print the legend
206 192 frmt0='(A26)'
207 193 frmt1='(1E10.2)'
208 194 frmt2='(F7.2)'
209   -use_legend_xpos=0.8
210   -use_legend_ypos=0.8
  195 +use_legend_xpos=0.59
  196 +use_legend_ypos=0.92
211 197 use_legend_offset=0.03 ;This is the offset between lines of the legend in normalized units
212 198 legend_charsize=1.0
  199 +k=0
213 200 IF keyword_set(legend_xpos) THEN use_legend_xpos=legend_xpos
214 201 IF keyword_set(legend_ypos) THEN use_legend_ypos=legend_ypos
215 202 IF keyword_set(legend_offset) THEN use_legend_offset=legend_offset
... ... @@ -226,20 +213,40 @@ IF keyword_set(res) THEN BEGIN
226 213 IF keyword_set(errors) THEN BEGIN
227 214 str=str+textoidl('\pm')+string(errors(i),format=frmt1)
228 215 ENDIF
229   - yypos=use_legend_ypos-(i+1)*use_legend_offset
230   - xyouts,use_legend_xpos,yypos,str,color=0,/normal,charsize=legend_charsize
  216 +
  217 + if STRUPCASE(strmid(strtrim(parameter_description,2),0,6)) eq 'DUSTEM' then begin
  218 + yypos=use_legend_ypos-(i-k)*use_legend_offset
  219 + xxpos=use_legend_xpos*0.52
  220 + endif else begin
  221 + yypos=use_legend_ypos-(i+1)*use_legend_offset
  222 + xxpos=use_legend_xpos
  223 + k=i
  224 + endelse
  225 + ;yypos=use_legend_ypos-(i+1)*use_legend_offset
  226 + xyouts,xxpos,yypos,str,color=0,/normal,charsize=legend_charsize
231 227 ENDFOR
232 228 ;stop
233 229 ENDIF
234 230 IF keyword_set(chi2) THEN BEGIN
235   - yypos=use_legend_ypos-(Npar+1)*use_legend_offset
236   - xyouts,use_legend_xpos,yypos,string('chi2=',format=frmt0)+string(chi2,format=frmt2),color=0,/normal,charsize=legend_charsize
  231 + xxpos=use_legend_xpos*0.11
  232 + yypos=use_legend_ypos-(1)*use_legend_offset
  233 + xyouts,xxpos,yypos,string('chi2=',format=frmt0)+string(chi2,format=frmt2),color=0,/normal,charsize=legend_charsize
237 234 ENDIF
238 235 IF keyword_set(rchi2) THEN BEGIN
239   - yypos=use_legend_ypos-(Npar+2)*use_legend_offset
240   - xyouts,use_legend_xpos,yypos,string('red. chi2=',format=frmt0)+string(rchi2,format=frmt2),color=0,/normal,charsize=legend_charsize
  236 + xxpos=use_legend_xpos*0.11
  237 + yypos=use_legend_ypos-(2)*use_legend_offset
  238 + xyouts,xxpos,yypos,string('red. chi2=',format=frmt0)+string(rchi2,format=frmt2),color=0,/normal,charsize=legend_charsize
241 239 ENDIF
242 240  
  241 +xtit=textoidl('\lambda (\mum)')
  242 +cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/sed,/nodata,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.17,0.14,0.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1
  243 +;plot the normalized data as well.
  244 +cgoplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/sed,psym=16,symsize=1,thick=2,color='Dodger Blue'
  245 +cgoplot,st.sed.wav,spec/spec,color='black'
  246 +cgerrplot,((*!dustem_data.sed).wav),((*!dustem_data.sed).values)/sed-3.*((*!dustem_data.sed).sigma)/2./sed,((*!dustem_data.sed).values)/sed+3.*((*!dustem_data.sed).sigma)/2./sed,color='Dodger Blue'
  247 +
  248 +
  249 +
243 250 IF keyword_set(ps) THEN BEGIN
244 251 device,/close
245 252 set_plot,'X'
... ...
src/idl/dustem_plot_fpol.pro 100644 → 100755
... ... @@ -6,7 +6,7 @@ if not keyword_set(multi) then multi = 0
6 6  
7 7 nsz_max = 200
8 8 ; Read DUSTEM_wrap fit files
9   -if not keyword_set(dir_in) then dir_in = !dustem_dat
  9 +if not keyword_set(dir_in) then dir_in = !dustem_soft_dir ;!dustem_dat because it was looking for a file only there. My might need to be revised I'm not sure
10 10  
11 11 if multi ne 2 then begin
12 12 if keyword_set(ps) then begin
... ... @@ -18,11 +18,14 @@ endif else begin
18 18 endelse
19 19 endif
20 20  
21   -;st_grains=dustem_read_grain(dir_in+'data/'+(*!dustem_inputs).grain,/silent)
22   -st_grains=dustem_read_grain(dir_in+'data/'+'GRAIN.DAT',/silent)
  21 +st_grains=dustem_read_grain(dir_in+'data/'+(*!dustem_inputs).grain,/silent)
  22 +;st_grains=dustem_read_grain(dir_in+'data/'+'GRAIN.DAT',/silent)
23 23  
24   -ngrains = st_grains.ngrains
25   -colors=dustem_grains_colors(Ngrains,ls,ps=ps)
  24 +;ngrains = st_grains.ngrains I WILL BE REPLACING THIS LINE
  25 +ngrains= ((*!dustem_params).ngrains)
  26 +
  27 +
  28 +colors=dustem_grains_colors(ngrains,ls,ps=ps)
26 29  
27 30 ;j=where(polfile ne '', count)
28 31 ;if count le 1 then begin
... ... @@ -41,15 +44,15 @@ amax = 1e8
41 44 nrad = 20000
42 45 a = amin * exp(alog(amax/amin)/(nrad-1)*indgen(nrad))
43 46  
44   -if multi eq 0 then plot,a,a,/nodata,xr=xr,yr=yr,/xs,/ys,xtit=xtit,ytit=ytit,tit='',/xlog
45   -if multi eq 2 then plot,a,a,/nodata,xr=xr,yr=yr,/xs,/ys,xtit=xtit,ytit=ytit,tit='',/xlog,position=[0.17,0.14,0.95,0.45],/noerase,ymino=2,xticklen=0.07
  47 +if multi eq 0 then cgplot,a,a,/nodata,xr=xr,yr=yr,/xs,/ys,xtit=xtit,ytit=ytit,tit=tit,/xlog
  48 +if multi eq 2 then cgplot,a,a,/nodata,xr=xr,yr=yr,/xs,/ys,xtit=xtit,ytit=ytit,tit='',/xlog,position=[0.17,0.14,0.95,0.45],/noerase,ymino=2,xticklen=0.07
46 49  
47 50 st_align=dustem_read_align(dir_in+'data/',st_grains,silent=silent)
48 51  
49 52 for i=0,ngrains-1 do begin
50 53 if st_align.grains(i).aligned then begin
51 54 f_pol = 0.5*st_align.grains(i).plev * (1 + TANH((alog(a*1e-3/(st_align.grains(i).athresh))) / st_align.grains(i).pstiff ) )
52   - oplot,a,f_pol, line=ls(i+1),col=colors(i+1)
  55 + cgoplot,a,f_pol, line=ls(i+1),color='black';,col=colors(i+1) ; this line should be tested with multiple polarizing grain species
53 56 endif
54 57 endfor
55 58  
... ...
src/idl/dustem_plot_legend.pro 100644 → 100755
src/idl/dustem_plot_nuinu_em.pro 100644 → 100755
src/idl/dustem_plot_polar.pro 100644 → 100755
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
  1 +PRO dustem_plot_polar,st,p_dim,ps=ps,_Extra=extra,help=help,UV=UV,SED=SED,Pfrac=Pfrac,print_ratio=print_ratio,win=win,xr=xr,yr=yr,donotclose=donotclose,aligned=aligned,noerrbars=noerrbars,multi=multi,almabands=almabands,nsmooth=nsmooth
2 2  
3 3 ;+
4 4 ; NAME:
... ... @@ -61,7 +61,9 @@ if multi ne 2 then begin
61 61 device, filename=ps, /color, /encapsulated
62 62 endif else begin
63 63 set_plot,'X'
64   - if keyword_set(win) then window,win
  64 + if keyword_set(win) then begin
  65 + if ptr_valid(!dustem_data.polsed) then window,win,title='DUSTEM WRAP (POLSED)' else window,win,title='DUSTEM WRAP (POLEXT)'
  66 + endif
65 67 endelse
66 68 endif
67 69  
... ... @@ -70,12 +72,11 @@ endif
70 72 use_col_sed_filt = 70
71 73 use_col_data_filt = 250
72 74  
73   -;Ngrains=n_tags(st.sed)-2
74   -Ngrains=(*!dustem_params).Ngrains
  75 +Ngrains=n_tags(st.sed)-2
75 76  
76 77 colors=dustem_grains_colors(Ngrains,ls,ps=ps)
77 78  
78   -xtit=textoidl('\lambda (\mum)')
  79 +if ptr_valid(!dustem_data.polsed) then xtit=textoidl('\lambda (\mum)') else xtit=textoidl('\lambda^{-1} (\mum^{-1})')
79 80  
80 81 ;ext_tot=st.ext.abs_grain(0)*0.
81 82 ;polext_tot=st.polext.abs_grain(0)*0.
... ... @@ -88,6 +89,7 @@ xtit=textoidl(&#39;\lambda (\mum)&#39;)
88 89 defsysv,'!dustem_data',exists=data_exists
89 90  
90 91 IF ptr_valid(!dustem_data.polext) THEN BEGIN
  92 + interpolext = interpol(st.polext.ext_tot,st.polext.wav,(*!dustem_data.polext).wav)
91 93 IF keyword_set(print_ratio) THEN BEGIN
92 94 ; Correction de couleur
93 95 cc = 1.11
... ... @@ -123,12 +125,12 @@ IF ptr_valid(!dustem_data.polext) THEN BEGIN
123 125 ; 1 W/m2/sr/1d20H = 1e-17 erg/s/sr/H
124 126 ENDIF
125 127 IF ptr_valid(!dustem_data.polsed) THEN BEGIN
126   - dustem_polsed = dustem_compute_polsed(p_dim,st.polsed,cont=cont)
  128 + dustem_polsed = dustem_compute_polsed(p_dim,st);dustem_compute_polsed(p_dim,st.polsed,cont=cont)
127 129 x=min(((*!dustem_data.polsed).wav-850)^2,jmin)
128 130 P353 = dustem_polsed(jmin)
129 131 P353_data = (*!dustem_data.polsed).values(jmin)
130 132  
131   - dustem_sed = dustem_compute_sed(p_dim,st,cont=cont)
  133 + dustem_sed = dustem_compute_sed(p_dim,st);dustem_compute_sed(p_dim,st,cont=cont)
132 134 x=min(((*!dustem_data.sed).wav-850)^2,jmin)
133 135 I353 = dustem_sed(jmin)
134 136 I353_data = (*!dustem_data.sed).values(jmin)
... ... @@ -189,16 +191,15 @@ IF ptr_valid(!dustem_data.polext) THEN BEGIN
189 191 ;print
190 192 ENDIF
191 193 ENDIF
192   -ENDIF ELSE BEGIN
  194 +
193 195 IF keyword_set(UV) THEN BEGIN
194   - if not keyword_set(ps) then tit='Dustem polarization cross-section'
  196 + if not keyword_set(ps) then tit='Polarization Cross-section'
195 197 ytit=textoidl('\sigma_{pol}/N_H (cm^2/H)')
196 198 if ptr_valid(!dustem_data.polext) then begin
197   -
198 199 if !dustem_show_plot EQ 2 or multi eq 2 then begin
199 200 ;normpol = polext_tot
200 201 normpol = st.polext.ext_tot
201   - dustem_polext = interpol(st.polext.ext_tot,st.polext.wav,(*!dustem_data.polext).wav)
  202 + dustem_polext = interpol(st.polext.ext_tot,st.polext.wav,(*!dustem_data.polext).wav)
202 203 yr=[0.,2.]
203 204 ylog = 0
204 205 endif else begin
... ... @@ -207,33 +208,51 @@ ENDIF ELSE BEGIN
207 208 ylog = 1
208 209 endelse
209 210  
210   - 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
211   - 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)'
212   - 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
213   -
214   -
  211 + if multi eq 0 then cgplot,st.polext.wav,st.polext.ext_tot/normpol,/nodata,xtit=xtit,ytit=ytit,tit=tit,xr=xr,yr=yr,ylog=ylog,xlog=1
  212 + if multi eq 1 then cgplot,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.93,0.93],xtickformat='(A1)' ,xlog=1
  213 + if multi eq 2 then cgplot,st.polext.wav,st.polext.ext_tot/normpol,/nodata,_Extra=extra,xtit=xtit,ytit='Normalized',tit='',xr=xr,yr=[0,2],ylog=ylog,position=[0.17,0.14,0.93,0.35],/noerase,yticks=3,ymino=5,xticklen=0.1,xlog=1
  214 +
  215 +
215 216 FOR i=0L,Ngrains-1 DO BEGIN
216   - 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)
  217 + if aligned(i) then cgoplot,1/st.polext.wav,(st.polext.abs_grain(i)+st.polext.sca_grain(i))/normpol,color=colors(i+1),psym=psym,line=ls(i+1)
217 218 ENDFOR
218 219  
219   - if total(aligned gt 0) then oplot,st.polext.wav,st.polext.ext_tot/normpol,psym=psym
220   -
  220 + if total(aligned gt 0) then cgoplot,1/st.polext.wav,st.polext.ext_tot/normpol,psym=psym
  221 +
221 222 if data_exists then BEGIN
222 223 if ptr_valid(!dustem_data.polext) then begin
223 224 plotsym,0,/FILL
224   - oplot, (*!dustem_data.polext).wav,(*!dustem_data.polext).values/dustem_polext,ps=4
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
  225 + cgoplot,1/(*!dustem_data.polext).wav,(*!dustem_data.polext).values/dustem_polext,psym=16,color='Olive Drab'
  226 + if multi ne 2 then cgoplot,1/(*!dustem_data.polext).wav,interpolext/dustem_polext,psym=6,syms=2,color='red'
  227 + 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,col=cgColor('Olive Drab')
226 228 endif
227 229 ENDIF
228 230 ENDIF
229 231 ENDIF
230   -ENDELSE
231 232  
  233 + if multi eq 0 then begin
  234 + if not keyword_set(ps) then tit='Polarization fraction in extinction'
  235 + ytit='Polarization fraction'
  236 + cgplot,st.polext.wav,st.polext.ext_tot/st.ext.ext_tot,/nodata,_Extra=extra,xtit=xtit,ytit=ytit,tit=tit,xr=xr,yr=yr
  237 + FOR i=0L,Ngrains-1 DO BEGIN
  238 + cgoplot,st.polext.wav,smooth((st.polext.abs_grain(i)+st.polext.sca_grain(i))/(st.ext.abs_grain(i)+st.ext.sca_grain(i)),nsmooth),color=colors(i+1),psym=psym,line=ls(i+1)
  239 + ENDFOR
  240 + if total(aligned gt 0) then cgoplot,st.polext.wav,st.polext.ext_tot/st.ext.ext_tot,psym=psym
  241 + cgoplot,st.polsed.wav,st.polsed.(Ngrains+1)/st.sed.(Ngrains+1),psym=psym,linestyle=2
  242 +
  243 + if data_exists and ptr_valid(!dustem_data.polext) and ptr_valid(!dustem_data.ext) then begin
  244 + values_ext_for_polext = interpol((*!dustem_data.ext).values,(*!dustem_data.ext).wav,(*!dustem_data.polext).wav)
  245 + sigma_ext_for_polext = interpol((*!dustem_data.ext).sigma ,(*!dustem_data.ext).wav,(*!dustem_data.polext).wav)
  246 + if not keyword_set(noerrbars) then err_bar,(*!dustem_data.polext).wav,(*!dustem_data.polext).values/values_ext_for_polext, $
  247 + yrms=1./values_ext_for_polext * (*!dustem_data.polext).sigma + (*!dustem_data.polext).values/values_ext_for_polext^2 * sigma_ext_for_polext
  248 + endif
  249 + endif
  250 +ENDIF
232 251 IF ptr_valid(!dustem_data.polsed) THEN BEGIN
233   - if not keyword_set(ps) then tit='Polarization intensity in emission'
234   - ytit='Polarization intensity'
  252 + if not keyword_set(ps) then tit='Polarization Intensity';' in emission'
  253 + ytit='Polarization Intensity' ;is this line used?
235 254 ytit=textoidl('\nuP_\nu/N_H (W/m^2/sr/H)')
236   -
  255 +
237 256 ; Convert from erg/s/H into W/m2/sr/H
238 257 fact=1./(4.*!pi)*1.e17/1d20
239 258  
... ... @@ -250,18 +269,22 @@ IF ptr_valid(!dustem_data.polsed) THEN BEGIN
250 269 if not keyword_set(yr) then yr=[1e-34,1e-30]
251 270 endelse
252 271  
253   - 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
254   - 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)'
255   - 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
256 272  
  273 + 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
  274 + 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.93,0.93],xtickformat='(A1)'
  275 + 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.93,0.35],/noerase,yticks=2,ymino=5,xticklen=0.1
  276 +
  277 +
257 278 for i = 0, Ngrains-1 do begin
258   - if aligned(i) then cgoplot,st.polsed.wav,st.polsed.(i+1)*fact/normpol,color=colors(i+1),psym=psym,line=ls(i+1)
  279 + if aligned(i) then cgoplot,st.polsed.wav,st.polsed.(i+1)*fact/normpol,color=colors(i+1),psym=psym,line=ls(i+1) ;investigate this line of the plotting of the different components in polarization
259 280 endfor
260   - if (total(aligned) gt 0) then cgoplot,st.polsed.wav,st.polsed.(Ngrains+1)*fact/normpol
261   -
  281 +
  282 + if total(aligned gt 0) && multi ne 2 then cgoplot,st.polsed.wav,st.polsed.(Ngrains+1)*fact/normpol,color='Teal'
  283 + if total(aligned gt 0) && multi eq 2 then cgoplot,st.polsed.wav,st.polsed.(Ngrains+1)*fact/normpol,color='black'
  284 +
262 285 if data_exists then BEGIN
263 286 if ptr_valid(!dustem_data.polsed) then begin
264   -
  287 +
265 288 ind_filt=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_filt)
266 289 ;=== Plot the data
267 290 IF count_filt NE 0 THEN BEGIN
... ... @@ -272,24 +295,42 @@ IF ptr_valid(!dustem_data.polsed) THEN BEGIN
272 295 ; Convert data which Pnu in MJy/sr/1d20 into nu*Pnu in W/m2/sr/1d20
273 296 fact = (3.e10/(1.e-4*(*!dustem_data.polsed).wav(ind_filt)))/1.d40
274 297 endelse
275   - dustem_polsed = dustem_compute_polsed(p_dim,st,cont=cont)*fact
  298 + dustem_polsed = dustem_compute_polsed(p_dim,st)*fact;dustem_compute_polsed(p_dim,st,cont=cont)*fact
  299 +
276 300 if !dustem_show_plot EQ 2 or multi eq 2 then begin
277 301 normpol = dustem_polsed
278 302 endif else begin
279 303 normpol = dustem_polsed * 0. + 1
280 304 endelse
281 305 plotsym,0,/FILL
282   - 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
283   - 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
  306 +
  307 + cgoplot,((*!dustem_data.polsed).wav)(ind_filt),((*!dustem_data.polsed).values)(ind_filt)*fact/normpol(ind_filt),psym=8,_extra=extra,color='Orchid',thick=2;use_col_sed_filt
  308 + if multi ne 2 then cgoplot,((*!dustem_data.polsed).wav),dustem_polsed/normpol,psym=6,color='red',symsize=2
  309 +
  310 +
  311 + 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=cgColor('Orchid');use_col_sed_filt
284 312  
285 313 ; PLot color corrected POLSED
286 314 if keyword_set(ps) then plotsym,8,thick=8 else plotsym,8
287 315 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='red'
288 316 ENDIF
  317 +
  318 + IF tag_exist(*!dustem_scope,'SYNCHROTRON') THEN BEGIN
  319 + indx=where((*(*!dustem_fit).param_descs) EQ 'dustem_create_stokes_3',counttx)
  320 +
  321 + if total(counttx) then frac_synch=p_dim(indx) else frac_synch=0.7 ;THIS LINE IS TO TAKE INTO ACCOUNT WHEN THE AMPLITUDE OF THE SYNCHROTRON IS FROZEN
  322 +
  323 + cgoplot,st.polsed.wav,(*(*!dustem_plugin).synchrotron)*frac_synch(0)*(3.e10/(1.e-4*(st.polsed.wav)))/1.d40,color='Crimson',linestyle=3
  324 + cgoplot,st.polsed.wav,st.polsed.em_tot*1./(4.*!pi)*1.e17/1d20+((*(*!dustem_plugin).synchrotron))*frac_synch(0)*(3.e10/(1.e-4*(st.polsed.wav)))/1.d40,color='black'
  325 +
  326 + ENDIF
  327 +
289 328 endif
290   - endif
291   -
292   -endif else BEGIN
  329 +
  330 +
  331 +
  332 + endif
  333 +endif ;else BEGIN
293 334 if keyword_set(Pfrac) then begin
294 335  
295 336 ; Also show prediction for polarized SED
... ... @@ -324,7 +365,7 @@ endif else BEGIN
324 365 if keyword_set(ps) then plotsym,8,thick=8 else plotsym,8
325 366  
326 367 dustem_polsed = dustem_compute_polsed(p_dim,st.polsed)
327   - dustem_sed_tmp = dustem_compute_sed(p_dim,st,cont=cont)
  368 + dustem_sed_tmp = dustem_compute_sed(p_dim,st);dustem_compute_sed(p_dim,st,cont=cont)
328 369 dustem_sed = dustem_polsed
329 370 for i = 0, count_filt-1 do begin
330 371 j=where((*!dustem_data.polsed).filt_names(ind_filt(i)) eq (*!dustem_data.sed).filt_names,count)
... ... @@ -335,32 +376,18 @@ endif else BEGIN
335 376  
336 377 endif
337 378  
338   - endif else begin
339   -
340   - if not keyword_set(ps) then tit='Polarization fraction in extinction'
341   - ytit='Polarization fraction'
342   - plot_oi,st.polext.wav,st.polext.ext_tot/st.ext.ext_tot,/nodata,_Extra=extra,xtit=xtit,ytit=ytit,tit=tit,xr=xr,yr=yr
343   - FOR i=0L,Ngrains-1 DO BEGIN
344   - oplot,st.polext.wav,smooth((st.polext.abs_grain(i)+st.polext.sca_grain(i))/(st.ext.abs_grain(i)+st.ext.sca_grain(i)),nsmooth),color=colors(i+1),psym=psym,line=ls(i+1)
345   - ENDFOR
346   - if (total(aligned) gt 0) then oplot,st.polext.wav,st.polext.ext_tot/st.ext.ext_tot,psym=psym
347   - oplot,st.polsed.wav,st.polsed.(Ngrains+1)/st.sed.(Ngrains+1),psym=psym,linestyle=2
  379 + endif ;else begin
348 380  
349   - if data_exists and ptr_valid(!dustem_data.polext) and ptr_valid(!dustem_data.ext) then begin
350   - values_ext_for_polext = interpol((*!dustem_data.ext).values,(*!dustem_data.ext).wav,(*!dustem_data.polext).wav)
351   - sigma_ext_for_polext = interpol((*!dustem_data.ext).sigma ,(*!dustem_data.ext).wav,(*!dustem_data.polext).wav)
352   - if not keyword_set(noerrbars) then err_bar,(*!dustem_data.polext).wav,(*!dustem_data.polext).values/values_ext_for_polext, $
353   - yrms=1./values_ext_for_polext * (*!dustem_data.polext).sigma + (*!dustem_data.polext).values/values_ext_for_polext^2 * sigma_ext_for_polext
354   - endif
355 381  
356   - endelse
357   -endelse
  382 +; endelse
  383 +; endelse
358 384  
359   -if keyword_set(ps) and not keyword_set(donotclose) then begin
  385 +if keyword_set(ps) && not keyword_set(donotclose) then begin
360 386 device,/close
361 387 set_plot,'X'
362 388 endif
363 389  
  390 +
364 391 the_end:
365 392  
366 393 END
... ...
src/idl/dustem_plot_sdist.pro 100644 → 100755
src/idl/dustem_radiance.pro 100644 → 100755
src/idl/dustem_rap_xisrf.pro 100644 → 100755
src/idl/dustem_read_align.pro 100644 → 100755
src/idl/dustem_read_all.pro 100644 → 100755
... ... @@ -14,7 +14,7 @@ FUNCTION dustem_read_all,dir_in,silent=silent,help=help
14 14 ; OPTIONAL INPUT PARAMETERS:
15 15 ; None
16 16 ; OUTPUTS:
17   -; st: input structure
  17 +; st: input structures
18 18 ; OPTIONAL OUTPUT PARAMETERS:
19 19 ; None
20 20 ; ACCEPTED KEY-WORDS:
... ... @@ -45,7 +45,7 @@ IF keyword_set(help) THEN BEGIN
45 45 ENDIF
46 46  
47 47 ;IF getenv('DUSTEM_WHICH') EQ 'WEB3p8' THEN BEGIN
48   -IF !dustem_which EQ 'WEB3p8' THEN BEGIN
  48 +IF !dustem_which EQ 'WEB3p8' THEN BEGIN ; I CHANGED THIS AS WELL
49 49 st=dustem_read_all_web3p8(dir_in,silent=silent)
50 50 goto,ionfrac
51 51 ENDIF
... ... @@ -62,7 +62,7 @@ st_grains=dustem_read_grain(file,silent=silent)
62 62  
63 63 ;stop
64 64  
65   -file=dir_in+'TAILLE.DAT'
  65 +file=dir_in+'TAILLE.DAT' ;I put this back to its original form
66 66 st_taille=dustem_read_taille(file,silent=silent)
67 67  
68 68 file=dir_in+'PROPMASS.DAT'
... ...
src/idl/dustem_read_all_lv.pro 100644 → 100755
src/idl/dustem_read_all_res.pro 100644 → 100755
src/idl/dustem_read_all_web3p8.pro 100644 → 100755
src/idl/dustem_read_calor.pro 100644 → 100755
src/idl/dustem_read_calor_lv.pro 100644 → 100755
src/idl/dustem_read_chrg.pro 100644 → 100755
src/idl/dustem_read_circ_vg.pro 100644 → 100755
src/idl/dustem_read_dustem.pro 100644 → 100755
src/idl/dustem_read_dustem_lv.pro 100644 → 100755
src/idl/dustem_read_dustem_lv.pro.~1.1.~ 100644 → 100755
src/idl/dustem_read_ext_lv.pro 100644 → 100755
src/idl/dustem_read_filters.pro 100644 → 100755
src/idl/dustem_read_gas.pro 100644 → 100755
src/idl/dustem_read_gemissiv.pro 100644 → 100755
src/idl/dustem_read_grain.pro 100644 → 100755
1 1 FUNCTION dustem_read_grain,file,silent=silent,key_str=key_str,G0=G0,help=help
2 2  
  3 +
  4 +
3 5 ;+
4 6 ; NAME:
5 7 ; dustem_read_grain
... ... @@ -86,7 +88,7 @@ ENDIF
86 88 defsysv,'!run_tls',0. ; Global flag to know if TLS is being run
87 89  
88 90 CASE !dustem_which OF
89   - 'WEB3p8': BEGIN
  91 + 'WEB3p8': BEGIN
90 92 frmt='(A,D,D,D,D,D,I,A)'
91 93 ;count # of lines
92 94 openr,unit,file,/get_lun
... ...
src/idl/dustem_read_ionfrac.pro 100644 → 100755
src/idl/dustem_read_isrf.pro 100644 → 100755
src/idl/dustem_read_lambda.pro 100644 → 100755
src/idl/dustem_read_mix.pro 100644 → 100755
src/idl/dustem_read_pol.pro 100644 → 100755
src/idl/dustem_read_propmass.pro 100644 → 100755
src/idl/dustem_read_qabs.pro 100644 → 100755
src/idl/dustem_read_qabs_desert.pro 100644 → 100755
src/idl/dustem_read_qabs_lv.pro 100644 → 100755
src/idl/dustem_read_secteffir.pro 100644 → 100755
src/idl/dustem_read_secteffuv.pro 100644 → 100755