OV-GSO-DC / dustem-wrapper_idl

Commit e7938fa312161cabc983099b9c3e2f5d3c5e7793

Authored by Ilyes Choubani
1 parent 0e608856
Exists in master

Corrected02: Implementation of plugin formalism

Showing 27 changed files with 1291 additions and 310 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_activate_plugins.pro 0 → 100755
  Diff comments   View file @e7938fa
... ... @@ -0,0 +1,105 @@
  1 +PRO dustem_activate_plugins,p_min,res=res,st;dustem_create_functions,p_min,cont=cont,freefree=freefree,synchrotron=synchrotron,st,powerlaw=powerlaw,res=res
  2 +
  3 +p_dim = p_min * (*(*!dustem_fit).param_init_values)
  4 +
  5 +;pl_dim = p_min * pliv
  6 +
  7 +f=1
  8 +FOR i=0L,n_elements(p_min)-1 DO BEGIN
  9 + ;stop
  10 + parameter_type=dustem_parameter_description2type((*(*!dustem_fit).param_descs)[i],string_name=string_name)
  11 + ;IF ((*(*!dustem_fit).param_func)(i) NE 0.) THEN BEGIN ;test if this is a parameter for a pluggin
  12 + IF parameter_type EQ 'PLUGIN' THEN BEGIN ;test if this is a parameter for a plugin
  13 +
  14 + ftn = strmid((*(*!dustem_fit).param_descs)(i), 0, strlen((*(*!dustem_fit).param_descs)(i))-2) ;caution, this would not allow more than 9 parameters. -WHY?
  15 + tmp = (where(*(*!dustem_fit).param_func eq f, count))
  16 + PDO_tmp = (*(*!dustem_fit).param_descs)(tmp)
  17 + p_dim_tmp = p_dim(tmp)
  18 + index = fltarr(count) & value = fltarr(count)
  19 +
  20 + FOR k=0, count-1 DO BEGIN
  21 + index[k] = strmid(PDO_tmp[k],0,1,/REV )
  22 + value[k] = p_dim_tmp[k]
  23 +
  24 + ENDFOR
  25 +
  26 + ;ASSIGNING DATA TO THE INITIALIZED PLUGIN DATA ARRAYS--------------
  27 + if strmid(ftn,0,13) eq 'dustem_create' then begin
  28 +
  29 + k=where(tag_names(*!dustem_scope) eq strupcase(strmid(ftn,14)))
  30 +
  31 + if total(strsplit((*!dustem_scope).(k),'+',/extract) eq 'ADD_SED') or total(strsplit((*!dustem_scope).(k),'+',/extract) eq 'ADD_POLSED') then begin
  32 + str='((*!dustem_plugin).('+strtrim(k,2)+'))=ptr_new('+ftn+'(key=index,val=value)'+')' & str=str(0) ;because for some reason, this initially returned a string array and we need a string
  33 + toto=execute(str) & IF !dustem_verbose NE 0 THEN message,strupcase(strmid(ftn,7)),/info
  34 + endif
  35 +
  36 + if total(strsplit((*!dustem_scope).(k),'+',/extract) eq 'ADD_EXT') or total(strsplit((*!dustem_scope).(k),'+',/extract) eq 'ADD_POLEXT') or total(strsplit((*!dustem_scope).(k),'+',/extract) eq 'MODIFY') then begin
  37 + str='((*!dustem_plugin).('+strtrim(k,2)+'))=ptr_new('+ftn+'(st,key=index,val=value)'+')' & str=str(0)
  38 + toto=execute(str) & IF !dustem_verbose NE 0 THEN message,strupcase(strmid(ftn,7)),/info
  39 + endif
  40 +
  41 + endif else begin
  42 + str='res='+ftn+'(key=index,val=value)'
  43 + toto=execute(str) & IF !dustem_verbose NE 0 THEN message,strupcase(strmid(ftn,7)),/info
  44 + endelse
  45 + ;-----------------------------------------------------------
  46 +
  47 +
  48 + ;I don't think this bloc is used anywhere!!! I commented it out and the wrapper works just fine
  49 +; IF ftn EQ 'dustem_create_ionfrac' THEN BEGIN
  50 +; key_ionfrac=index
  51 +; val_ionfrac=value
  52 +; ENDIF
  53 +;
  54 +
  55 +
  56 +; IF ftn NE 'dustem_create_continuum' and ftn NE 'dustem_create_freefree' and ftn NE 'dustem_create_synchrotron' and ftn NE 'dustem_create_powerlaw' THEN BEGIN
  57 +; ;str='res='+ftn+'(key=index,val=value)'
  58 +; str='res='+ftn+'(key=index,val=value)'
  59 +; toto=execute(str)
  60 +; IF !dustem_verbose NE 0 THEN message,str,/info
  61 +; ENDIF ELSE BEGIN
  62 +;
  63 +; IF ftn EQ 'dustem_create_continuum' THEN BEGIN
  64 +;
  65 +; ;str='cont='+ftn+'(key=index,val=value)' ; original line
  66 +;
  67 +; str='((*!dustem_plugin).continuum)=ptr_new('+ftn+'(key=index,val=value)'+')'
  68 +;
  69 +; toto=execute(str)
  70 +; IF !dustem_verbose NE 0 THEN message,str,/info
  71 +; ENDIF
  72 +;
  73 +;
  74 +; IF ftn EQ 'dustem_create_freefree' THEN BEGIN
  75 +; ;str='freefree='+ftn+'(key=index,val=value)'
  76 +; str='((*!dustem_plugin).freefree)=ptr_new('+ftn+'(key=index,val=value)'+')'
  77 +; toto=execute(str)
  78 +; IF !dustem_verbose NE 0 THEN message,str,/info
  79 +; ;((*!dustem_plugin).freefree)=freefree
  80 +; ENDIF
  81 +;
  82 +; IF ftn EQ 'dustem_create_synchrotron' THEN BEGIN
  83 +; ;str='synchrotron='+ftn+'(key=index,val=value)'
  84 +; str='((*!dustem_plugin).synchrotron)=ptr_new('+ftn+'(key=index,val=value)'+')'
  85 +; toto=execute(str)
  86 +; IF !dustem_verbose NE 0 THEN message,str,/info
  87 +; ;((*!dustem_plugin).synchrotron)=synchrotron
  88 +; ENDIF
  89 +;
  90 +; IF ftn EQ 'dustem_create_powerlaw' THEN BEGIN
  91 +; ;str='powerlaw='+ftn+'(st,key=index,val=value)' ; the fact that st is called here makes it harder to give a common treatment to all the plugins
  92 +; str='((*!dustem_plugin).powerlaw)=ptr_new('+ftn+'(st,key=index,val=value)' +')'
  93 +; toto=execute(str)
  94 +; IF !dustem_verbose NE 0 THEN message,str,/info
  95 +; ;((*!dustem_plugin).powerlaw)=powerlaw
  96 +; ENDIF
  97 +;
  98 +
  99 +
  100 + ; ENDELSE
  101 + f=f+1 & i=i+count-1
  102 + ENDIF
  103 +ENDFOR
  104 +
  105 +END
0 106 \ No newline at end of file
... ...
src/idl/dustem_compute_polext.pro 0 → 100644
  Diff comments   View file @e7938fa
... ... @@ -0,0 +1,52 @@
  1 +FUNCTION dustem_compute_polext ,p_dim,sti,_extra=extra,out_st=out_st,dustem_qext,dustem_uext,Q_ext,U_ext
  2 +
  3 +;CREATED AS A CONSEQUENCE OF THE INCLUSION OF THE PLUGINS
  4 +
  5 +
  6 +;written to homogonize DustEM
  7 +
  8 +IF not keyword_set(sti) THEN BEGIN
  9 + sti=dustem_run(p_dim)
  10 + dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),res=res,sti
  11 +ENDIF
  12 +
  13 +
  14 + ;ADDING PLUGIN TO SPECTRUM----------------
  15 + ;if n_tags(!dustem_data.polext) gt 1 then begin
  16 + for i=0L,n_tags(*!dustem_scope)-1 do begin
  17 + if total(strsplit((*!dustem_scope).(i),'+',/extract) eq 'ADD_POLEXT') then sti.polext.ext_tot+=(*(*!dustem_plugin).(i))
  18 + endfor
  19 + ;endif
  20 + ;------------------------------------------
  21 +
  22 + ;MODIFYING VIA PLUGIN-------------- THIS BLOCK NEEDS TO BE UNDER THE DUSTEM_POLSED LINE
  23 +
  24 + for i=0L,n_tags(*!dustem_scope)-1 do begin
  25 +
  26 + if total(strsplit((*!dustem_scope).(i),'+',/extract) eq 'MODIFY') then begin
  27 +
  28 + ;making sure it's the stokes plugin that gets used here
  29 +
  30 + iidx=where(tag_names(*!dustem_scope) eq 'STEXT')
  31 +
  32 +
  33 + IF iidx eq i THEN BEGIN ;not sure about this condition but you get the spirit
  34 +
  35 + Q_ext=(*(*!dustem_plugin).(i))(0:n_elements((*(*!dustem_plugin).(i)))/2-1)
  36 + U_ext=(*(*!dustem_plugin).(i))(n_elements((*(*!dustem_plugin).(i)))/2:*)
  37 + dustem_qext=interpol(Q_ext,sti.polext.wav,(*!dustem_data.polext).wav)
  38 + dustem_uext=interpol(U_ext,sti.polext.wav,(*!dustem_data.polext).wav)
  39 +
  40 + endif
  41 + endif
  42 + endfor
  43 + ;endif
  44 + ;-------------------------------------------
  45 +
  46 +dustem_polext = interpol(sti.polext.ext_tot,sti.polext.wav,(*!dustem_data.polext).wav)
  47 +
  48 +
  49 +return, dustem_polext
  50 +
  51 +
  52 +END
0 53 \ No newline at end of file
... ...
src/idl/dustem_compute_polsed.pro 100644 → 100755
  Diff comments   View file @e7938fa
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
  Diff comments   View file @e7938fa
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
  Diff comments   View file @e7938fa
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
  Diff comments   View file @e7938fa
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_stext.pro 0 → 100644
  Diff comments   View file @e7938fa
... ... @@ -0,0 +1,32 @@
  1 +FUNCTION dustem_create_stext ,st, key=key, val=val, scope
  2 +
  3 +psi_ref_ext=0
  4 +
  5 +
  6 +out=['']
  7 +
  8 +IF keyword_set(key) THEN BEGIN
  9 + a=where(key EQ 1,count1)
  10 +
  11 +
  12 + IF count1 NE 0 then begin
  13 + psi_ref_ext=(val(a))(0)
  14 + ;test to not devide by zero
  15 + frac=st.polext.ext_tot/st.ext.ext_tot;st.polext.ext_tot*0+1
  16 + tes=where(finite(frac) eq 0)
  17 + frac(tes)=0.
  18 +
  19 + polar_ippsi2iqu,((st.polext).ext_tot),Q_ext,U_ext,frac,psi_ref_ext
  20 +
  21 + out=[out,Q_ext,U_ext]
  22 + endif
  23 +ENDIF
  24 +
  25 +scope=((*!dustem_scope).stext)
  26 +
  27 +return, out[1:*]
  28 +
  29 +end
  30 +
  31 +
  32 +END
0 33 \ No newline at end of file
... ...
src/idl/dustem_create_stokes.pro 0 → 100644
  Diff comments   View file @e7938fa
... ... @@ -0,0 +1,55 @@
  1 +Function dustem_create_stokes, st, key=key, val=val, scope
  2 +
  3 +psi_ref_dust=0
  4 +psi_ref_synch=0
  5 +
  6 +out=['']
  7 +
  8 +IF keyword_set(key) THEN BEGIN
  9 + a=where(key EQ 1,count1)
  10 + b=where(key EQ 2,count2)
  11 + c=where(key EQ 3,count3)
  12 +
  13 +
  14 + IF count3 NE 0 then begin
  15 +
  16 + frac_synch=(val(c))(0)
  17 +
  18 + endif else frac_synch=0.7 ;default value for synchrotron polarization fraction
  19 +
  20 +
  21 + IF count1 NE 0 then begin
  22 + psi_ref_dust=(val(a))(0)
  23 + ;test to not devide by zero
  24 +
  25 + frac=st.polsed.em_tot/st.sed.em_tot
  26 + tes=where(finite(frac) eq 0)
  27 + frac(tes)=0.
  28 +
  29 +
  30 + polar_ippsi2iqu,((st.sed).em_tot),Q_dust,U_dust,frac,psi_ref_dust
  31 +
  32 + out=[out,Q_dust,U_dust]
  33 + endif
  34 +
  35 + IF count2 NE 0 then begin
  36 +
  37 + psi_ref_synch=(val(b))(0)
  38 +
  39 +
  40 + frac=(*(*!dustem_plugin).synchrotron)/(*(*!dustem_plugin).synchrotron)
  41 + tes=where(finite(frac) eq 0)
  42 + frac(tes)=0.
  43 + frac*=frac_synch
  44 + polar_ippsi2iqu,(*(*!dustem_plugin).synchrotron),Q_synch,U_synch,frac,psi_ref_synch
  45 + out=[out,Q_synch,U_synch]
  46 + endif
  47 +
  48 +ENDIF
  49 +
  50 +scope=((*!dustem_scope).stokes)
  51 +
  52 +return, out[1:*]
  53 +
  54 +end
  55 +
... ...
src/idl/dustem_create_synchrotron.pro 100644 → 100755
  Diff comments   View file @e7938fa
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_fit_sed_extsed_readme.pro 100644 → 100755
  Diff comments   View file @e7938fa
... ... @@ -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
  Diff comments   View file @e7938fa
... ... @@ -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
  Diff comments   View file @e7938fa
... ... @@ -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_init.pro 100644 → 100755
  Diff comments   View file @e7938fa
... ... @@ -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'
... ... @@ -325,17 +345,17 @@ ENDIF
325 345 ;Initialize filters
326 346 dustem_filter_init,dir=dir
327 347  
328   -;=== test plugins
329   -plugins_desc_file=getenv('DUSTEM_SOFT_DIR')+'plugins_description.xcat'
330   -plugins_st=read_xcat(plugins_desc_file,/silent)
331   -Nplugins=n_elements(plugins_st)
332   -FOR i=0L,Nplugins-1 DO BEGIN
333   - str='res='+'dustem_'+plugins_st[i].PLUGIN_NAME+'()'
334   - IF !dustem_verbose NE 0 THEN BEGIN
335   - message,'Executing:'+str,/continue
336   - ENDIF
337   - toto=execute(str)
338   -ENDFOR
  348 +; ;=== test plugins; This part isn't needed. It was here before the implemetation of the plugins
  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
  350 +; plugins_st=read_xcat(plugins_desc_file,/silent)
  351 +; Nplugins=n_elements(plugins_st)
  352 +; FOR i=0L,Nplugins-1 DO BEGIN
  353 +; str='res='+'dustem_'+plugins_st[i].PLUGIN_NAME+'()'
  354 +; IF !dustem_verbose NE 0 THEN BEGIN
  355 +; message,'Executing:'+str,/continue
  356 +; ENDIF
  357 +; toto=execute(str)
  358 +; ENDFOR
339 359 ;stop
340 360  
341 361 the_end:
... ...
src/idl/dustem_init_fixed_params.pro 100644 → 100755
  Diff comments   View file @e7938fa
... ... @@ -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_plugins.pro 0 → 100644
  Diff comments   View file @e7938fa
... ... @@ -0,0 +1,97 @@
  1 +PRO dustem_init_plugins, pd
  2 +
  3 +
  4 +;DEFINING THE SCOPES--------------------
  5 +
  6 +varvar=create_struct('varvar',ptr_new())
  7 +
  8 +
  9 +for i=0L,n_elements(pd)-1 do begin
  10 +
  11 +fi=strtrim(strmid(pd(i),0,6),2)
  12 +
  13 +test=(n_tags(varvar) eq 1 && tag_exist(varvar,'varvar'))
  14 +
  15 +if fi eq 'dustem' then begin
  16 +
  17 +fii=strtrim(strmid(pd(i),14,4),2)
  18 +
  19 +j=0L
  20 +if i eq 0 then j+=1
  21 +
  22 +trl=fii eq strtrim(strmid(pd(i+j-1),14,4),2)
  23 +
  24 +
  25 +case fii of ;TESTING ON THE 4 FIRST LETTERS OF THE PLUGIN'S NAME
  26 +
  27 +;starting a fresh structure or appending to a pre-existing one
  28 +'cont': begin
  29 +
  30 +plugin_name='CONTINUUM' ; THIS STRING SHOULD MATCH THE NAME OF THE PLUGIN IN THE SCOPE LINE OF THE PLUGIN FUNCTION
  31 +if test && (~trl) then varvar=create_struct(plugin_name,'ADD_SED') else begin
  32 +if (~trl) then varvar=create_struct(varvar,plugin_name,'ADD_SED')
  33 +endelse
  34 +end
  35 +
  36 +'sync': begin
  37 +
  38 +plugin_name='SYNCHROTRON'
  39 +if test && (~trl) then varvar=create_struct(plugin_name,'ADD_SED+ADD_POLSED') else begin
  40 +if (~trl) then varvar=create_struct(varvar,plugin_name,'ADD_SED+ADD_POLSED')
  41 +endelse
  42 +end
  43 +
  44 +'free': begin
  45 +
  46 +plugin_name='FREE_FREE'
  47 +if test && (~trl) then varvar=create_struct(plugin_name,'ADD_SED') else begin
  48 +if (~trl) then varvar=create_struct(varvar,plugin_name,'ADD_SED')
  49 +endelse
  50 +end
  51 +
  52 +'stok': begin
  53 +
  54 +plugin_name='STOKES'
  55 +if test && (~trl) then varvar=create_struct(plugin_name,'MODIFY') else begin ; plugin that allows the computation of Q and U and fitting of PSI_dust and PSI_synch
  56 +if (~trl) then varvar=create_struct(varvar,plugin_name,'MODIFY')
  57 +endelse
  58 +end
  59 +
  60 +'stex':begin
  61 +
  62 +print,'found'
  63 +plugin_name='STEXT' ;'stokes'+extinction
  64 +if test && (~trl) then varvar=create_struct(plugin_name,'MODIFY') else begin ; plugin that allows the computation of Q_ext = sigext_q and U_ext=sigext_u and fitting of PSI_ext
  65 +if (~trl) then varvar=create_struct(varvar,plugin_name,'MODIFY')
  66 +endelse
  67 +print, varvar
  68 +end
  69 +
  70 +;NEW PLUGIN TO BE ADDED HERE
  71 +endcase
  72 +
  73 +endif
  74 +
  75 +endfor
  76 +
  77 +
  78 +test=(n_tags(varvar) eq 1 && tag_exist(varvar,'varvar'))
  79 +if test then varvar=create_struct('varvar','NULL')
  80 +
  81 +defsysv, '!dustem_scope', ptr_new(varvar)
  82 +;---------------------------------------------
  83 +
  84 +
  85 +;INITIALIZING THE PLUGIN DATA ARRAYS----------
  86 +tgnms=tag_names(varvar)
  87 +
  88 +for i=0L, n_tags(varvar)-1 do begin
  89 +
  90 +if i eq 0 then strct= create_struct(tgnms(i),ptr_new())
  91 +if i ge 1 then strct= create_struct(strct, tgnms(i),ptr_new())
  92 +endfor
  93 +
  94 +defsysv, '!dustem_plugin', ptr_new(strct) ; structure that will contain the spectra associated to each plugin
  95 +;---------------------------------------------
  96 +
  97 +end
0 98 \ No newline at end of file
... ...
src/idl/dustem_mpfit_run.pro 100644 → 100755
  Diff comments   View file @e7938fa
... ... @@ -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_plot_ext.pro 100644 → 100755
  Diff comments   View file @e7938fa
... ... @@ -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_fit_sed.pro 100644 → 100755
  Diff comments   View file @e7938fa
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
  Diff comments   View file @e7938fa
... ... @@ -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_polar.pro 100644 → 100755
  Diff comments   View file @e7938fa
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_read_align.pro 100644 → 100755
  Diff comments   View file @e7938fa
src/idl/dustem_read_dustem_lv.pro 100644 → 100755
  Diff comments   View file @e7938fa
src/idl/dustem_read_ext_lv.pro 100644 → 100755
  Diff comments   View file @e7938fa
src/idl/dustem_sed_plot.pro 100644 → 100755
  Diff comments   View file @e7938fa
... ... @@ -5,9 +5,13 @@ IF not keyword_set(function_name) THEN BEGIN
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,freefree=freefree,synchrotron=synchrotron,out_st=out_st)
  8 +
  9 + ;dustem_sed=dustem_compute_sed(p_min,st=st,cont=cont,freefree=freefree,synchrotron=synchrotron,out_st=out_st)
  10 + dustem_sed=dustem_compute_sed(p_min,st=st,out_st=out_st)
  11 +
9 12 st=out_st
10   - dustem_plot_fit_sed,st,dustem_sed,cont,freefree,synchrotron,_extra=_extra,legend_xpos=legend_xpos,legend_ypos=legend_ypos
  13 + ;dustem_plot_fit_sed,st,dustem_sed,cont,freefree,synchrotron,_extra=_extra,legend_xpos=legend_xpos,legend_ypos=legend_ypos
  14 + dustem_plot_fit_sed,st,dustem_sed,_extra=_extra,legend_xpos=legend_xpos,legend_ypos=legend_ypos
11 15 IF keyword_set(pol) THEN BEGIN
12 16 stop
13 17 ;dustem_polsed=dustem_compute_polsed(p_min,st=pst,cont=cont,freefree=freefree,synchrotron=synchrotron,out_st=out_st)
... ...
src/idl/dustem_set_data.pro 100644 → 100755
  Diff comments   View file @e7938fa
1   -FUNCTION dustem_set_data,help=help,ext=ext,polsed=polsed,polext=polext,sed=sed,polfrac=polfrac,rchi2_weight=rchi2_weight,f_HI=f_HI
  1 +FUNCTION dustem_set_data,help=help,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
2 2  
3 3 ;+
4 4 ; NAME:
... ... @@ -56,6 +56,12 @@ ENDIF
56 56 if keyword_set(polext) then !dustem_data.polext = ptr_new()
57 57 if keyword_set(polsed) then !dustem_data.polsed = ptr_new()
58 58 if keyword_set(polfrac) then !dustem_data.polfrac = ptr_new()
  59 +if keyword_set(qsed) then !dustem_data.qsed = ptr_new()
  60 +if keyword_set(used) then !dustem_data.used = ptr_new()
  61 +if keyword_set(qext) then !dustem_data.qext = ptr_new()
  62 +if keyword_set(uext) then !dustem_data.uext = ptr_new()
  63 +
  64 +
59 65  
60 66 IF keyword_set(sed) THEN BEGIN
61 67 wavs=sed.wave
... ... @@ -95,6 +101,84 @@ IF keyword_set(polsed) THEN BEGIN
95 101 ;VGe
96 102 ENDIF
97 103  
  104 +
  105 +IF keyword_set(qsed) THEN BEGIN
  106 + wavs=qsed.wave
  107 + ;=== Impose central wavelengths for photometric channels
  108 + ind=where(qsed.filter NE 'SPECTRUM',count)
  109 + IF count NE 0 THEN BEGIN
  110 + wavs(ind)=dustem_filter2wav(qsed(ind).filter)
  111 + ENDIF
  112 + ;=== define observations structure
  113 + obs_st={instru_names:qsed.instru,filt_names:qsed.filter,wav:wavs, values:qsed.spec,sigma:qsed.error}
  114 + ;=== fill !dustem_data
  115 + !dustem_data.qsed = ptr_new(obs_st)
  116 + ;VGe
  117 + ; If f_HI is specified, multiply the data by f_HI
  118 + if keyword_set(f_HI) then if f_HI gt 0 then $
  119 + (*!dustem_data.qsed).values *= f_HI
  120 + ;VGe
  121 +ENDIF
  122 +
  123 +IF keyword_set(used) THEN BEGIN
  124 + wavs=used.wave
  125 + ;=== Impose central wavelengths for photometric channels
  126 + ind=where(used.filter NE 'SPECTRUM',count)
  127 + IF count NE 0 THEN BEGIN
  128 + wavs(ind)=dustem_filter2wav(used(ind).filter)
  129 + ENDIF
  130 + ;=== define observations structure
  131 + obs_st={instru_names:used.instru,filt_names:used.filter,wav:wavs, values:used.spec,sigma:used.error}
  132 + ;=== fill !dustem_data
  133 + !dustem_data.used = ptr_new(obs_st)
  134 + ;VGe
  135 + ; If f_HI is specified, multiply the data by f_HI
  136 + if keyword_set(f_HI) then if f_HI gt 0 then $
  137 + (*!dustem_data.used).values *= f_HI
  138 + ;VGe
  139 +ENDIF
  140 +
  141 +
  142 +
  143 +IF keyword_set(qext) THEN BEGIN
  144 + wavs=qext.wave
  145 + ;=== Impose central wavelengths for photometric channels
  146 + ind=where(qext.filter NE 'SPECTRUM',count)
  147 + IF count NE 0 THEN BEGIN
  148 + wavs(ind)=dustem_filter2wav(qext(ind).filter)
  149 + ENDIF
  150 + ;=== define observations structure
  151 + obs_st={instru_names:qext.instru,filt_names:qext.filter,wav:wavs, values:qext.spec,sigma:qext.error}
  152 + ;=== fill !dustem_data
  153 + !dustem_data.qext = ptr_new(obs_st)
  154 + ;VGe
  155 + ; If f_HI is specified, multiply the data by f_HI
  156 + if keyword_set(f_HI) then if f_HI gt 0 then $
  157 + (*!dustem_data.qext).values *= f_HI
  158 + ;VGe
  159 +ENDIF
  160 +
  161 +IF keyword_set(uext) THEN BEGIN
  162 + wavs=uext.wave
  163 + ;=== Impose central wavelengths for photometric channels
  164 + ind=where(uext.filter NE 'SPECTRUM',count)
  165 + IF count NE 0 THEN BEGIN
  166 + wavs(ind)=dustem_filter2wav(uext(ind).filter)
  167 + ENDIF
  168 + ;=== define observations structure
  169 + obs_st={instru_names:uext.instru,filt_names:uext.filter,wav:wavs, values:uext.spec,sigma:uext.error}
  170 + ;=== fill !dustem_data
  171 + !dustem_data.uext = ptr_new(obs_st)
  172 + ;VGe
  173 + ; If f_HI is specified, multiply the data by f_HI
  174 + if keyword_set(f_HI) then if f_HI gt 0 then $
  175 + (*!dustem_data.uext).values *= f_HI
  176 + ;VGe
  177 +ENDIF
  178 +
  179 +
  180 +
  181 +
98 182 IF keyword_set(polfrac) THEN BEGIN
99 183 wavs=polfrac.wave
100 184 ;=== Impose central wavelengths for photometric channels
... ... @@ -135,11 +219,19 @@ IF keyword_set(ext) THEN BEGIN
135 219 !dustem_data.ext = ptr_new(obs_st)
136 220 ENDIF
137 221  
  222 +
  223 +
  224 +
  225 +
138 226 !fit_rchi2_weight.sed=1.
139 227 !fit_rchi2_weight.ext=1.
140 228 if keyword_set(polext) then !fit_rchi2_weight.polext=1.
141 229 if keyword_set(polsed) then !fit_rchi2_weight.polsed=1.
142 230 if keyword_set(polfrac) then !fit_rchi2_weight.polfrac=1.
  231 +if keyword_set(qsed) then !fit_rchi2_weight.qsed=1.
  232 +if keyword_set(used) then !fit_rchi2_weight.used=1.
  233 +if keyword_set(qext) then !fit_rchi2_weight.qext=1.
  234 +if keyword_set(uext) then !fit_rchi2_weight.uext=1.
143 235  
144 236 IF keyword_set(rchi2_weight) THEN begin
145 237 !fit_rchi2_weight.sed = rchi2_weight.sed
... ... @@ -147,6 +239,10 @@ IF keyword_set(rchi2_weight) THEN begin
147 239 if keyword_set(polext) then !fit_rchi2_weight.polext = rchi2_weight.polext
148 240 if keyword_set(polsed) then !fit_rchi2_weight.polsed = rchi2_weight.polsed
149 241 if keyword_set(polfrac) then !fit_rchi2_weight.polfrac = rchi2_weight.polfrac
  242 + if keyword_set(qsed) then !fit_rchi2_weight.qsed = rchi2_weight.qsed
  243 + if keyword_set(used) then !fit_rchi2_weight.used = rchi2_weight.used
  244 + if keyword_set(qext) then !fit_rchi2_weight.qext = rchi2_weight.qext
  245 + if keyword_set(uext) then !fit_rchi2_weight.uext = rchi2_weight.uext
150 246 endif
151 247  
152 248 return,ptr_new(obs_st)
... ...
src/idl/dustem_set_params.pro 100644 → 100755
  Diff comments   View file @e7938fa
... ... @@ -30,8 +30,18 @@ FOR i=0L,Nparams-1 DO BEGIN
30 30 'PLUGIN': BEGIN
31 31 length=strlen((*(*!dustem_fit).param_descs)(i))
32 32 strr=strmid((*(*!dustem_fit).param_descs)(i),14,length-16) ;This is not very clean programming ....
  33 +
  34 + ;strr= strupcase(strmid(strmid((*(*!dustem_fit).param_descs)(i),0,strlen((*(*!dustem_fit).param_descs)(i))-2),14))
  35 +
  36 + k=where(tag_names(*!dustem_scope) eq strupcase(strr), countt)
  37 +
  38 +
33 39 ; stop
34   - IF strr NE 'continuum' AND strr NE 'freefree' AND strr NE 'synchrotron' THEN BEGIN ;This is not very clean programming ....
  40 +
  41 + ;IF strr NE 'continuum' AND strr NE 'freefree' AND strr NE 'synchrotron' AND strr NE 'powerlaw' THEN BEGIN ;This is not very clean programming ....
  42 +
  43 + IF countt EQ 0 THEN BEGIN
  44 +
35 45 IF strr eq 'isrf2' then strr = 'isrf'
36 46 IF !dustem_which EQ 'VERSTRAETE' THEN BEGIN
37 47 ddir=!dustem_dat+'les_DAT/'
... ...
src/idl/dustem_write_all.pro 100644 → 100755
  Diff comments   View file @e7938fa
... ... @@ -52,7 +52,7 @@ IF !dustem_which EQ &#39;VERSTRAETE&#39; THEN BEGIN
52 52 ENDIF
53 53  
54 54 ;IF getenv('DUSTEM_WHICH') EQ 'WEB3p8' THEN BEGIN
55   -IF !dustem_which EQ 'WEB3p8' THEN BEGIN
  55 +IF !dustem_which EQ 'WEB3p8' THEN BEGIN ; changed this as well
56 56 dustem_write_all_web3p8,st,dir_out
57 57 goto,sortie
58 58 ENDIF
... ...