OV-GSO-DC / dustem-wrapper_idl

Commit 83b3ddeecfddf348ed9cfa136da997e80447bb89

Authored by Jean-Philippe Bernard
1 parent e22e01bf
Exists in master

modified with Ilyes for new polarization fitting startegy

Showing 5 changed files with 148 additions and 124 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_compute_polsed.pro
  Diff comments   View file @83b3dde
1   -FUNCTION dustem_compute_polsed,p_dim,stp,_extra=extra,out_st=out_st,dustem_qsed,dustem_used,Q_sed,U_sed
  1 +FUNCTION dustem_compute_polsed,p_dim,stp,dustem_qsed,dustem_used,Q_sed,U_sed,_extra=extra,out_st=out_st
2 2  
3 3 ;+
4 4 ; NAME:
... ... @@ -12,11 +12,13 @@ FUNCTION dustem_compute_polsed,p_dim,stp,_extra=extra,out_st=out_st,dustem_qsed,
12 12 ; INPUTS:
13 13 ; p_dim = parameter values
14 14 ; OPTIONAL INPUT PARAMETERS:
15   -; None
  15 +; stp = Dustem structure
16 16 ; OUTPUTS:
17   -; sed = computed SED for filters in !dustem_data
  17 +; sed : computed SED for filters in !dustem_data
  18 +; dustem_qsed : SED in Q
  19 +; dustem_used : SED in U
18 20 ; OPTIONAL OUTPUT PARAMETERS:
19   -; st = Dustem output structure
  21 +; stp = Dustem output structure
20 22 ; ACCEPTED KEY-WORDS:
21 23 ; help = If set, print this help
22 24 ; COMMON BLOCKS:
... ... @@ -33,7 +35,6 @@ FUNCTION dustem_compute_polsed,p_dim,stp,_extra=extra,out_st=out_st,dustem_qsed,
33 35 ; Written by J.-Ph. Bernard
34 36 ; see evolution details on the dustem cvs maintained at CESR
35 37 ; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
36   -;-
37 38  
38 39 IF keyword_set(help) THEN BEGIN
39 40 doc_library,'dustem_compute_polsed'
... ... @@ -43,65 +44,95 @@ ENDIF
43 44  
44 45 IF not keyword_set(stp) THEN BEGIN
45 46 stp=dustem_run(p_dim)
46   - dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),res=res,stp
  47 + dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),stp,res=res
47 48 ENDIF
48 49  
49 50 If keyword_set(result) THEN result=stp
50 51 fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/(stp.polsed).wav)*1.e20/1.e7
51   -spec = stp.polsed.em_tot * fact
  52 +StokesI = stp.sed.em_tot * fact ;This is Total intensity I
  53 +largeP = stp.polsed.em_tot * fact ;This is Polarized intensity P
  54 +
  55 +Nwaves=(size(StokesI))[1]
  56 +polar_ippsi2iqu,StokesI,StokesQ,StokesU,largeP/StokesI,replicate(0.,Nwaves)
  57 +
52 58 ;stop
53 59  
  60 +;MUST INCLUDE HERE ACTION OF PLUGIN WITH SCOPE MODIFY_POLAR_ANGLE
  61 +FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
  62 + IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'MODIFY_POLAR_ANGLE') THEN BEGIN
  63 + angle=0. ;A TROUVER !!!
  64 + StokesQ=largeP*cos(2.*angle)
  65 + StokesU=largeP*sin(2.*angle)
  66 + ENDIF
  67 +ENDFOR
  68 +
54 69 ;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 spec+=(*(*!dustem_plugin).(i))[*,0]
57   -endfor
  70 +FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
  71 + IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'ADD_POLSED') THEN BEGIN
  72 + StokesI+=(*(*!dustem_plugin).(i))[*,0]
  73 + StokesQ+=(*(*!dustem_plugin).(i))[*,1]
  74 + StokesU+=(*(*!dustem_plugin).(i))[*,2]
  75 + ENDIF
  76 +ENDFOR
58 77 ;-----------------------------------------
59 78  
  79 +largeP=sqrt(StokesQ^2+StokesU^2)
60 80  
61   -if not isarray(spec) THEN stop
62   -;stop
63   -;COMPUTE THE MODEL SED
  81 +dustem_qsed = (*!dustem_data.polsed).values * 0.
  82 +dustem_used = dustem_qsed
64 83  
65 84 dustem_polsed = (*!dustem_data.polsed).values * 0.
66 85  
  86 +goto,on_saute
  87 +
  88 +Q_sed=StokesQ
  89 +U_sed=StokesU
  90 +
  91 +if not isarray(StokesI) THEN stop
  92 +;stop
  93 +;COMPUTE THE MODEL SED
  94 +
  95 +;This will be the output of the fucntion
67 96  
68 97 ;MODIFYING VIA PLUGIN-------------- THIS BLOCK NEEDS TO BE UNDER THE DUSTEM_POLSED LINE
69   -for i=0L,n_tags(*!dustem_scope)-1 do begin
  98 +FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
70 99  
71   -if total(strsplit((*(*!dustem_scope).(i)),'+',/extract) eq 'MODIFY_POLSED') then begin
  100 + IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'MODIFY_POLSED') THEN BEGIN
72 101  
73   -;making sure it's the stokes plugin that gets used here
  102 + ;making sure it's the stokes plugin that gets used here
74 103  
75   -iidx=where(tag_names(*!dustem_scope) eq 'SYNCHROTRON',countiidx)
  104 + iidx=where(tag_names(*!dustem_scope) eq 'SYNCHROTRON',countiidx)
76 105  
77   -IF countiidx ne 0 THEN BEGIN
  106 + IF countiidx ne 0 THEN BEGIN
78 107  
79   -Q_dust=(*(*!dustem_plugin).(i))(0:n_elements((*(*!dustem_plugin).(i)))/4-1)
80   -U_dust=(*(*!dustem_plugin).(i))(n_elements((*(*!dustem_plugin).(i)))/4:n_elements((*(*!dustem_plugin).(i)))/2-1)
81   -Q_synch=(*(*!dustem_plugin).(i))(n_elements((*(*!dustem_plugin).(i)))/2:3*n_elements((*(*!dustem_plugin).(i)))/4-1)
82   -U_synch=(*(*!dustem_plugin).(i))(3*n_elements((*(*!dustem_plugin).(i)))/4:*)
  108 + Q_dust=(*(*!dustem_plugin).(i))(0:n_elements((*(*!dustem_plugin).(i)))/4-1)
  109 + U_dust=(*(*!dustem_plugin).(i))(n_elements((*(*!dustem_plugin).(i)))/4:n_elements((*(*!dustem_plugin).(i)))/2-1)
  110 + Q_synch=(*(*!dustem_plugin).(i))(n_elements((*(*!dustem_plugin).(i)))/2:3*n_elements((*(*!dustem_plugin).(i)))/4-1)
  111 + U_synch=(*(*!dustem_plugin).(i))(3*n_elements((*(*!dustem_plugin).(i)))/4:*)
83 112  
84   -ENDIF ELSE BEGIN
85   -Q_dust=(*(*!dustem_plugin).(i))(0:n_elements((*(*!dustem_plugin).(i)))/2-1)
86   -U_dust=(*(*!dustem_plugin).(i))(n_elements((*(*!dustem_plugin).(i)))/2:*)
87   -Q_synch=Q_dust*0
88   -U_synch=U_dust*0
89   -ENDELSE
  113 + ENDIF ELSE BEGIN
  114 + Q_dust=(*(*!dustem_plugin).(i))(0:n_elements((*(*!dustem_plugin).(i)))/2-1)
  115 + U_dust=(*(*!dustem_plugin).(i))(n_elements((*(*!dustem_plugin).(i)))/2:*)
  116 + Q_synch=Q_dust*0
  117 + U_synch=U_dust*0
  118 + ENDELSE
90 119  
91   -;THE ABOVE VARIABLES SHOULD BE ACCESSED DIRECTLY BY THE MPFIT_RUN.PRO
92   -Q_sed=Q_dust*fact+Q_synch
93   -U_sed=U_dust*fact+U_synch
  120 + ;THE ABOVE VARIABLES SHOULD BE ACCESSED DIRECTLY BY THE MPFIT_RUN.PRO
  121 + Q_sed=Q_dust*fact+Q_synch
  122 + U_sed=U_dust*fact+U_synch
94 123  
95   -;data the dustem_wrap will use
96   -dustem_qsed = dustem_polsed
97   -dustem_used = dustem_polsed
  124 + ;data the dustem_wrap will use
  125 + dustem_qsed = dustem_polsed
  126 + dustem_used = dustem_polsed
98 127  
99   -endif
100   -endfor
101   -;-------------------------------------------
  128 + ENDIF
  129 +ENDFOR
102 130  
  131 +on_saute:
103 132  
104   -;stop
  133 +;-------------------------------------------
  134 +
  135 + ;stop
105 136 IF !dustem_do_cc NE 0 AND !dustem_never_do_cc EQ 0 THEN BEGIN
106 137 message,'DOING color correction calculations',/info
107 138 ENDIF ELSE BEGIN
... ... @@ -109,20 +140,22 @@ ENDIF ELSE BEGIN
109 140 ENDELSE
110 141  
111 142 ind_polsed=where((*!dustem_data.polsed).filt_names NE 'SPECTRUM',count_polsed)
112   -
  143 +
113 144 IF count_polsed NE 0 THEN BEGIN ; INCLUDE COLOR CORRECTIONS OF (Q,U) for dust & synch
114 145 filter_names=((*!dustem_data.polsed).filt_names)(ind_polsed)
115 146  
116   - spolsed=dustem_cc(stp.sed.wav,spec,filter_names,cc=cc) & dustem_polsed(ind_polsed)=spolsed
  147 + spolsed=dustem_cc(stp.sed.wav,largeP,filter_names,cc=cc)
  148 + dustem_polsed[ind_polsed]=spolsed
117 149  
118   -
119 150 ;PLUGIN 'STOKES' W/ scope='MODIFY' PART------------------------------------
120 151  
121 152 if isa(dustem_qsed) then begin
122   - sqsed=dustem_cc(stp.sed.wav,Q_sed,filter_names,cc=cc) & dustem_qsed(ind_polsed)=sqsed
  153 + sqsed=dustem_cc(stp.sed.wav,StokesQ,filter_names,cc=cc)
  154 + dustem_qsed(ind_polsed)=sqsed
123 155 endif
124 156 if isa(dustem_used) then begin
125   - sused=dustem_cc(stp.sed.wav,U_sed,filter_names,cc=cc) & dustem_used(ind_polsed)=sused
  157 + sused=dustem_cc(stp.sed.wav,StokesU,filter_names,cc=cc)
  158 + dustem_used(ind_polsed)=sused
126 159 endif
127 160 ;--------------------------------------------------------------------------
128 161  
... ... @@ -133,12 +166,10 @@ ENDIF
133 166 ;wavelengths points exist in the model (long wavelengths).
134 167 ind_spec=where((*!dustem_data.polsed).filt_names EQ 'SPECTRUM',count_spec)
135 168 IF count_spec NE 0 THEN BEGIN
136   - dustem_polsed(ind_spec)=interpol(spec,stp.wav,(((*!dustem_data.polsed).wav)(ind_spec)))
137   -
  169 + dustem_polsed(ind_spec)=interpol(largeP,stp.wav,(((*!dustem_data.polsed).wav)(ind_spec)))
138 170 ;STOKES PARAMS TREATMENT
139   - if isa(dustem_qsed) then dustem_qsed(ind_spec)=interpol(Q_sed,stp.wav,(((*!dustem_data.polsed).wav)(ind_spec)))
140   - if isa(dustem_used) then dustem_used(ind_spec)=interpol(U_sed,stp.wav,(((*!dustem_data.polsed).wav)(ind_spec)))
141   -
  171 + if isa(dustem_qsed) then dustem_qsed(ind_spec)=interpol(StokesQ,stp.wav,(((*!dustem_data.polsed).wav)(ind_spec)))
  172 + if isa(dustem_used) then dustem_used(ind_spec)=interpol(StokesU,stp.wav,(((*!dustem_data.polsed).wav)(ind_spec)))
142 173 ENDIF
143 174 out_st=stp
144 175  
... ...
src/idl/dustem_fit_sed_polsed_readme.pro
  Diff comments   View file @83b3dde
... ... @@ -73,11 +73,11 @@ polsed=sed ; comment if polext is to be used
73 73  
74 74 ;UNCOMMENT THIS BLOCK AFTER THE TESTING
75 75 ; ;stokes fake parameters - comment if stokes extinction parameters are to be used
76   -; qsed=sed
77   -; used=sed
78   -
79   -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)
  76 +qsed=sed
  77 +used=sed
80 78  
  79 +;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)
  80 +st=dustem_set_data(ext=ext,polsed=polsed,polext=polext,sed=sed,polfrac=polfrac,rchi2_weight=rchi2_weight,f_HI=f_HI,qext=qext,uext=uext)
81 81  
82 82 ;=== Set which parameters you want to fit
83 83 pd = [ $
... ... @@ -85,37 +85,24 @@ pd = [ $
85 85 '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
86 86 '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
87 87 '(*!dustem_params).grains(2).mdust_o_mh',$ ;amCBEx
88   - 'dustem_plugin_synchrotron_2' $ ;Synchrotron amplitude
89   -; 'dustem_create_stokes_1',$ ;Polarization angle psi_ref_dust
  88 + 'dustem_plugin_synchrotron_2', $ ;Synchrotron amplitude
  89 + 'dustem_plugin_synchrotron_4', $ ;Synchrotron angle
  90 + 'dustem_create_stokes_1' $ ;Polarization angle appliued to the Fortran dust model
90 91 ; 'dustem_create_stokes_2',$;Polarization angle psi_ref_synch
91 92 ; 'dustem_create_stokes_3'$;Polarization fraction of synchrotron spectrum
92 93 ;'dustem_create_stext_1'$;Polarization angle in extinction PSI_ext
93 94 ]
94 95  
95   -p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04,5e-3]
96   -; ;p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04,25]
97   -
98   -
99   -;UNCOMMENT THIS BLOCK AFTER THE TESTING
100   -;p_truth=[1.,7.8000E-04,7.8000E-04,2E-03,5e-3,70,20,0.2]
101   -
102   -;TESTING THE ISRF PLUGIN
103   -
104   -;p_truth=[1.,7.8000E-04,7.8000E-04,2E-03,2e-3,13.,3.,9.,2.]
105   -;p_truth=[7.8000E-04,7.8000E-04,2E-03,2e-3,13.,3.];,9.,2.]
106   -;p_truth=[7.8000E-04,7.8000E-04,2E-03,10.,8.]
107   -;p_truth=[10.,8.]
108   -
  96 +p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04,5e-3,0.,0.]
109 97  
110 98 ;=== set initial value of parameters you want to fit
111   -iv=p_truth ;exact solution
  99 +;iv=p_truth ;exact solution
112 100  
113 101  
114 102  
115 103 ;UNCOMMENT THIS BLOCK AFTER THE TESTING
116 104  
117   - iv=p_truth+[0.,3.e-5,1.e-5,-2.e-5,5e-2,5,5,0.2];shifted from solution
118   -; ;iv=p_truth+[0.,3.e-5,1.e-5,-2.e-5,5];shifted from solution
  105 + iv=p_truth+[0.,3.e-5,1.e-5,-2.e-5,5e-4,5.,-5.] ;shifted from solution
119 106  
120 107  
121 108 ;TESTING THE ISRF PLUGIN
... ... @@ -128,26 +115,24 @@ llims=replicate(0,Npar)
128 115 ;== SET THE FITTED PARAMETERS (now done in dustem_init_parinfo)
129 116 dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
130 117  
131   -;dustem_init_fixed_params,fpd,fiv
132   -
133   -;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-chrg-spin'
134   -;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-chrg-spin'
135   -;(*!dustem_params).grains[4].TYPE_KEYWORDS='pol-plaw-ed'
136   -
137   -
138   -dustem_init_plugins, pd
  118 +dustem_init_plugins,pd
139 119  
140 120 help,*!dustem_fit
141 121  
142 122 help,!run_pol
143 123  
  124 +;stop
  125 +toto=dustem_compute_sed(p_truth,sst,_extra=extra)
  126 +
  127 +titi=dustem_compute_polsed(p_truth,sstp,dustem_qsed,dustem_used,_extra=extra)
  128 +
144 129 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
145 130 sed.error=(sed.spec*0.1);>0.1
146 131  
147 132 ;UNCOMMENT THIS BLOCK AFTER THE TESTING
148   -;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
  133 +polsed.spec=dustem_compute_polsed(p_truth,sstp,dustem_qsed,dustem_used,_extra=extra) ;sst is not set on purpose, for dustem_compute_sed to compute the sed using fortran code
149 134  
150   -polsed.spec=dustem_compute_polsed(p_truth,sstp,_extra=extra)
  135 +;polsed.spec=dustem_compute_polsed(p_truth,sstp,_extra=extra)
151 136  
152 137 polsed.error=(polsed.spec*0.1)
153 138  
... ... @@ -162,12 +147,12 @@ polsed.error=(polsed.spec*0.1)
162 147 ;UNCOMMENT THIS BLOCK AFTER THE TESTING
163 148  
164 149 ; ;SET qsed and used data - comment if extinction polarization data is to be used
165   -; qsed.spec=dustem_qsed
166   -; qsed.error=qsed.spec*0.1
  150 +qsed.spec=dustem_qsed
  151 +qsed.error=qsed.spec*0.1
167 152 ;
168 153 ;
169   -; used.spec=dustem_used
170   -; used.error=used.spec*0.1
  154 +used.spec=dustem_used
  155 +used.error=used.spec*0.1
171 156  
172 157  
173 158 ; ;SET qext and uext data - comment if polarizarion data is to be used
... ... @@ -184,15 +169,10 @@ polsed.error=(polsed.spec*0.1)
184 169 ind=where(sed.wave GE 0.1)
185 170  
186 171 ;UNCOMMENT THIS BLOCK AFTER THE TESTING
187   -; st=dustem_set_data(polsed=polsed[ind],sed=sed[ind],qsed=qsed[ind],used=used[ind]) ; comment if second line is to be used
  172 +;st_bidon=dustem_set_data(polsed=polsed[ind],sed=sed[ind],qsed=qsed[ind],used=used[ind]) ; comment if second line is to be used
  173 +st_bidon=dustem_set_data(polsed=polsed[ind],sed=sed[ind]) ; comment if second line is to be used
188 174 ; ;st=dustem_set_data(sed=sed[ind],polext=polext[ind],qext=qext[ind],uext=uext[ind])
189   -
190   -
191   -st=dustem_set_data(polsed=polsed[ind],sed=sed[ind])
192   -
193   -
194   -
195   -
  175 +;st=dustem_set_data(polsed=polsed[ind],sed=sed[ind])
196 176  
197 177  
198 178 ;=== RUN fit
... ... @@ -222,7 +202,7 @@ file_out='/tmp/DUSTEM_polsed_fit_example.sav'
222 202 dustem_save_system_variables,file_out
223 203 message,'Saved '+file_out,/continue
224 204  
225   -stop
  205 +;stop
226 206  
227 207 ;======================================
228 208 ;====You can exit IDL here and re-enter - THIS PART APPEARS TO BE FAULTY. NEEDS INVESTIGATION
... ...
src/idl/dustem_init.pro
  Diff comments   View file @83b3dde
... ... @@ -72,6 +72,12 @@ defsysv, '!run_lin', 0.
72 72  
73 73 defsysv, '!dustem_redshift', 0. ;in case the SED needs to be shifted for a non zero redshift
74 74  
  75 +;This is just to make sure that plugin system variables are created even if no plugins are used
  76 +pd_bidon=['NONE']
  77 +dustem_init_plugins,pd_bidon
  78 +
  79 +;stop
  80 +
75 81 ;==== define a stellar modyfier to the ISRF
76 82  
77 83 ;toto=dustem_define_isrf_stellar_modifyer_variable()
... ...
src/idl/dustem_init_plugins.pro
  Diff comments   View file @83b3dde
... ... @@ -8,6 +8,7 @@ varvar=create_struct('varvar',ptr_new())
8 8 plugin_names=['']
9 9 plugind_detect_string='dustem_plugin'
10 10  
  11 +Nplugins=0L
11 12 for i=0L,n_elements(pd)-1 do begin
12 13 fi=strtrim(strmid(pd(i),0,strlen(plugind_detect_string)),2)
13 14 if fi eq plugind_detect_string then begin
... ... @@ -16,20 +17,26 @@ for i=0L,n_elements(pd)-1 do begin
16 17 ii = ii(countx-1)-1 ; Locating the last underscore to automate the extraction of the plugin's keyword
17 18 ftn = strmid(ftn,14,ii-14) ; String containing the name of the plugin without the associated keyword
18 19 plugin_names=[plugin_names,ftn]
  20 + Nplugins=Nplugins+1
19 21 ENDIF
20 22 ENDFOR
21 23  
22   -plugin_names=plugin_names[1:*]
23   -order=sort(plugin_names)
24   -plugin_names=plugin_names[order]
25   -un=uniq(plugin_names)
26   -plugin_names=plugin_names[un]
27   -Nplugins=n_elements(plugin_names)
28   -
29   -print, 'plugin_names is:'
  24 +IF Nplugins EQ 0 THEN BEGIN
  25 + plugin_names=['NONE']
  26 + Nplugins=1
  27 +ENDIF ELSE BEGIN
  28 + plugin_names=plugin_names[1:*]
  29 + order=sort(plugin_names)
  30 + plugin_names=plugin_names[order]
  31 + un=uniq(plugin_names)
  32 + plugin_names=plugin_names[un]
  33 + Nplugins=n_elements(plugin_names)
  34 +ENDELSE
  35 +
  36 +message, 'plugin_names is:',/continue
30 37 print, plugin_names
31 38  
32   -print, 'Nplugins is:'
  39 +message, 'Nplugins is:',/continue
33 40 print, Nplugins
34 41  
35 42  
... ... @@ -43,7 +50,7 @@ IF Nplugins NE 0 THEN BEGIN
43 50 str='varvar=create_struct('+str_fin+')'
44 51 toto=execute(str)
45 52  
46   - print, 'varvar is:'
  53 + message, 'varvar is:',/continue
47 54 print, varvar
48 55  
49 56 defsysv, '!dustem_scope', ptr_new(varvar)
... ... @@ -52,13 +59,13 @@ IF Nplugins NE 0 THEN BEGIN
52 59  
53 60 varvar=create_struct(plugin_names,ptr_new())
54 61 defsysv, '!dustem_scope', ptr_new(varvar)
55   - ENDELSE
  62 + ENDELSE
56 63  
57 64 ENDIF
58 65  
59 66  
60   -test=(n_tags(varvar) eq 1 && tag_exist(varvar,'varvar'))
61   -if test then varvar=create_struct('varvar','ERROR IN DUSTEM_INIT_PLUGINS!')
  67 +;test=(n_tags(varvar) eq 1 && tag_exist(varvar,'varvar'))
  68 +;if test then varvar=create_struct('varvar','ERROR IN DUSTEM_INIT_PLUGINS!')
62 69  
63 70 defsysv, '!dustem_scope', ptr_new(varvar)
64 71 ;---------------------------------------------
... ...
src/idl/dustem_set_data.pro
  Diff comments   View file @83b3dde
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
  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,qext=qext,uext=uext
  2 +; qsed=qsed,used=used
2 3  
3 4 ;+
4 5 ; NAME:
... ... @@ -56,8 +57,12 @@ ENDIF
56 57 if keyword_set(polext) then !dustem_data.polext = ptr_new()
57 58 if keyword_set(polsed) then !dustem_data.polsed = ptr_new()
58 59 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()
  60 +indQ=where(sed.StokesQ NE la_undef(),NQ)
  61 +if NQ NE 0 then !dustem_data.qsed = ptr_new()
  62 +;if keyword_set(qsed) then !dustem_data.qsed = ptr_new()
  63 +indU=where(sed.StokesU NE la_undef(),NU)
  64 +if NU NE 0 then !dustem_data.used = ptr_new()
  65 +;if keyword_set(used) then !dustem_data.used = ptr_new()
61 66 if keyword_set(qext) then !dustem_data.qext = ptr_new()
62 67 if keyword_set(uext) then !dustem_data.uext = ptr_new()
63 68  
... ... @@ -102,15 +107,17 @@ IF keyword_set(polsed) THEN BEGIN
102 107 ENDIF
103 108  
104 109  
105   -IF keyword_set(qsed) THEN BEGIN
  110 +;IF keyword_set(qsed) THEN BEGIN
  111 +IF NQ NE 0 THEN BEGIN
  112 + qsed=sed[indQ]
106 113 wavs=qsed.wave
107 114 ;=== Impose central wavelengths for photometric channels
108 115 ind=where(qsed.filter NE 'SPECTRUM',count)
109 116 IF count NE 0 THEN BEGIN
110   - wavs(ind)=dustem_filter2wav(qsed(ind).filter)
  117 + wavs[ind]=dustem_filter2wav(qsed[ind].filter)
111 118 ENDIF
112 119 ;=== define observations structure
113   - obs_st={instru_names:qsed.instru,filt_names:qsed.filter,wav:wavs, values:qsed.spec,sigma:qsed.error}
  120 + obs_st={instru_names:qsed.instru,filt_names:qsed.filter,wav:wavs, values:qsed.StokesQ,sigma:sqrt(qsed.SIGMAQQ)}
114 121 ;=== fill !dustem_data
115 122 !dustem_data.qsed = ptr_new(obs_st)
116 123 ;VGe
... ... @@ -120,7 +127,9 @@ IF keyword_set(qsed) THEN BEGIN
120 127 ;VGe
121 128 ENDIF
122 129  
123   -IF keyword_set(used) THEN BEGIN
  130 +IF NU NE 0 THEN BEGIN
  131 +;IF keyword_set(used) THEN BEGIN
  132 + used=sed[indU]
124 133 wavs=used.wave
125 134 ;=== Impose central wavelengths for photometric channels
126 135 ind=where(used.filter NE 'SPECTRUM',count)
... ... @@ -128,7 +137,7 @@ IF keyword_set(used) THEN BEGIN
128 137 wavs(ind)=dustem_filter2wav(used(ind).filter)
129 138 ENDIF
130 139 ;=== define observations structure
131   - obs_st={instru_names:used.instru,filt_names:used.filter,wav:wavs, values:used.spec,sigma:used.error}
  140 + obs_st={instru_names:used.instru,filt_names:used.filter,wav:wavs, values:used.StokesU,sigma:sqrt(qsed.SIGMAUU)}
132 141 ;=== fill !dustem_data
133 142 !dustem_data.used = ptr_new(obs_st)
134 143 ;VGe
... ... @@ -138,8 +147,6 @@ IF keyword_set(used) THEN BEGIN
138 147 ;VGe
139 148 ENDIF
140 149  
141   -
142   -
143 150 IF keyword_set(qext) THEN BEGIN
144 151 wavs=qext.wave
145 152 ;=== Impose central wavelengths for photometric channels
... ... @@ -176,9 +183,6 @@ IF keyword_set(uext) THEN BEGIN
176 183 ;VGe
177 184 ENDIF
178 185  
179   -
180   -
181   -
182 186 IF keyword_set(polfrac) THEN BEGIN
183 187 wavs=polfrac.wave
184 188 ;=== Impose central wavelengths for photometric channels
... ... @@ -194,7 +198,6 @@ ENDIF
194 198  
195 199 IF keyword_set(polext) THEN BEGIN
196 200 wavs=polext.wave
197   -
198 201 ;=== Impose central wavelengths for photometric channels
199 202 ind=where(polext.filter NE 'SPECTRUM',count)
200 203 IF count NE 0 THEN BEGIN
... ... @@ -220,9 +223,6 @@ IF keyword_set(ext) THEN BEGIN
220 223 ENDIF
221 224  
222 225  
223   -
224   -
225   -
226 226 !fit_rchi2_weight.sed=1.
227 227 !fit_rchi2_weight.ext=1.
228 228 if keyword_set(polext) then !fit_rchi2_weight.polext=1.
... ...