OV-GSO-DC / dustem-wrapper_idl

Commit 88f842727aafb0d6c88ffef44e9086a7879f6e15

Authored by Jean-Philippe Bernard
1 parent 5eebdedb
Exists in master

changed

Showing 8 changed files with 101 additions and 88 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_brute_force_fit.pro
  Diff comments   View file @88f8427
1 FUNCTION dustem_brute_force_fit,sed,table_name,filters,normalize=normalize,fact=fact,chi2=chi2,rchi2=rchi2,show_sed=show_sed 1 FUNCTION dustem_brute_force_fit,sed,table_name,filters,normalize=normalize,fact=fact,chi2=chi2,rchi2=rchi2,show_sed=show_sed
2 2
3 ;stop 3 ;stop
4 -;sed must be normalized to NH=1e21 H/cm2 4 +;sed must be normalized to NH=1e20 H/cm2
5 5
6 defsysv,'!dustem_grid',0,exist=exist 6 defsysv,'!dustem_grid',0,exist=exist
7 IF exist EQ 0 THEN BEGIN 7 IF exist EQ 0 THEN BEGIN
src/idl/dustem_compute_sed.pro
  Diff comments   View file @88f8427
1 FUNCTION dustem_compute_sed,p_dim,$ 1 FUNCTION dustem_compute_sed,p_dim,$
2 st=st,$ 2 st=st,$
3 sed_spec=sed_spec, $ 3 sed_spec=sed_spec, $
  4 + use_previous_fortran=use_previous_fortran, $
4 help=help 5 help=help
5 6
6 ;+ 7 ;+
@@ -20,7 +21,7 @@ FUNCTION dustem_compute_sed,p_dim,$ @@ -20,7 +21,7 @@ FUNCTION dustem_compute_sed,p_dim,$
20 ; p_dim = parameter values 21 ; p_dim = parameter values
21 ; 22 ;
22 ; OPTIONAL INPUT PARAMETERS: 23 ; OPTIONAL INPUT PARAMETERS:
23 -; st = Dustem output structure 24 +; st = Dustem Fortran output structure
24 ; 25 ;
25 ; OUTPUTS: 26 ; OUTPUTS:
26 ; sed = computed SED for filters in !dustem_data 27 ; sed = computed SED for filters in !dustem_data
@@ -30,6 +31,7 @@ FUNCTION dustem_compute_sed,p_dim,$ @@ -30,6 +31,7 @@ FUNCTION dustem_compute_sed,p_dim,$
30 ; sed_spec = dustem emission spectrum 31 ; sed_spec = dustem emission spectrum
31 ; 32 ;
32 ; ACCEPTED KEY-WORDS: 33 ; ACCEPTED KEY-WORDS:
  34 +; use_previous_fortran = if set, uses the output of the previous fortran run, and does not run the fortran.
33 ; help = If set, print this help 35 ; help = If set, print this help
34 ; 36 ;
35 ; COMMON BLOCKS: 37 ; COMMON BLOCKS:
@@ -59,7 +61,7 @@ IF keyword_set(help) THEN BEGIN @@ -59,7 +61,7 @@ IF keyword_set(help) THEN BEGIN
59 ENDIF 61 ENDIF
60 62
61 IF not keyword_set(st) THEN BEGIN 63 IF not keyword_set(st) THEN BEGIN
62 - dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st=st 64 + dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st=st,use_previous_fortran=use_previous_fortran
63 ENDIF 65 ENDIF
64 66
65 ; Convert into MJy/sr/1d20 67 ; Convert into MJy/sr/1d20
@@ -67,20 +69,22 @@ fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/st.sed.wav)*1.e20/1.e7 ;this is @@ -67,20 +69,22 @@ fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/st.sed.wav)*1.e20/1.e7 ;this is
67 SED_spec = st.sed.em_tot * fact 69 SED_spec = st.sed.em_tot * fact
68 70
69 ;ADDING PLUGIN TO SPECTRUM---------------- 71 ;ADDING PLUGIN TO SPECTRUM----------------
70 -scopes=tag_names((*!dustem_plugin)) 72 +Nplugins=n_elements(*!dustem_plugin)
  73 +scopes=(*!dustem_plugin).scope
71 74
72 IF scopes[0] NE 'NONE' THEN BEGIN 75 IF scopes[0] NE 'NONE' THEN BEGIN
73 -;IF ptr_valid(!dustem_plugin) THEN BEGIN  
74 - for i=0L,n_tags(*!dustem_plugin)-1 do begin  
75 - if total(strsplit((*(*!dustem_plugin).(i).scope),'+',/extract) eq 'ADD_SED') then SED_spec+=(*(*!dustem_plugin).(i).spec)[*,0]  
76 - endfor 76 + FOR i=0L,Nplugins-1 DO BEGIN
  77 + IF total(strsplit((*!dustem_plugin)[i].scope,'+',/extract) eq 'ADD_SED') THEN BEGIN
  78 + SED_spec+=(*(*!dustem_plugin)[i].spec)[*,0]
  79 + ENDIF
  80 + ENDFOR
77 ENDIF 81 ENDIF
78 ;------------------------------------------ 82 ;------------------------------------------
79 83
80 -if not isarray(SED_spec) THEN stop 84 +IF not isarray(SED_spec) THEN stop
81 85
82 ;COMPUTE THE MODEL SED 86 ;COMPUTE THE MODEL SED
83 - 87 +;JPB: problem there: we should be able to compute an SED even when there is no data in !dustem_data
84 dustem_sed = (*(*!dustem_data).sed).values * 0. 88 dustem_sed = (*(*!dustem_data).sed).values * 0.
85 89
86 IF !dustem_do_cc NE 0 AND !dustem_never_do_cc EQ 0 THEN BEGIN 90 IF !dustem_do_cc NE 0 AND !dustem_never_do_cc EQ 0 THEN BEGIN
src/idl/dustem_fit_intensity_example.pro
  Diff comments   View file @88f8427
@@ -257,7 +257,8 @@ if keyword_set(wait) then begin @@ -257,7 +257,8 @@ if keyword_set(wait) then begin
257 end 257 end
258 258
259 ;== INITIALISE DUSTEM 259 ;== INITIALISE DUSTEM
260 -dustem_init,model=use_model,polarization=use_polarization 260 +dustem_init,model=use_model,polarization=use_polarization,/show_plot
  261 +;stop
261 !dustem_nocatch=1 262 !dustem_nocatch=1
262 !dustem_verbose=use_verbose 263 !dustem_verbose=use_verbose
263 IF keyword_set(noobj) THEN !dustem_noobj=1 264 IF keyword_set(noobj) THEN !dustem_noobj=1
@@ -303,7 +304,7 @@ tit='FIT INTENSITY EXAMPLE' ; plot title @@ -303,7 +304,7 @@ tit='FIT INTENSITY EXAMPLE' ; plot title
303 ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') ; y-axis title 304 ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') ; y-axis title
304 xtit=textoidl('\lambda (\mum)') ; x-axis title 305 xtit=textoidl('\lambda (\mum)') ; x-axis title
305 306
306 - 307 +;stop
307 ;=== RUN THE FIT 308 ;=== RUN THE FIT
308 t1=systime(0,/sec) 309 t1=systime(0,/sec)
309 res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $ 310 res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $
src/idl/dustem_fit_intensity_phangs_example.pro
  Diff comments   View file @88f8427
@@ -10,7 +10,7 @@ PRO dustem_fit_intensity_phangs_example,model=model $ @@ -10,7 +10,7 @@ PRO dustem_fit_intensity_phangs_example,model=model $
10 10
11 ;+ 11 ;+
12 ; NAME: 12 ; NAME:
13 -; dustem_fit_intensity_example 13 +; dustem_fit_intensity_phangs_example
14 ; 14 ;
15 ; PURPOSE: 15 ; PURPOSE:
16 ; This routine is an example of how to fit an observational SED 16 ; This routine is an example of how to fit an observational SED
@@ -169,7 +169,9 @@ if keyword_set(wait) then begin @@ -169,7 +169,9 @@ if keyword_set(wait) then begin
169 end 169 end
170 170
171 ;== INITIALISE DUSTEM 171 ;== INITIALISE DUSTEM
172 -dustem_init,model=use_model,polarization=use_polarization 172 +;!dustem_show_plot=1
  173 +dustem_init,model=use_model,polarization=use_polarization,/show_plot
  174 +;stop
173 !dustem_nocatch=1 175 !dustem_nocatch=1
174 !dustem_verbose=use_verbose 176 !dustem_verbose=use_verbose
175 IF keyword_set(noobj) THEN !dustem_noobj=1 177 IF keyword_set(noobj) THEN !dustem_noobj=1
src/idl/dustem_plugin_phangs_stellar_continuum.pro
  Diff comments   View file @88f8427
@@ -31,6 +31,7 @@ FUNCTION dustem_plugin_phangs_stellar_continuum,key=key $ @@ -31,6 +31,7 @@ FUNCTION dustem_plugin_phangs_stellar_continuum,key=key $
31 ; paramtag = plugin parameter names as strings 31 ; paramtag = plugin parameter names as strings
32 ; ACCEPTED KEY-WORDS: 32 ; ACCEPTED KEY-WORDS:
33 ; help = if set, print this help 33 ; help = if set, print this help
  34 +; method = SSP used. can be EMILES or CB19. default='EMILES'
34 ; COMMON BLOCKS: 35 ; COMMON BLOCKS:
35 ; None 36 ; None
36 ; SIDE EFFECTS: 37 ; SIDE EFFECTS:
@@ -60,31 +61,25 @@ IF keyword_set(scope) THEN BEGIN @@ -60,31 +61,25 @@ IF keyword_set(scope) THEN BEGIN
60 goto, the_scope 61 goto, the_scope
61 ENDIF 62 ENDIF
62 63
63 -;IF keyword_set(paramtag) THEN BEGIN  
64 -; out=0  
65 -; goto, the_paramtag  
66 -;ENDIF  
67 -  
68 -;stop  
69 -  
70 -;default values of input parameters 64 +;use_method='EMILES'
  65 +use_method='CB19'
71 66
72 age_values=[0.03, 0.05, 0.08, 0.15, 0.25, 0.40, 0.60, 1.0, 1.75, 3.0, 5.0, 8.5, 13.5] 67 age_values=[0.03, 0.05, 0.08, 0.15, 0.25, 0.40, 0.60, 1.0, 1.75, 3.0, 5.0, 8.5, 13.5]
73 metalicity_values=[-1.48630, -0.961400, -0.351200, +0.0600000, +0.255900, +0.397100] 68 metalicity_values=[-1.48630, -0.961400, -0.351200, +0.0600000, +0.255900, +0.397100]
74 69
75 ;read template info 70 ;read template info
76 -info_st=dustem_read_emiles_stellar_templates(age_values,metalicity_values,/info_only)  
77 -;stop  
78 -age_values_emiles=info_st.age  
79 -order=sort(age_values_emiles)  
80 -age_values_emiles=age_values_emiles[order]  
81 -un=uniq(age_values_emiles)  
82 -age_values_emiles=age_values_emiles[un]  
83 -met_values_emiles=info_st.z  
84 -order=sort(met_values_emiles)  
85 -met_values_emiles=met_values_emiles[order]  
86 -un=uniq(met_values_emiles)  
87 -met_values_emiles=met_values_emiles[un] 71 +;==== This is actually not used
  72 +;info_st=dustem_read_emiles_stellar_templates(age_values,metalicity_values,/info_only)
  73 +;age_values_emiles=info_st.age
  74 +;order=sort(age_values_emiles)
  75 +;age_values_emiles=age_values_emiles[order]
  76 +;un=uniq(age_values_emiles)
  77 +;age_values_emiles=age_values_emiles[un]
  78 +;met_values_emiles=info_st.z
  79 +;order=sort(met_values_emiles)
  80 +;met_values_emiles=met_values_emiles[order]
  81 +;un=uniq(met_values_emiles)
  82 +;met_values_emiles=met_values_emiles[un]
88 83
89 ;==== This does not work because the info files has many more entries than there are available ssps. 84 ;==== This does not work because the info files has many more entries than there are available ssps.
90 ;age_values=age_values_emiles 85 ;age_values=age_values_emiles
@@ -132,8 +127,18 @@ IF exist EQ 0 THEN BEGIN @@ -132,8 +127,18 @@ IF exist EQ 0 THEN BEGIN
132 defsysv,'!dustem_plugin_phangs_stellar_continuum',st 127 defsysv,'!dustem_plugin_phangs_stellar_continuum',st
133 ENDIF 128 ENDIF
134 IF not ptr_valid(!dustem_plugin_phangs_stellar_continuum.ssps) THEN BEGIN 129 IF not ptr_valid(!dustem_plugin_phangs_stellar_continuum.ssps) THEN BEGIN
135 - message,'Initializing SSP templates',/continue  
136 - template_flux=dustem_read_emiles_stellar_templates(age_values,metalicity_values,template_wav=template_wav,Mstars=Mstars) 130 + message,'Initializing SSP templates',/info
  131 + CASE use_method of
  132 + 'EMILES': BEGIN
  133 + template_flux=dustem_read_emiles_stellar_templates(age_values,metalicity_values,template_wav=template_wav,Mstars=Mstars)
  134 + END
  135 + 'CB19': BEGIN
  136 + ;== this is just to get the Mstars
  137 + bidon=dustem_read_emiles_stellar_templates(age_values,metalicity_values,template_wav=template_wav,Mstars=Mstars)
  138 + ;== This is to read in the CB19 templates
  139 + template_flux=dustem_read_cb19_stellar_templates(age_values=age_values,metalicity_values=metalicity_values,template_wav=template_wav,Mstars=Mstars)
  140 + END
  141 + ENDCASE
137 !dustem_plugin_phangs_stellar_continuum.ssps=ptr_new(template_flux) 142 !dustem_plugin_phangs_stellar_continuum.ssps=ptr_new(template_flux)
138 !dustem_plugin_phangs_stellar_continuum.wavs=ptr_new(template_wav) 143 !dustem_plugin_phangs_stellar_continuum.wavs=ptr_new(template_wav)
139 !dustem_plugin_phangs_stellar_continuum.mstars=ptr_new(Mstars) 144 !dustem_plugin_phangs_stellar_continuum.mstars=ptr_new(Mstars)
@@ -157,7 +162,7 @@ FOR i=2L,Nparam-1 DO BEGIN @@ -157,7 +162,7 @@ FOR i=2L,Nparam-1 DO BEGIN
157 IF count NE 0 THEN BEGIN 162 IF count NE 0 THEN BEGIN
158 ssp=*(*!dustem_plugin_phangs_stellar_continuum.ssps)[i-2] 163 ssp=*(*!dustem_plugin_phangs_stellar_continuum.ssps)[i-2]
159 output[ind,0]=output[ind,0]+paramvalues[i]*Mstar*interpol(ssp,wavs,lambir[ind]) 164 output[ind,0]=output[ind,0]+paramvalues[i]*Mstar*interpol(ssp,wavs,lambir[ind])
160 - message,'Mstar='+strtrim(Mstar,2)+' paramvalue='+strtrim(paramvalues[i],2),/continue 165 + message,'Mstar='+strtrim(Mstar,2)+' paramvalue='+strtrim(paramvalues[i],2),/info
161 ENDIF 166 ENDIF
162 ENDIF 167 ENDIF
163 ENDFOR 168 ENDFOR
@@ -178,12 +183,12 @@ fact3=1./NHref ;[cm2/H] to go from ergs/s/cm2 to ergs/s/H @@ -178,12 +183,12 @@ fact3=1./NHref ;[cm2/H] to go from ergs/s/cm2 to ergs/s/H
178 fact4=1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/lambir)*1.e20/1.e7 ;This goes from ergs/s/H to MJy/sr (see dustem_plot_dataset) 183 fact4=1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/lambir)*1.e20/1.e7 ;This goes from ergs/s/H to MJy/sr (see dustem_plot_dataset)
179 fact=fact2*fact3*fact4 184 fact=fact2*fact3*fact4
180 185
181 -print,'flux(1mic)=',interpol(flux,lambir,1.),' ergs/s/cm2/AA'  
182 -print,'fact2=',interpol(fact2,lambir,1.)  
183 -print,'fact3=',fact3  
184 -print,'fact4=',interpol(fact4,lambir,1.)  
185 -print,'fact=',interpol(fact,lambir,1.)  
186 -print,'flux(1mic)=',interpol(flux,lambir,1.)*interpol(fact,lambir,1.),' MJy/sr' 186 +;print,'flux(1mic)=',interpol(flux,lambir,1.),' ergs/s/cm2/AA'
  187 +;print,'fact2=',interpol(fact2,lambir,1.)
  188 +;print,'fact3=',fact3
  189 +;print,'fact4=',interpol(fact4,lambir,1.)
  190 +;print,'fact=',interpol(fact,lambir,1.)
  191 +;print,'flux(1mic)=',interpol(flux,lambir,1.)*interpol(fact,lambir,1.),' MJy/sr'
187 192
188 ;goto,no_unred 193 ;goto,no_unred
189 194
@@ -213,7 +218,7 @@ output[*,0]=output[*,0]*fact @@ -213,7 +218,7 @@ output[*,0]=output[*,0]*fact
213 output[*,1]=output[*,1]*fact 218 output[*,1]=output[*,1]*fact
214 output[*,2]=output[*,2]*fact 219 output[*,2]=output[*,2]*fact
215 220
216 -print,'output=',interpol(output[*,0],lambir,1.),' MJy/sr' 221 +;print,'output=',interpol(output[*,0],lambir,1.),' MJy/sr'
217 ;stop 222 ;stop
218 223
219 the_scope: 224 the_scope:
src/idl/dustem_read_emiles_stellar_templates.pro
  Diff comments   View file @88f8427
@@ -39,7 +39,7 @@ FUNCTION dustem_read_emiles_stellar_templates,age_values, $ @@ -39,7 +39,7 @@ FUNCTION dustem_read_emiles_stellar_templates,age_values, $
39 ; PROCEDURES AND SUBROUTINES USED: 39 ; PROCEDURES AND SUBROUTINES USED:
40 ; 40 ;
41 ; EXAMPLES 41 ; EXAMPLES
42 -; 42 +; st=dustem_read_emiles_stellar_templates(age_values,Z_values)
43 ; MODIFICATION HISTORY: 43 ; MODIFICATION HISTORY:
44 ; Written by JPB June-2023 44 ; Written by JPB June-2023
45 ; Evolution details on the DustEMWrap gitlab. 45 ; Evolution details on the DustEMWrap gitlab.
@@ -91,8 +91,8 @@ metalicity_values_4str=string(metalicity_values_4str,format='(F5.2)') @@ -91,8 +91,8 @@ metalicity_values_4str=string(metalicity_values_4str,format='(F5.2)')
91 FOR i=0L,n_elements(metalicity_values_4str)-1 DO metalicity_values_4str[i]=textoidl_str_replace(metalicity_values_4str[i],' ','p') 91 FOR i=0L,n_elements(metalicity_values_4str)-1 DO metalicity_values_4str[i]=textoidl_str_replace(metalicity_values_4str[i],' ','p')
92 FOR i=0L,n_elements(metalicity_values_4str)-1 DO metalicity_values_4str[i]=textoidl_str_replace(metalicity_values_4str[i],'-','m') 92 FOR i=0L,n_elements(metalicity_values_4str)-1 DO metalicity_values_4str[i]=textoidl_str_replace(metalicity_values_4str[i],'-','m')
93 93
94 -print,age_values  
95 -print,metalicity_values 94 +;print,age_values
  95 +;print,metalicity_values
96 ;stop 96 ;stop
97 97
98 ;Ntemplates=n_elements(paramvalues)-1 98 ;Ntemplates=n_elements(paramvalues)-1
@@ -119,11 +119,11 @@ FOR i=1L,Ntemplates DO BEGIN @@ -119,11 +119,11 @@ FOR i=1L,Ntemplates DO BEGIN
119 met=metalicity_values[imet] 119 met=metalicity_values[imet]
120 met_str=metalicity_values_4str[imet] 120 met_str=metalicity_values_4str[imet]
121 age_str=age_values_4str[iage] 121 age_str=age_values_4str[iage]
122 - message,'reading template for Age '+age_str+' Met '+met_str,/continue 122 + message,'reading template for Age '+age_str+' Met '+met_str,/info
123 file=dir_templates+file_start+'Z'+met_str+'T'+age_str+file_end 123 file=dir_templates+file_start+'Z'+met_str+'T'+age_str+file_end
124 fi=file_info(file) 124 fi=file_info(file)
125 IF fi.exists THEN BEGIN 125 IF fi.exists THEN BEGIN
126 - message,'Reading SSP '+file,/continue 126 + message,'Reading SSP '+file,/info
127 spec=readfits(file,h,2) 127 spec=readfits(file,h,2)
128 templates[ii]=ptr_new(spec) 128 templates[ii]=ptr_new(spec)
129 make_axis,h,wavs 129 make_axis,h,wavs
src/idl/dustem_sed_extractor.pro
  Diff comments   View file @88f8427
1 -FUNCTION dustem_sed_extractor,maps $  
2 - ,index_roi $  
3 - ,filters, $ 1 +FUNCTION dustem_sed_extractor,maps, $
  2 + index_roi, $
  3 + filters, $
4 comments=comments, $ 4 comments=comments, $
5 maps_order=maps_order, $ 5 maps_order=maps_order, $
6 method=method, $ 6 method=method, $
@@ -25,7 +25,7 @@ FUNCTION dustem_sed_extractor,maps $ @@ -25,7 +25,7 @@ FUNCTION dustem_sed_extractor,maps $
25 ; sed=sed_extractor(maps,index,filters[,comments=][,maps_order=][,method=][,mask_I_undef=][,/total_intensity_only]) 25 ; sed=sed_extractor(maps,index,filters[,comments=][,maps_order=][,method=][,mask_I_undef=][,/total_intensity_only])
26 ; 26 ;
27 ; INPUTS: 27 ; INPUTS:
28 -; maps : a set of Nmaps maps of same dimension [Nx,Ny,Nmaps,Ndata_sets] 28 +; maps : a set of Nmaps maps of same dimension [Nx,Ny,Nmaps,Ndata_sets]. Caution, maps are changed if set_undefined is set !
29 ; index_roi : spatial index in Nx,Ny where the SED is to be extracted from (ROI=region of interest) 29 ; index_roi : spatial index in Nx,Ny where the SED is to be extracted from (ROI=region of interest)
30 ; filters : DustEMWrap filters corresponding to each dataset [Ndata_sets] 30 ; filters : DustEMWrap filters corresponding to each dataset [Ndata_sets]
31 ; 31 ;
@@ -117,9 +117,10 @@ ENDIF ELSE BEGIN @@ -117,9 +117,10 @@ ENDIF ELSE BEGIN
117 ENDIF 117 ENDIF
118 ENDELSE 118 ENDELSE
119 119
120 -used_maps=maps 120 +;used_maps=maps
121 121
122 IF keyword_set(set_undefined) THEN BEGIN 122 IF keyword_set(set_undefined) THEN BEGIN
  123 + used_maps=maps
123 ;===== do a map of saturated stokesI pixels 124 ;===== do a map of saturated stokesI pixels
124 mask_I_undef=maps[*,*,0,0] 125 mask_I_undef=maps[*,*,0,0]
125 mask_I_undef[*]=0 126 mask_I_undef[*]=0
@@ -132,8 +133,6 @@ IF keyword_set(set_undefined) THEN BEGIN @@ -132,8 +133,6 @@ IF keyword_set(set_undefined) THEN BEGIN
132 mask_I_undef[ind]=1 133 mask_I_undef[ind]=1
133 ENDIF 134 ENDIF
134 ENDFOR 135 ENDFOR
135 -  
136 - ;used_maps=maps  
137 ;set saturated stokesI pixels to la_undef in 136 ;set saturated stokesI pixels to la_undef in
138 indbad=where(mask_I_undef EQ 1,count) 137 indbad=where(mask_I_undef EQ 1,count)
139 IF count NE 0 THEN BEGIN 138 IF count NE 0 THEN BEGIN
@@ -146,34 +145,36 @@ IF keyword_set(set_undefined) THEN BEGIN @@ -146,34 +145,36 @@ IF keyword_set(set_undefined) THEN BEGIN
146 ENDFOR 145 ENDFOR
147 ENDFOR 146 ENDFOR
148 ENDIF 147 ENDIF
  148 + maps=used_maps
  149 + used_maps=0
149 ENDIF 150 ENDIF
150 151
151 CASE use_method OF 152 CASE use_method OF
152 'MEAN':BEGIN 153 'MEAN':BEGIN
153 - FOR j=0L,Nfilt-1 DO BEGIN  
154 - this_map=used_maps[*,*,indI,j] 154 + FOR j=0L,Nfilt-1 DO BEGIN
  155 + this_map=maps[*,*,indI,j]
155 sed[j].stokesI=la_mean(this_map[index_roi]) 156 sed[j].stokesI=la_mean(this_map[index_roi])
156 - this_map=used_maps[*,*,indII,j] 157 + this_map=maps[*,*,indII,j]
157 sed[j].sigmaII=la_div(la_mean(this_map[index_roi]),count_roi) 158 sed[j].sigmaII=la_div(la_mean(this_map[index_roi]),count_roi)
158 IF not keyword_set(total_intensity_only) THEN BEGIN 159 IF not keyword_set(total_intensity_only) THEN BEGIN
159 - this_map=used_maps[*,*,indQ,j] 160 + this_map=maps[*,*,indQ,j]
160 sed[j].stokesQ=la_mean(this_map[index_roi]) 161 sed[j].stokesQ=la_mean(this_map[index_roi])
161 - this_map=used_maps[*,*,indU,j] 162 + this_map=maps[*,*,indU,j]
162 sed[j].stokesU=la_mean(this_map[index_roi]) 163 sed[j].stokesU=la_mean(this_map[index_roi])
163 - this_map=used_maps[*,*,indQQ,j] 164 + this_map=maps[*,*,indQQ,j]
164 sed[j].sigmaQQ=la_div(la_mean(this_map[index_roi]),count_roi) 165 sed[j].sigmaQQ=la_div(la_mean(this_map[index_roi]),count_roi)
165 - this_map=used_maps[*,*,indUU,j] 166 + this_map=maps[*,*,indUU,j]
166 sed[j].sigmaUU=la_div(la_mean(this_map[index_roi]),count_roi) 167 sed[j].sigmaUU=la_div(la_mean(this_map[index_roi]),count_roi)
167 168
168 - this_map=used_maps[*,*,indIQ,j] 169 + this_map=maps[*,*,indIQ,j]
169 sed[j].sigmaIQ=la_div(la_mean(this_map[index_roi]),count_roi) 170 sed[j].sigmaIQ=la_div(la_mean(this_map[index_roi]),count_roi)
170 - this_map=used_maps[*,*,indIU,j] 171 + this_map=maps[*,*,indIU,j]
171 sed[j].sigmaIU=la_div(la_mean(this_map[index_roi]),count_roi) 172 sed[j].sigmaIU=la_div(la_mean(this_map[index_roi]),count_roi)
172 - this_map=used_maps[*,*,indQU,j] 173 + this_map=maps[*,*,indQU,j]
173 sed[j].sigmaQU=la_div(la_mean(this_map[index_roi]),count_roi) 174 sed[j].sigmaQU=la_div(la_mean(this_map[index_roi]),count_roi)
174 ENDIF 175 ENDIF
175 - ENDFOR  
176 - END 176 + ENDFOR
  177 + END
177 'MEAN_ONOFF': BEGIN ;This is mean ON - mean MEAN_ONOFF 178 'MEAN_ONOFF': BEGIN ;This is mean ON - mean MEAN_ONOFF
178 IF not keyword_set(index_off) THEN BEGIN 179 IF not keyword_set(index_off) THEN BEGIN
179 message,'index_off keyword must be set with method '+use_method,/continue 180 message,'index_off keyword must be set with method '+use_method,/continue
@@ -182,36 +183,36 @@ CASE use_method OF @@ -182,36 +183,36 @@ CASE use_method OF
182 count_off=n_elements(index_off) 183 count_off=n_elements(index_off)
183 comments=[comments,'OFF has '+strtrim(count_off,2)+' map pixels'] 184 comments=[comments,'OFF has '+strtrim(count_off,2)+' map pixels']
184 FOR j=0L,Nfilt-1 DO BEGIN 185 FOR j=0L,Nfilt-1 DO BEGIN
185 - this_map=used_maps[*,*,indI,j] 186 + this_map=maps[*,*,indI,j]
186 sed[j].stokesI=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off])) 187 sed[j].stokesI=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off]))
187 - this_map=used_maps[*,*,indII,j] 188 + this_map=maps[*,*,indII,j]
188 sig_on=la_div(la_mean(this_map[index_roi]),count_roi) 189 sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
189 sig_off=la_div(la_mean(this_map[index_off]),count_off) 190 sig_off=la_div(la_mean(this_map[index_off]),count_off)
190 sed[j].sigmaII=la_add(sig_on,sig_off) 191 sed[j].sigmaII=la_add(sig_on,sig_off)
191 IF not keyword_set(total_intensity_only) THEN BEGIN 192 IF not keyword_set(total_intensity_only) THEN BEGIN
192 - this_map=used_maps[*,*,indQ,j] 193 + this_map=maps[*,*,indQ,j]
193 sed[j].stokesQ=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off])) 194 sed[j].stokesQ=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off]))
194 - this_map=used_maps[*,*,indU,j] 195 + this_map=maps[*,*,indU,j]
195 sed[j].stokesU=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off])) 196 sed[j].stokesU=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off]))
196 ;variances 197 ;variances
197 - this_map=used_maps[*,*,indQQ,j] 198 + this_map=maps[*,*,indQQ,j]
198 sig_on=la_div(la_mean(this_map[index_roi]),count_roi) 199 sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
199 sig_off=la_div(la_mean(this_map[index_off]),count_off) 200 sig_off=la_div(la_mean(this_map[index_off]),count_off)
200 sed[j].sigmaQQ=la_add(sig_on,sig_off) 201 sed[j].sigmaQQ=la_add(sig_on,sig_off)
201 - this_map=used_maps[*,*,indUU,j] 202 + this_map=maps[*,*,indUU,j]
202 sig_on=la_div(la_mean(this_map[index_roi]),count_roi) 203 sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
203 sig_off=la_div(la_mean(this_map[index_off]),count_off) 204 sig_off=la_div(la_mean(this_map[index_off]),count_off)
204 sed[j].sigmaUU=la_add(sig_on,sig_off) 205 sed[j].sigmaUU=la_add(sig_on,sig_off)
205 ;co-variances 206 ;co-variances
206 - this_map=used_maps[*,*,indIQ,j] 207 + this_map=maps[*,*,indIQ,j]
207 sig_on=la_div(la_mean(this_map[index_roi]),count_roi) 208 sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
208 sig_off=la_div(la_mean(this_map[index_off]),count_off) 209 sig_off=la_div(la_mean(this_map[index_off]),count_off)
209 sed[j].sigmaIQ=la_add(sig_on,sig_off) 210 sed[j].sigmaIQ=la_add(sig_on,sig_off)
210 - this_map=used_maps[*,*,indIU,j] 211 + this_map=maps[*,*,indIU,j]
211 sig_on=la_div(la_mean(this_map[index_roi]),count_roi) 212 sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
212 sig_off=la_div(la_mean(this_map[index_off]),count_off) 213 sig_off=la_div(la_mean(this_map[index_off]),count_off)
213 sed[j].sigmaIU=la_add(sig_on,sig_off) 214 sed[j].sigmaIU=la_add(sig_on,sig_off)
214 - this_map=used_maps[*,*,indQU,j] 215 + this_map=maps[*,*,indQU,j]
215 sig_on=la_div(la_mean(this_map[index_roi]),count_roi) 216 sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
216 sig_off=la_div(la_mean(this_map[index_off]),count_off) 217 sig_off=la_div(la_mean(this_map[index_off]),count_off)
217 sed[j].sigmaQU=la_add(sig_on,sig_off) 218 sed[j].sigmaQU=la_add(sig_on,sig_off)
@@ -232,21 +233,21 @@ CASE use_method OF @@ -232,21 +233,21 @@ CASE use_method OF
232 message,'reference filter '+reference_filter+' not found in provided filter list',/continue 233 message,'reference filter '+reference_filter+' not found in provided filter list',/continue
233 stop 234 stop
234 ENDIF 235 ENDIF
235 - mapI_ref=used_maps[*,*,indI,indf]  
236 - mapQ_ref=used_maps[*,*,indQ,indf]  
237 - mapU_ref=used_maps[*,*,indU,indf]  
238 - mapII_ref=used_maps[*,*,indII,indf]  
239 - mapQQ_ref=used_maps[*,*,indQQ,indf]  
240 - mapUU_ref=used_maps[*,*,indUU,indf] 236 + mapI_ref=maps[*,*,indI,indf]
  237 + mapQ_ref=maps[*,*,indQ,indf]
  238 + mapU_ref=maps[*,*,indU,indf]
  239 + mapII_ref=maps[*,*,indII,indf]
  240 + mapQQ_ref=maps[*,*,indQQ,indf]
  241 + mapUU_ref=maps[*,*,indUU,indf]
241 Iref=la_mean(mapI_ref[index_roi]) 242 Iref=la_mean(mapI_ref[index_roi])
242 Qref=la_mean(mapQ_ref[index_roi]) 243 Qref=la_mean(mapQ_ref[index_roi])
243 Uref=la_mean(mapU_ref[index_roi]) 244 Uref=la_mean(mapU_ref[index_roi])
244 FOR j=0L,Nfilt-1 DO BEGIN 245 FOR j=0L,Nfilt-1 DO BEGIN
245 ;==== do the correlation in Stokes I 246 ;==== do the correlation in Stokes I
246 x=mapI_ref[index_roi] 247 x=mapI_ref[index_roi]
247 - y=(used_maps[*,*,indI,j])[index_roi] 248 + y=(maps[*,*,indI,j])[index_roi]
248 sx=la_power(mapII_ref[index_roi],0.5) 249 sx=la_power(mapII_ref[index_roi],0.5)
249 - sy=la_power((used_maps[*,*,indI,j])[index_roi],0.5) 250 + sy=la_power((maps[*,*,indI,j])[index_roi],0.5)
250 ind_good=where(x NE la_undef() AND y NE la_undef() AND sx NE la_undef() and sy NE la_undef(),count_good) 251 ind_good=where(x NE la_undef() AND y NE la_undef() AND sx NE la_undef() and sy NE la_undef(),count_good)
251 comments=[comments,'I corelation has '+strtrim(count_good,2)+' map pixels'] 252 comments=[comments,'I corelation has '+strtrim(count_good,2)+' map pixels']
252 IF count_good NE 0 THEN BEGIN 253 IF count_good NE 0 THEN BEGIN
@@ -268,9 +269,9 @@ CASE use_method OF @@ -268,9 +269,9 @@ CASE use_method OF
268 ENDIF 269 ENDIF
269 ;Do the correlation between Q maps 270 ;Do the correlation between Q maps
270 x=[mapQ_ref[index_roi],mapU_ref[index_roi]] 271 x=[mapQ_ref[index_roi],mapU_ref[index_roi]]
271 - y=[(used_maps[*,*,indQ,j])[index_roi],(used_maps[*,*,indU,j])[index_roi]] 272 + y=[(maps[*,*,indQ,j])[index_roi],(maps[*,*,indU,j])[index_roi]]
272 sx=la_power([mapQQ_ref[index_roi],mapUU_ref[index_roi]],0.5) 273 sx=la_power([mapQQ_ref[index_roi],mapUU_ref[index_roi]],0.5)
273 - sy=la_power([(used_maps[*,*,indQ,j])[index_roi],(used_maps[*,*,indU,j])[index_roi]],0.5) 274 + sy=la_power([(maps[*,*,indQ,j])[index_roi],(maps[*,*,indU,j])[index_roi]],0.5)
274 ind_good=where(x NE la_undef() AND y NE la_undef() AND sx NE la_undef() and sy NE la_undef(),count_good) 275 ind_good=where(x NE la_undef() AND y NE la_undef() AND sx NE la_undef() and sy NE la_undef(),count_good)
275 comments=[comments,'QU corelation has '+strtrim(count_good,2)+' map pixels'] 276 comments=[comments,'QU corelation has '+strtrim(count_good,2)+' map pixels']
276 IF count_good NE 0 THEN BEGIN 277 IF count_good NE 0 THEN BEGIN
src/idl/dustem_write_chrg.pro
  Diff comments   View file @88f8427
@@ -94,7 +94,7 @@ FOR i=0L,Nst-1 DO BEGIN @@ -94,7 +94,7 @@ FOR i=0L,Nst-1 DO BEGIN
94 fv=str_sep(ffile,'/') 94 fv=str_sep(ffile,'/')
95 file=dir+fv(n_elements(fv)-1) 95 file=dir+fv(n_elements(fv)-1)
96 ;Open file 96 ;Open file
97 - message,'Will write CHRG file '+file,/continue 97 + message,'Will write CHRG file '+file,/info
98 openw,unit,file,/get_lun 98 openw,unit,file,/get_lun
99 ;Print comments 99 ;Print comments
100 IF Ncomments NE 0 THEN BEGIN 100 IF Ncomments NE 0 THEN BEGIN