OV-GSO-DC / dustem-wrapper_idl

Commit ecfb8fb3839a73d255d943d4de9c858086245a6b

Authored by Annie Hughes
1 parent 26f12f12
Exists in master

updated comments for user example

Showing 2 changed files with 217 additions and 94 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_extract_sed_example.pro
  Diff comments   View file @ecfb8fb
1   -PRO test_dustem_sed_extractor
  1 +PRO dustem_extract_sed_example
  2 +
  3 +;+
  4 +; NAME:
  5 +; dustem_extract_sed_example
  6 +;
  7 +; PURPOSE:
  8 +; This routine contains a simple example of how to extract an SED to use with DustEMWrap.
  9 +;
  10 +; CATEGORY:
  11 +; DustEMWrap, Distributed, High-Level, User Example
  12 +;
  13 +; CALLING SEQUENCE:
  14 +; dustem_extract_sed_example
  15 +;
  16 +; INPUTS:
  17 +; None
  18 +;
  19 +; OPTIONAL INPUT PARAMETERS:
  20 +; None
  21 +;
  22 +; OUTPUTS:
  23 +; None
  24 +;
  25 +; OPTIONAL OUTPUT PARAMETERS:
  26 +;
  27 +; ACCEPTED KEY-WORDS:
  28 +; None
  29 +;
  30 +; COMMON BLOCKS:
  31 +; None
  32 +;
  33 +; SIDE EFFECTS:
  34 +; None
  35 +;
  36 +; RESTRICTIONS:
  37 +; The DustEM fortran code must be installed
  38 +; The DustEMWrap IDL code must be installed
  39 +;
  40 +; PROCEDURES AND SUBROUTINES USED:
  41 +;
  42 +; EXAMPLES
  43 +;
  44 +; MODIFICATION HISTORY:
  45 +; Written by JPB Jan 2022
  46 +; Evolution details on the DustEMWrap gitlab.
  47 +; See http://dustemwrap.irap.omp.eu/ for FAQ and help.
  48 +;-
  49 +
  50 +IF keyword_set(help) THEN BEGIN
  51 + doc_library,'dustem_extract_sed_example'
  52 + goto,the_end
  53 +END
2 54  
3 55 define_la_common
4   -dustem_init
  56 +use_model='DBP90'
  57 +use_polarization=0
  58 +dustem_init,model=use_model,polarization=use_polarization
5 59  
6   -Nmaps=2 ;This is total intensity and variances
  60 +Nmaptypes=2 ; number of map types in the structure. Here 2 (Stokes I, variance in I)
7 61 filters=['IRAS3','IRAS4','PACS1','PACS2','PACS3','PILOT1','SPIRE1','SPIRE2','SPIRE3','HFI3','HFI4','HFI5','HFI6']
8 62 Ndata_sets=n_elements(filters)
9 63  
10   -maps=fltarr(100,100,Nmaps,Ndata_sets)
  64 +maps=fltarr(100,100,Nmaptypes,Ndata_sets)
11 65 maps[*]=la_undef()
12 66  
  67 +;; Define the region of interest (ON) and OFF region
13 68 mask=maps[*,*,0,0]
14   -mask[5:35,5:35]=1
15   -mask[10:30,10:30]=mask[10:30,10:30]+1
  69 +mask[5:35,5:35]=1 ; OFF will be a guard band of 5 pixel around the ON defined below
  70 +mask[10:30,10:30]=mask[10:30,10:30]+1 ; ON source
16 71  
17 72 index_on=where(mask EQ 2,count_on)
18 73 index_off=where(mask EQ 1,count_off)
... ... @@ -21,34 +76,82 @@ indI=0
21 76 indII=1
22 77 maps_order=[indI,indII]
23 78  
24   -temp_on=40. ;K
  79 +;; Loop through the maps corresponding to different filters
  80 +;; Generate some fake data using a BB function
  81 +;; ON pixels will be set to blackbody emission at 40K
  82 +;; OFF pixels will be set to blackbody emission at 20K
  83 +;; Loop through the maps corresponding to different filters
  84 +temp_on=40. ;Kelvin
25 85 temp_off=20.
26 86 FOR i=0L,Ndata_sets-1 DO BEGIN
27   - ;fill in intensity data
28 87 wave=dustem_filter2wav(filters[i])
29 88 this_map0=maps[*,*,indI,i]
30 89 this_map=this_map0
31 90 this_map[index_on]=dustem_planck_function(temp_on,wave)
32 91 this_map[index_off]=dustem_planck_function(temp_off,wave)
33 92 maps[*,*,indI,i]=this_map
34   - ;fill in intensity uncertainties
  93 +;; assign a variance value to each valid pixel
35 94 maps[*,*,indII,i]=la_power(la_mul(maps[*,*,indI,i],0.3),2.)
36 95 ENDFOR
37 96  
  97 +;; extract the default type of Stokes I only SED for the ON pixels (method = 'MEAN')
38 98 sed=dustem_sed_extractor(maps,index_on,filters,maps_order=maps_order,/total_intensity_only)
39 99  
40   -nsigma=10
  100 +;; show SED with 3-sigma error bars
  101 +nsigma=3
41 102 cgplot,sed.wave,sed.stokesI,/xlog,/ylog,xtit='Wavelength [mic]',ytit='SED [MJy/sr]'
42 103 cgoplot,sed.wave,sed.stokesI,psym='Filled Circle',color='red'
43 104 cgerrplot,sed.wave,sed.stokesI-nsigma*sqrt(sed.sigmaII),sed.stokesI+nsigma*sqrt(sed.sigmaII)
44 105  
  106 +;; extract a MEAN_ONOFF Stokes I only SED for the ON pixels (method = 'MEAN_ONOFF')
45 107 sed2=dustem_sed_extractor(maps,index_on,filters,maps_order=maps_order,/total_intensity_only,index_off=index_off,method='MEAN_ONOFF')
46 108  
47 109 cgoplot,sed2.wave,sed2.stokesI
48 110 cgoplot,sed2.wave,sed2.stokesI,psym='Filled Circle',color='blue'
49 111  
50   -stop
  112 +;; Try a simple BB fit of a dust model to SED
  113 +
  114 +;=== INFORMATION TO RUN THE FITS
  115 +tol=1.e-16 ;fit tolerence
  116 +use_Nitermax=10
  117 +
  118 +;=== INFORMATION TO MAKE THE PLOTS
  119 +yr=[1.00e1,8.00E8] ; y-axis limits
  120 +xr=[5.00E0,2.00e4] ; x-axis limits
  121 +tit='DUSTEM SED EXTRACT EXAMPLE' ; plot title
  122 +ytit=textoidl('I_\nu (MJy/sr)') ; y-axis title
  123 +xtit=textoidl('\lambda (\mum)') ; x-axis title
  124 +
  125 +pd = [ $
  126 + '(*!dustem_params).G0', $ ;G0
  127 + 'dustem_plugin_continuum_1', $ ;Temperature of a BB
  128 + 'dustem_plugin_continuum_2'] ;Peak amplitude of a BB
  129 +
  130 +iv = [1., 10., 1.e6] ; initial guess
  131 +
  132 +Npar=n_elements(pd)
  133 +ulimed=replicate(0,Npar)
  134 +llimed=replicate(1,Npar)
  135 +llims=replicate(1.e-15,Npar)
  136 +
  137 +
  138 +fpd = [ $
  139 + '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH
  140 + '(*!dustem_params).grains(1).mdust_o_mh',$ ;VSG
  141 + '(*!dustem_params).grains(2).mdust_o_mh'] ;BG
  142 +
  143 +fiv = [1.e-12,1.e-12,1.e-12]
  144 +
  145 +ind=where(sed.sigmaII LT (0.2*sed.StokesI)^2,count)
  146 +IF count NE 0 THEN sed[ind].sigmaII=(0.2*sed[ind].StokesI)^2
  147 +dustem_set_data,m_fit=sed,m_show=sed
  148 +dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims
  149 +
  150 +res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $
  151 + ,/xlog,/ylog,xr=xr,yr=yr,xtit=xtit,ytit=ytit,title=tit $
  152 + ,legend_xpos=legend_xpos,legend_ypos=legend_ypos $
  153 + ,errors=errors,chi2=chi2,rchi2=rchi2,show_plot=show_plot)
51 154  
52 155 the_end:
53 156  
54   -END
55 157 \ No newline at end of file
  158 +END
... ...
src/idl/dustem_sed_extractor.pro
  Diff comments   View file @ecfb8fb
1   -FUNCTION dustem_sed_extractor,maps,index_roi,filters, $
2   - comments=comments, $
3   - maps_order=maps_order, $
4   - method=method, $
5   - index_off=index_off, $
6   - reference_filter=reference_filter, $
7   - mask_I_undef=mask_I_undef, $
8   - total_intensity_only=total_intensity_only, $
9   - _extra=_extra
  1 +FUNCTION dustem_sed_extractor,maps $
  2 + ,index_roi $
  3 + ,filters, $
  4 + comments=comments, $
  5 + maps_order=maps_order, $
  6 + method=method, $
  7 + index_off=index_off, $
  8 + reference_filter=reference_filter, $
  9 + mask_I_undef=mask_I_undef, $
  10 + total_intensity_only=total_intensity_only, $
  11 + _extra=_extra
10 12  
11 13 ;+
12 14 ; NAME:
13 15 ; sed_extractor
  16 +;
14 17 ; PURPOSE:
15 18 ; Extract SEDs from a set of maps
  19 +;
16 20 ; CATEGORY:
17   -; Dustemwrap
  21 +; DustEMWrap, Distributed, High-Level, Analysis Helper
  22 +;
18 23 ; CALLING SEQUENCE:
19 24 ; sed=sed_extractor(maps,index,filters[,comments=][,maps_order=][,method=][,mask_I_undef=][,/total_intensity_only])
  25 +;
20 26 ; INPUTS:
21   -; maps : a set of Nmaps maps of same dimension [Nx,Ny,Nmaps,Ndata_sets]
  27 +; maps : a set of Nmaps maps of same dimension [Nx,Ny,Nmaps,Ndata_sets]
22 28 ; index_roi : spatial index in Nx,Ny where the SED is to be extracted from (ROI=region of interest)
23   -; filters : dustemwrap filters corresponding to each maps [Nmaps]
  29 +; filters : DustEMWrap filters corresponding to each dataset [Ndata_sets]
  30 +;
24 31 ; OPTIONAL INPUT PARAMETERS:
25 32 ; comments = a vector of strings to be attached to the SED. Also an output.
26 33 ; maps_order = a vector of 9 integers giving the 0-based index for I,Q,U,II,QQ,UU,IQ,IU,QU in that order
... ... @@ -28,27 +35,40 @@ FUNCTION dustem_sed_extractor,maps,index_roi,filters, $
28 35 ; method = extraction method (default = 'MEAN')/ Possible values : 'MEAN','MEAN_ONOFF','QU_CORREL'
29 36 ; index_off = same as index for the OFF region (needed only for method='MEAN_ONOFF')
30 37 ; reference_filter = reference filter to be correlated with each map (needed only for method='QU_CORREL')
  38 +;
31 39 ; OUTPUTS:
32 40 ; SED : a dustemwrap SED structure containing te extracted SED
  41 +;
33 42 ; OPTIONAL OUTPUT PARAMETERS:
34 43 ; mask_I_undef = a map of saturated pixels found in Stokes I
35   -; comments = a vector of strings to be writtena s comments into the SED.
  44 +; comments = a vector of strings to be written as comments into the SED.
  45 +;
36 46 ; ACCEPTED KEY-WORDS:
37 47 ; total_intensity_only = if set, data is assumed to contain only total intensity (no polarization)
38 48 ; comments = comments to be added to the sed. sed_extractor adds comments to this
  49 +;
39 50 ; COMMON BLOCKS:
40 51 ; None
  52 +;
41 53 ; SIDE EFFECTS:
42 54 ; None
  55 +;
43 56 ; RESTRICTIONS:
44 57 ; None
  58 +;
45 59 ; PROCEDURE:
46 60 ; None
  61 +;
47 62 ; EXAMPLES
  63 +;
48 64 ; MODIFICATION HISTORY:
49   -; Written by JPB
  65 +; Written by JPB 2022
  66 +; Evolution details on the DustEMWrap gitlab.
  67 +; See http://dustemwrap.irap.omp.eu/ for FAQ and help.
50 68 ;-
51 69  
  70 +
  71 +
52 72 IF keyword_set(help) THEN BEGIN
53 73 doc_library,'dustem_sed_extractor'
54 74 sed=0.
... ... @@ -126,73 +146,73 @@ ENDIF
126 146 CASE use_method OF
127 147 'MEAN':BEGIN
128 148 FOR j=0L,Nfilt-1 DO BEGIN
129   - this_map=used_maps[*,*,indI,j]
130   - sed[j].stokesI=la_mean(this_map[index_roi])
131   - this_map=used_maps[*,*,indII,j]
132   - sed[j].sigmaII=la_div(la_mean(this_map[index_roi]),count_roi)
133   - IF not keyword_set(total_intensity_only) THEN BEGIN
134   - this_map=used_maps[*,*,indQ,j]
135   - sed[j].stokesQ=la_mean(this_map[index_roi])
136   - this_map=used_maps[*,*,indU,j]
137   - sed[j].stokesU=la_mean(this_map[index_roi])
138   - this_map=used_maps[*,*,indQQ,j]
139   - sed[j].sigmaQQ=la_div(la_mean(this_map[index_roi]),count_roi)
140   - this_map=used_maps[*,*,indUU,j]
141   - sed[j].sigmaUU=la_div(la_mean(this_map[index_roi]),count_roi)
142   -
143   - this_map=used_maps[*,*,indIQ,j]
144   - sed[j].sigmaIQ=la_div(la_mean(this_map[index_roi]),count_roi)
145   - this_map=used_maps[*,*,indIU,j]
146   - sed[j].sigmaIU=la_div(la_mean(this_map[index_roi]),count_roi)
147   - this_map=used_maps[*,*,indQU,j]
148   - sed[j].sigmaQU=la_div(la_mean(this_map[index_roi]),count_roi)
149   - ENDIF
150   - ENDFOR
151   - END
152   - 'MEAN_ONOFF': BEGIN ;This is mean ON - mean MEAN_ONOFF
153   - IF not keyword_set(index_off) THEN BEGIN
154   - message,'index_off keyword must be set with method '+use_method,/continue
155   - stop
156   - ENDIF
157   - count_off=n_elements(index_off)
158   - comments=[comments,'OFF has '+strtrim(count_off,2)+' map pixels']
159   - FOR j=0L,Nfilt-1 DO BEGIN
160   - this_map=used_maps[*,*,indI,j]
161   - sed[j].stokesI=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off]))
162   - this_map=used_maps[*,*,indII,j]
163   - sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
164   - sig_off=la_div(la_mean(this_map[index_off]),count_off)
165   - sed[j].sigmaII=la_add(sig_on,sig_off)
166   - IF not keyword_set(total_intensity_only) THEN BEGIN
167   - this_map=used_maps[*,*,indQ,j]
168   - sed[j].stokesQ=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off]))
169   - this_map=used_maps[*,*,indU,j]
170   - sed[j].stokesU=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off]))
171   - ;variances
172   - this_map=used_maps[*,*,indQQ,j]
173   - sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
174   - sig_off=la_div(la_mean(this_map[index_off]),count_off)
175   - sed[j].sigmaQQ=la_add(sig_on,sig_off)
176   - this_map=used_maps[*,*,indUU,j]
177   - sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
178   - sig_off=la_div(la_mean(this_map[index_off]),count_off)
179   - sed[j].sigmaUU=la_add(sig_on,sig_off)
180   - ;co-variances
181   - this_map=used_maps[*,*,indIQ,j]
182   - sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
183   - sig_off=la_div(la_mean(this_map[index_off]),count_off)
184   - sed[j].sigmaIQ=la_add(sig_on,sig_off)
185   - this_map=used_maps[*,*,indIU,j]
186   - sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
187   - sig_off=la_div(la_mean(this_map[index_off]),count_off)
188   - sed[j].sigmaIU=la_add(sig_on,sig_off)
189   - this_map=used_maps[*,*,indQU,j]
190   - sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
191   - sig_off=la_div(la_mean(this_map[index_off]),count_off)
192   - sed[j].sigmaQU=la_add(sig_on,sig_off)
193   - ENDIF
194   - ENDFOR
195   - ;stop
  149 + this_map=used_maps[*,*,indI,j]
  150 + sed[j].stokesI=la_mean(this_map[index_roi])
  151 + this_map=used_maps[*,*,indII,j]
  152 + sed[j].sigmaII=la_div(la_mean(this_map[index_roi]),count_roi)
  153 + IF not keyword_set(total_intensity_only) THEN BEGIN
  154 + this_map=used_maps[*,*,indQ,j]
  155 + sed[j].stokesQ=la_mean(this_map[index_roi])
  156 + this_map=used_maps[*,*,indU,j]
  157 + sed[j].stokesU=la_mean(this_map[index_roi])
  158 + this_map=used_maps[*,*,indQQ,j]
  159 + sed[j].sigmaQQ=la_div(la_mean(this_map[index_roi]),count_roi)
  160 + this_map=used_maps[*,*,indUU,j]
  161 + sed[j].sigmaUU=la_div(la_mean(this_map[index_roi]),count_roi)
  162 +
  163 + this_map=used_maps[*,*,indIQ,j]
  164 + sed[j].sigmaIQ=la_div(la_mean(this_map[index_roi]),count_roi)
  165 + this_map=used_maps[*,*,indIU,j]
  166 + sed[j].sigmaIU=la_div(la_mean(this_map[index_roi]),count_roi)
  167 + this_map=used_maps[*,*,indQU,j]
  168 + sed[j].sigmaQU=la_div(la_mean(this_map[index_roi]),count_roi)
  169 + ENDIF
  170 + ENDFOR
  171 + END
  172 + 'MEAN_ONOFF': BEGIN ;This is mean ON - mean MEAN_ONOFF
  173 + IF not keyword_set(index_off) THEN BEGIN
  174 + message,'index_off keyword must be set with method '+use_method,/continue
  175 + stop
  176 + ENDIF
  177 + count_off=n_elements(index_off)
  178 + comments=[comments,'OFF has '+strtrim(count_off,2)+' map pixels']
  179 + FOR j=0L,Nfilt-1 DO BEGIN
  180 + this_map=used_maps[*,*,indI,j]
  181 + sed[j].stokesI=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off]))
  182 + this_map=used_maps[*,*,indII,j]
  183 + sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
  184 + sig_off=la_div(la_mean(this_map[index_off]),count_off)
  185 + sed[j].sigmaII=la_add(sig_on,sig_off)
  186 + IF not keyword_set(total_intensity_only) THEN BEGIN
  187 + this_map=used_maps[*,*,indQ,j]
  188 + sed[j].stokesQ=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off]))
  189 + this_map=used_maps[*,*,indU,j]
  190 + sed[j].stokesU=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off]))
  191 + ;variances
  192 + this_map=used_maps[*,*,indQQ,j]
  193 + sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
  194 + sig_off=la_div(la_mean(this_map[index_off]),count_off)
  195 + sed[j].sigmaQQ=la_add(sig_on,sig_off)
  196 + this_map=used_maps[*,*,indUU,j]
  197 + sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
  198 + sig_off=la_div(la_mean(this_map[index_off]),count_off)
  199 + sed[j].sigmaUU=la_add(sig_on,sig_off)
  200 + ;co-variances
  201 + this_map=used_maps[*,*,indIQ,j]
  202 + 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 + sed[j].sigmaIQ=la_add(sig_on,sig_off)
  205 + this_map=used_maps[*,*,indIU,j]
  206 + sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
  207 + sig_off=la_div(la_mean(this_map[index_off]),count_off)
  208 + sed[j].sigmaIU=la_add(sig_on,sig_off)
  209 + this_map=used_maps[*,*,indQU,j]
  210 + sig_on=la_div(la_mean(this_map[index_roi]),count_roi)
  211 + sig_off=la_div(la_mean(this_map[index_off]),count_off)
  212 + sed[j].sigmaQU=la_add(sig_on,sig_off)
  213 + ENDIF
  214 + ENDFOR
  215 + ;stop
196 216 END
197 217 'QU_CORREL': BEGIN ;This is to compute SED of polarization correlated to a reference
198 218 IF keyword_set(total_intensity_only) THEN BEGIN
... ... @@ -298,4 +318,4 @@ sed=ssed
298 318 the_end:
299 319 RETURN,sed
300 320  
301   -END
302 321 \ No newline at end of file
  322 +END
... ...