OV-GSO-DC / dustem-wrapper_idl

Commit b7dca888ba4c971de4b88b0d6e9626c3c42abe6e

Authored by Jean-Philippe Bernard
1 parent 2fb72ad4
Exists in master

improved to recover all system variables

Showing 6 changed files with 109 additions and 69 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_fillup_systvar_from_fits.pro
  Diff comments   View file @b7dca88
1   -PRO dustem_fillup_systvar_from_fits,syst_var,sed,str_input_SED,used_pol
  1 +PRO dustem_fillup_systvar_from_fits,syst_var,str_input_SED,used_pol
  2 +
  3 + Nsed=n_elements(str_input_SED)
  4 + sed=dustem_initialize_internal_sed(Nsed)
2 5  
3 6 (syst_var).sed=ptr_new(sed)
4 7  
... ... @@ -10,17 +13,49 @@ PRO dustem_fillup_systvar_from_fits,syst_var,sed,str_input_SED,used_pol
10 13  
11 14 IF used_pol THEN BEGIN
12 15 ;stop
13   - syst_var.qsed=ptr_new(sed)
14   - syst_var.used=ptr_new(sed)
15   - (*syst_var.qsed).instru_names=dustem_filter2instru(str_input_SED.filter)
16   - (*syst_var.qsed).filt_names=strtrim(str_input_SED.filter,2)
17   - (*syst_var.qsed).wav=str_input_SED.wavelength
18   - (*syst_var.qsed).values=str_input_SED.STOKESQ
19   - (*syst_var.qsed).sigma=la_power(str_input_SED.varianceQQ,0.5)
20   - (*syst_var.used).instru_names=dustem_filter2instru(str_input_SED.filter)
21   - (*syst_var.used).filt_names=str_input_SED.filter
22   - (*syst_var.used).wav=str_input_SED.wavelength
23   - (*syst_var.used).values=str_input_SED.STOKESU
24   - (*syst_var.used).sigma=la_power(str_input_SED.varianceUU,0.5)
  16 + ;===== note: those should include only points with polarization data, not all sed points
  17 + indq=where(str_input_SED.STOKESQ NE la_undef(),Nqsed)
  18 + indu=where(str_input_SED.STOKESU NE la_undef(),Nused)
  19 + qsed=dustem_initialize_internal_sed(Nqsed)
  20 + used=dustem_initialize_internal_sed(Nused)
  21 + indpolar=where(str_input_SED.STOKESU NE la_undef() AND str_input_SED.STOKESQ NE la_undef(),Npolar)
  22 + psi_sed=dustem_initialize_internal_sed(Npolar)
  23 + smallp_sed=dustem_initialize_internal_sed(Npolar)
  24 + largep_sed=dustem_initialize_internal_sed(Npolar)
  25 +
  26 + syst_var.qsed=ptr_new(qsed)
  27 + syst_var.used=ptr_new(used)
  28 + (*syst_var.qsed).instru_names=dustem_filter2instru(str_input_SED[indq].filter)
  29 + (*syst_var.qsed).filt_names=strtrim(str_input_SED[indq].filter,2)
  30 + (*syst_var.qsed).wav=str_input_SED[indq].wavelength
  31 + (*syst_var.qsed).values=str_input_SED[indq].STOKESQ
  32 + (*syst_var.qsed).sigma=la_power(str_input_SED[indq].varianceQQ,0.5)
  33 + ;now U
  34 + (*syst_var.used).instru_names=dustem_filter2instru(str_input_SED[indu].filter)
  35 + (*syst_var.used).filt_names=str_input_SED[indu].filter
  36 + (*syst_var.used).wav=str_input_SED[indu].wavelength
  37 + (*syst_var.used).values=str_input_SED[indu].STOKESU
  38 + (*syst_var.used).sigma=la_power(str_input_SED[indu].varianceUU,0.5)
  39 + ;=== add up dependent fields
  40 + syst_var.psi_em=ptr_new(psi_sed)
  41 + syst_var.polfrac=ptr_new(smallp_sed)
  42 + syst_var.polsed=ptr_new(largep_sed)
  43 + polar_iqu2ippsi,str_input_SED[indpolar].STOKESI,str_input_SED[indpolar].STOKESQ,str_input_SED[indpolar].STOKESU,p,psi,largep=largep
  44 + (*syst_var.psi_em).values=psi
  45 + (*syst_var.polfrac).values=p
  46 + (*syst_var.polsed).values=largep
  47 + CIQ=str_input_SED[indpolar].VARIANCEII & CIQ[*]=0. ;because input SEDs don't include covariances, yet
  48 + CIU=str_input_SED[indpolar].VARIANCEII & CIU[*]=0. ;because input SEDs don't include covariances, yet
  49 + CQU=str_input_SED[indpolar].VARIANCEII & CQU[*]=0. ;because input SEDs don't include covariances, yet
  50 + polar_variance_iqu2ippsi,str_input_SED[indpolar].STOKESI $
  51 + ,str_input_SED[indpolar].STOKESQ $
  52 + ,str_input_SED[indpolar].STOKESU $
  53 + ,str_input_SED[indpolar].VARIANCEII $
  54 + ,str_input_SED[indpolar].VARIANCEQQ $
  55 + ,str_input_SED[indpolar].VARIANCEUU $
  56 + ,CIQ,CIU,CQU,variance_smallp,variance_psi,variance_largep,/lac,snr_largeP=snr_largeP
  57 + (*syst_var.psi_em).sigma=la_power(variance_psi,0.5)
  58 + (*syst_var.polfrac).sigma=la_power(variance_smallp,0.5)
  59 + (*syst_var.polsed).sigma=la_power(variance_largep,0.5)
25 60 ENDIF
26 61 END
27 62 \ No newline at end of file
... ...
src/idl/dustem_init.pro
  Diff comments   View file @b7dca88
... ... @@ -179,45 +179,9 @@ dustem_fit_st={data:ptr_new(), $ ;,wavelength:ptr_new(), ;because wavelength arr
179 179 }
180 180 defsysv, '!dustem_fit', ptr_new(dustem_fit_st) ;Data to fit
181 181  
182   -if keyword_set(polarization) then begin
183   -
184   - dustem_data_st ={ $ ;Structure contaning the data to fit
185   - sed: ptr_new(), $ ;
186   - ext: ptr_new(), $ ;
187   - qsed: ptr_new(), $ ;
188   - used: ptr_new(), $ ;
189   - polsed: ptr_new(), $ ;
190   - polfrac: ptr_new(), $ ;
191   - psi_em: ptr_new(), $ ;
192   - qext: ptr_new(), $ ;
193   - uext: ptr_new(), $ ;
194   - polext: ptr_new(), $ ;
195   - fpolext: ptr_new(), $ ;
196   - psi_ext: ptr_new() $ ;
197   -
198   - }
199   -
200   - defsysv, '!dustem_data', ptr_new(dustem_data_st)
201   -
202   - defsysv,'!dustem_show', ptr_new(dustem_data_st)
203   -
204   - rchi2_weight = {sed: 0.,ext: 0.,qsed:0.,used:0.,polsed:0.,polfrac:0.,psi_em:0.,qext:0.,uext:0.,polext:0.,psi_ext:0.}
205   -
206   -endif else begin
207   -
208   - dustem_data_st = { $
209   - sed: ptr_new(), $
210   - ext: ptr_new() $
211   - }
212   -
213   - defsysv, '!dustem_data', ptr_new(dustem_data_st)
214   -
215   -
216   - defsysv, '!dustem_show', ptr_new(dustem_data_st)
217   -
218   - rchi2_weight = {sed: 0.,ext: 0.}
219   -
220   -endelse
  182 +dustem_data_st = dustem_define_dustem_data(polarization=polarization,rchi2_weight =rchi2_weight)
  183 +defsysv, '!dustem_data', ptr_new(dustem_data_st)
  184 +defsysv,'!dustem_show', ptr_new(dustem_data_st)
221 185  
222 186 if !run_pol then begin
223 187 tagnames=tag_names(*!dustem_data)
... ...
src/idl/dustem_make_fits_predicted_ext.pro
  Diff comments   View file @b7dca88
1   -FUNCTION dustem_make_fits_predicted_ext,syst_var,dustem_predicted_ext,str_input_ext=str_input_ext,help=help
  1 +FUNCTION dustem_make_fits_predicted_ext ,syst_var $
  2 + ,dustem_predicted_ext $
  3 + ,str_input_ext=str_input_ext $
  4 + ,dustem_predicted_Qext=dustem_predicted_Qext, $
  5 + ,dustem_predicted_Uext=dustem_predicted_Uext, $
  6 + ,help=help
2 7  
3 8 ;+
4 9 ; NAME:
... ... @@ -14,7 +19,8 @@ FUNCTION dustem_make_fits_predicted_ext,syst_var,dustem_predicted_ext,str_input_
14 19 ; syst_var : dustem system variable (*!dustem_data or *!dustem_show)
15 20 ; dustem_predicted_ext : dustemwrap predicted extinction SED as output by dustem_compute_ext.pro
16 21 ; OPTIONAL INPUT PARAMETERS:
17   -; None
  22 +; dustem_predicted_Qext : Stokes Q extinction SED
  23 +; dustem_predicted_Uext : Stokes U extinction SED
18 24 ; OUTPUTS:
19 25 ; str_predicted_EXT : dustemwrap predicted extinction SED as used to store into fits file
20 26 ; OPTIONAL OUTPUT PARAMETERS:
... ...
src/idl/dustem_make_fits_predicted_sed.pro
  Diff comments   View file @b7dca88
1   -FUNCTION dustem_make_fits_predicted_sed,syst_var,dustem_predicted_sed,str_input_sed=str_input_sed,help=help
  1 +FUNCTION dustem_make_fits_predicted_sed ,syst_var $
  2 + ,dustem_predicted_sed $
  3 + ,str_input_sed=str_input_sed $
  4 + ,dustem_predicted_Qsed=dustem_predicted_Qsed $
  5 + ,dustem_predicted_Used=dustem_predicted_Used $
  6 + ,help=help
2 7  
3 8 ;+
4 9 ; NAME:
... ... @@ -12,9 +17,11 @@ FUNCTION dustem_make_fits_predicted_sed,syst_var,dustem_predicted_sed,str_input_
12 17 ; see dustem_write_fits_table.pro
13 18 ; INPUTS:
14 19 ; syst_var : dustem system variable (*!dustem_data or *!dustem_show)
  20 +; what = which tag of syst_var must be treated (can be sed, qsed, used, ...)
15 21 ; dustem_predicted_sed : dustemwrap predicted SED as output by dustem_compute_sed.pro
16 22 ; OPTIONAL INPUT PARAMETERS:
17   -; None
  23 +; dustem_predicted_Qsed : Stokes Q SED
  24 +; dustem_predicted_Used : Stokes U SED
18 25 ; OUTPUTS:
19 26 ; str_predicted_SED : dustemwrap predicted SED as used to store into fits file
20 27 ; OPTIONAL OUTPUT PARAMETERS:
... ... @@ -50,6 +57,11 @@ ENDIF
50 57 one_str_input_SED={FILTER:'',Wavelength:la_undef(4),STOKESI:la_undef(5),STOKESQ:la_undef(5),STOKESU:la_undef(5),varianceII:la_undef(5),varianceQQ:la_undef(5),varianceUU:la_undef(5)}
51 58 one_str_predicted_SED={FILTER:'',Wavelength:la_undef(4),STOKESI:la_undef(5),STOKESQ:la_undef(5),STOKESU:la_undef(5)}
52 59  
  60 +;tags=tag_names(syst_var)
  61 +;ind_what=where(tags EQ strupcase(what),count_what)
  62 +
  63 +;stop
  64 +
53 65 Nsed=n_elements((*syst_var.sed).filt_names)
54 66 str_input_SED=replicate(one_str_input_SED,Nsed)
55 67 str_input_SED.FILTER=(*syst_var.sed).filt_names
... ...
src/idl/dustem_read_fits_table.pro
  Diff comments   View file @b7dca88
... ... @@ -217,6 +217,8 @@ ENDIF ELSE BEGIN
217 217 dustem_init,model=0
218 218 ENDELSE
219 219  
  220 +defsysv,'!dustem_hcd',ptr_new(float(sxpar(header_input_sed,'NH')))
  221 +
220 222 !dustem_model=used_model
221 223 IF !dustem_version.version NE used_version THEN BEGIN
222 224 message,'Caution: fits file was not created with the same Dustemwrap version as your current version',/continue
... ... @@ -254,9 +256,9 @@ ENDFOR
254 256 ;==== This is to create (*!dustem_data).sed in case it does not exist already
255 257 ; it doesn't exist at all times because dustem_init has been run
256 258 IF unit_input_emission_sed NE -1 THEN BEGIN
257   - Nsed=n_elements(str_input_SED)
258   - sed=dustem_initialize_internal_sed(Nsed)
259   - dustem_fillup_systvar_from_fits,*!dustem_data,sed,str_input_SED,used_pol
  259 + ;Nsed=n_elements(str_input_SED)
  260 + ;sed=dustem_initialize_internal_sed(Nsed)
  261 + dustem_fillup_systvar_from_fits,*!dustem_data,str_input_SED,used_pol
260 262 ENDIF
261 263  
262 264 ;help,str_shown_SED
... ... @@ -264,15 +266,19 @@ ENDIF
264 266  
265 267 ;==== This is to create (*!dustem_show).sed in case it does not exist already
266 268 IF unit_dustem_shown_sed NE -1 THEN BEGIN
267   - Nsed=n_elements(str_shown_SED)
268   - sed=dustem_initialize_internal_sed(Nsed)
269   - dustem_fillup_systvar_from_fits,*!dustem_show,sed,str_shown_SED,used_pol
  269 + ;Nsed=n_elements(str_shown_SED)
  270 + ;sed=dustem_initialize_internal_sed(Nsed)
  271 + dustem_fillup_systvar_from_fits,*!dustem_show,str_shown_SED,used_pol
270 272 ENDIF
271 273  
272 274 IF unit_input_extinction_sed NE -1 THEN BEGIN
273   - Next=n_elements(str_input_EXT)
274   - ext=dustem_initialize_internal_sed(Next)
275   - dustem_fillup_systvar_from_fits,*!dustem_data,ext,str_input_EXT,used_pol
  275 + ;Next=n_elements(str_input_EXT)
  276 + ;ext=dustem_initialize_internal_sed(Next)
  277 + dustem_fillup_systvar_from_fits,*!dustem_data,str_input_EXT,used_pol
  278 +ENDIF
  279 +
  280 +IF unit_dustem_shown_extinction_sed NE -1 THEN BEGIN
  281 + dustem_fillup_systvar_from_fits,*!dustem_show,str_shown_EXT,used_pol
276 282 ENDIF
277 283  
278 284 ;=== Below is only to get the same output form as dustem_compute_sed
... ...
src/idl/dustem_write_fits_table.pro
  Diff comments   View file @b7dca88
... ... @@ -60,6 +60,8 @@ extinction_sed_present=ptr_valid((*!dustem_data).ext)
60 60 ;fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/st.sed.wav)*1.e20/1.e7 ;this is correct despite the way it is presented (1.e4/1e.7*1e.20=1.e17 which is the factor to convert from ergs to Mega Janskys)
61 61 ;SED_spec = st.sed.em_tot * fact
62 62  
  63 +;stop
  64 +
63 65 ;===== run dustem to get predicted spectrum and sed
64 66 IF emission_sed_present THEN BEGIN
65 67 dustem_predicted_sed=dustem_compute_sed(*(*!dustem_fit).current_param_values,st=dustem_st)
... ... @@ -71,8 +73,14 @@ IF emission_sed_present THEN BEGIN
71 73 dustem_predicted_Qsed=dustem_predicted_polsed[0]
72 74 dustem_predicted_Used=dustem_predicted_polsed[1]
73 75 ENDIF
74   - str_predicted_SED=dustem_make_fits_predicted_sed(*!dustem_data,dustem_predicted_sed,str_input_sed=str_input_sed)
75   - str_predicted_shown_SED=dustem_make_fits_predicted_sed(*!dustem_show,dustem_predicted_sed,str_input_sed=str_shown_sed)
  76 + str_predicted_SED=dustem_make_fits_predicted_sed(*!dustem_data,dustem_predicted_sed $
  77 + ,str_input_sed=str_input_sed $
  78 + ,dustem_predicted_Qsed=dustem_predicted_Qsed $
  79 + ,dustem_predicted_Used=dustem_predicted_Used)
  80 + str_predicted_shown_SED=dustem_make_fits_predicted_sed(*!dustem_show,dustem_predicted_sed $
  81 + ,str_input_sed=str_shown_sed $
  82 + ,dustem_predicted_Qsed=dustem_predicted_Qsed $
  83 + ,dustem_predicted_Used=dustem_predicted_Used)
76 84 ENDIF
77 85  
78 86 IF extinction_sed_present THEN BEGIN
... ... @@ -84,8 +92,14 @@ IF extinction_sed_present THEN BEGIN
84 92 dustem_predicted_Qext=dustem_predicted_polext[0]
85 93 dustem_predicted_Uext=dustem_predicted_polext[1]
86 94 ENDIF
87   - str_predicted_EXT=dustem_make_fits_predicted_ext(*!dustem_data,dustem_predicted_ext,str_input_ext=str_input_ext)
88   - str_predicted_shown_EXT=dustem_make_fits_predicted_ext(*!dustem_show,dustem_predicted_ext,str_input_ext=str_shown_ext)
  95 + str_predicted_EXT=dustem_make_fits_predicted_ext(*!dustem_data,dustem_predicted_ext $
  96 + ,dustem_predicted_Qext=dustem_predicted_Qext $
  97 + ,dustem_predicted_Uext=dustem_predicted_Uext $
  98 + ,str_input_ext=str_input_ext)
  99 + str_predicted_shown_EXT=dustem_make_fits_predicted_ext(*!dustem_show,dustem_predicted_ext $
  100 + ,dustem_predicted_Qext=dustem_predicted_Qext $
  101 + ,dustem_predicted_Uext=dustem_predicted_Uext $
  102 + ,str_input_ext=str_shown_ext)
89 103 ENDIF
90 104  
91 105 Ngrains=n_tags(dustem_st.sed)-2
... ... @@ -255,6 +269,9 @@ sxaddpar,header_input_sed,'TTYPE5','U','Observed Stokes U Intensity [??]'
255 269 sxaddpar,header_input_sed,'TTYPE6','VARIANCEII','Observed Total Stokes I Variance [??]'
256 270 sxaddpar,header_input_sed,'TTYPE7','VARIANCEQQ','Observed Stokes Q Variance [??]'
257 271 sxaddpar,header_input_sed,'TTYPE8','VARIANCEUU','Observed Stokes U Variance [??]'
  272 +sxaddpar,header_input_sed,'NH',*!dustem_hcd,'Column density [H/cm2]'
  273 +
  274 +;stop
258 275  
259 276 ;===== complement extension header for predicted SED
260 277 sxaddpar,header_predicted_SED,'TITLE','OUTPUT BEST FIT SED BY DUSTEMWRAP'
... ...