OV-GSO-DC / dustem-wrapper_idl

Commit 5c381fd96d8427dc28f60b67905677c67f711a98

Authored by Jean-Philippe Bernard
1 parent 7dfd7844
Exists in master

first commit

Showing 1 changed file with 301 additions and 0 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_sed_extractor.pro 0 → 100644
  Diff comments   View file @5c381fd
... ... @@ -0,0 +1,301 @@
  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
  10 +
  11 +;+
  12 +; NAME:
  13 +; sed_extractor
  14 +; PURPOSE:
  15 +; Extract SEDs from a set of maps
  16 +; CATEGORY:
  17 +; Dustemwrap
  18 +; CALLING SEQUENCE:
  19 +; sed=sed_extractor(maps,index,filters[,comments=][,maps_order=][,method=][,mask_I_undef=][,/total_intensity_only])
  20 +; INPUTS:
  21 +; maps : a set of Nmaps maps of same dimension [Nx,Ny,Nmaps,Ndata_sets]
  22 +; 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]
  24 +; OPTIONAL INPUT PARAMETERS:
  25 +; comments = a vector of strings to be attached to the SED. Also an output.
  26 +; 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
  27 +; (default [0,1,2,3,4,5,6,7,8]). For total intensity only, use something like [0,-1,-1,1,-1,-1,-1,-1,-1].
  28 +; method = extraction method (default = 'MEAN')/ Possible values : 'MEAN','MEAN_ONOFF','QU_CORREL'
  29 +; index_off = same as index for the OFF region (needed only for method='MEAN_ONOFF')
  30 +; reference_filter = reference filter to be correlated with each map (needed only for method='QU_CORREL')
  31 +; OUTPUTS:
  32 +; SED : a dustemwrap SED structure containing te extracted SED
  33 +; OPTIONAL OUTPUT PARAMETERS:
  34 +; 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.
  36 +; ACCEPTED KEY-WORDS:
  37 +; total_intensity_only = if set, data is assumed to contain only total intensity (no polarization)
  38 +; comments = comments to be added to the sed. sed_extractor adds comments to this
  39 +; COMMON BLOCKS:
  40 +; None
  41 +; SIDE EFFECTS:
  42 +; None
  43 +; RESTRICTIONS:
  44 +; None
  45 +; PROCEDURE:
  46 +; None
  47 +; EXAMPLES
  48 +; MODIFICATION HISTORY:
  49 +; Written by JPB
  50 +;-
  51 +
  52 +IF keyword_set(help) THEN BEGIN
  53 + doc_library,'dustem_sed_extractor'
  54 + sed=0.
  55 + goto,the_end
  56 +ENDIF
  57 +
  58 +Nfilt=n_elements(filters)
  59 +sed=dustem_initialize_sed(Nfilt,comments=comments)
  60 +sed.filter=filters
  61 +sed.instru=dustem_filter2instru(sed.filter)
  62 +sed.wave=dustem_filter2wav(filters)
  63 +comments=[comments,'\Computed using sed_extractor on '+systime(0)]
  64 +order=sort(sed.wave) ;wavelength order
  65 +
  66 +count_roi=n_elements(index_roi)
  67 +
  68 +use_method='MEAN' ;This is the default extraction method
  69 +IF keyword_set(method) THEN BEGIN
  70 + use_method=method
  71 +ENDIF
  72 +comments=[comments,'SED extracted using method '+use_method]
  73 +comments=[comments,'ROI has '+strtrim(count_roi,2)+' map pixels']
  74 +N1=(size(maps))[1]
  75 +N2=(size(maps))[2]
  76 +Nplanes=(size(maps))[3]
  77 +
  78 +IF keyword_set(total_intensity_only) THEN BEGIN
  79 + indI=0L & indII=1L
  80 + IF keyword_set(maps_order) THEN BEGIN
  81 + indI=maps_order[0]
  82 + indII=maps_order[1]
  83 + ENDIF
  84 +ENDIF ELSE BEGIN
  85 + indI=0L & indQ=1L & indU=2L & indII=3L & indQQ=4L & indUU=5L &indIQ=6L & indIU=7L &indQU=8L
  86 + IF keyword_set(maps_order) THEN BEGIN
  87 + indI=maps_order[0]
  88 + indQ=maps_order[1]
  89 + indU=maps_order[2]
  90 + indII=maps_order[3]
  91 + indQQ=maps_order[4]
  92 + indUU=maps_order[5]
  93 + indIQ=maps_order[6]
  94 + indIU=maps_order[7]
  95 + indQU=maps_order[8]
  96 + ENDIF
  97 +ENDELSE
  98 +
  99 +;===== do a map of saturated stokesI pixels
  100 +mask_I_undef=maps[*,*,0,0]
  101 +mask_I_undef[*]=0
  102 +FOR j=0L,Nfilt-1 DO BEGIN
  103 + this_map=maps[*,*,indI,j]
  104 + ind=where(this_map EQ la_undef(),count)
  105 + IF count NE 0 AND count NE N1*N2 THEN BEGIN
  106 + message,'Found '+strtrim(count,2)+' undefined values for map of '+filters[j],/continue
  107 + ;stop
  108 + mask_I_undef[ind]=1
  109 + ENDIF
  110 +ENDFOR
  111 +
  112 +used_maps=maps
  113 +;set saturated stokesI pixels to la_undef in
  114 +indbad=where(mask_I_undef EQ 1,count)
  115 +IF count NE 0 THEN BEGIN
  116 + message,'Setting Stokes I saturated values to undef in the whole cube',/continue
  117 + FOR j=0L,Nfilt-1 DO BEGIN
  118 + FOR k=0L,Nplanes-1 DO BEGIN
  119 + this_map=used_maps[*,*,k,j]
  120 + this_map[indbad]=la_undef()
  121 + used_maps[*,*,k,j]=this_map
  122 + ENDFOR
  123 + ENDFOR
  124 +ENDIF
  125 +
  126 +CASE use_method OF
  127 + 'MEAN':BEGIN
  128 + 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
  196 + END
  197 + 'QU_CORREL': BEGIN ;This is to compute SED of polarization correlated to a reference
  198 + IF keyword_set(total_intensity_only) THEN BEGIN
  199 + message,'mode QU_CORREL can only be used with polarization data'
  200 + ENDIF
  201 + IF not keyword_set(reference_filter) THEN BEGIN
  202 + message,'reference_filter keyword must be set with method '+use_method,/continue
  203 + stop
  204 + ENDIF
  205 + indf=where(filters EQ reference_filter,countf)
  206 + IF countf EQ 0 THEN BEGIN
  207 + message,'reference filter '+reference_filter+' not found in provided filter list',/continue
  208 + stop
  209 + ENDIF
  210 + mapI_ref=used_maps[*,*,indI,indf]
  211 + mapQ_ref=used_maps[*,*,indQ,indf]
  212 + mapU_ref=used_maps[*,*,indU,indf]
  213 + mapII_ref=used_maps[*,*,indII,indf]
  214 + mapQQ_ref=used_maps[*,*,indQQ,indf]
  215 + mapUU_ref=used_maps[*,*,indUU,indf]
  216 + Iref=la_mean(mapI_ref[index_roi])
  217 + Qref=la_mean(mapQ_ref[index_roi])
  218 + Uref=la_mean(mapU_ref[index_roi])
  219 + FOR j=0L,Nfilt-1 DO BEGIN
  220 + ;==== do the correlation in Stokes I
  221 + x=mapI_ref[index_roi]
  222 + y=(used_maps[*,*,indI,j])[index_roi]
  223 + sx=la_power(mapII_ref[index_roi],0.5)
  224 + sy=la_power((used_maps[*,*,indI,j])[index_roi],0.5)
  225 + ind_good=where(x NE la_undef() AND y NE la_undef() AND sx NE la_undef() and sy NE la_undef(),count_good)
  226 + comments=[comments,'I corelation has '+strtrim(count_good,2)+' map pixels']
  227 + IF count_good NE 0 THEN BEGIN
  228 + x=x[ind_good]
  229 + y=y[ind_good]
  230 + sx=sx[ind_good]
  231 + sy=sy[ind_good]
  232 + start=[1.,1.]
  233 + try=linear_mpfit(x,y,sx,sy,start,status=status $
  234 + ,perror=perror,bestnorm=bestnorm,_extra=_extra)
  235 + ;stop
  236 + ENDIF ELSE BEGIN
  237 + status=0
  238 + try=[la_undef(),la_undef()]
  239 + ENDELSE
  240 + IF status NE 0 AND status NE -16 THEN BEGIN
  241 + sed[j].stokesI=try[0]*Iref
  242 + sed[j].sigmaII=la_power(la_mul(perror[0],Iref),2.)
  243 + ENDIF
  244 + ;Do the correlation between Q maps
  245 + x=[mapQ_ref[index_roi],mapU_ref[index_roi]]
  246 + y=[(used_maps[*,*,indQ,j])[index_roi],(used_maps[*,*,indU,j])[index_roi]]
  247 + sx=la_power([mapQQ_ref[index_roi],mapUU_ref[index_roi]],0.5)
  248 + sy=la_power([(used_maps[*,*,indQ,j])[index_roi],(used_maps[*,*,indU,j])[index_roi]],0.5)
  249 + ind_good=where(x NE la_undef() AND y NE la_undef() AND sx NE la_undef() and sy NE la_undef(),count_good)
  250 + comments=[comments,'QU corelation has '+strtrim(count_good,2)+' map pixels']
  251 + IF count_good NE 0 THEN BEGIN
  252 + x=x[ind_good]
  253 + y=y[ind_good]
  254 + sx=sx[ind_good]
  255 + sy=sy[ind_good]
  256 + start=[1.,1.]
  257 + try=linear_mpfit(x,y,sx,sy,start,status=status $
  258 + ,perror=perror,bestnorm=bestnorm,_extra=_extra)
  259 + ;stop
  260 + ENDIF ELSE BEGIN
  261 + status=0
  262 + try=[la_undef(),la_undef()]
  263 + ENDELSE
  264 + IF status NE 0 AND status NE -16 THEN BEGIN
  265 + sed[j].stokesQ=la_mul(try[0],Qref)
  266 + sed[j].stokesU=la_mul(try[0],Uref)
  267 + sed[j].sigmaQQ=la_power(la_mul(perror[0],Qref),2.)
  268 + sed[j].sigmaUU=la_power(la_mul(perror[0],Uref),2.)
  269 + ;for now assume no co-variances. But we could aclculate this, right ?
  270 + sed[j].sigmaIQ=0.
  271 + sed[j].sigmaIU=0.
  272 + sed[j].sigmaQU=0.
  273 + ;IF filters[j] EQ reference_filter THEN stop
  274 + ENDIF
  275 + ENDFOR
  276 + comments=[comments,'Reference filter for QU correlation '+reference_filter]
  277 + END
  278 + ELSE: BEGIN
  279 + message,'Method '+use_method+' not implemented',/continue
  280 + sed=0L
  281 + goto,the_end
  282 + END
  283 +ENDCASE
  284 +
  285 +polar_iqu2ippsi,sed.stokesI,sed.stokesQ,sed.stokesU,p_values,psi_values,/lac,largep=largep
  286 +sed.smallp=p_values
  287 +sed.psi=psi_values
  288 +sed.largep=largep
  289 +polar_variance_iqu2ippsi,sed.stokesI,sed.stokesQ,sed.stokesU,sed.sigmaII,sed.sigmaQQ,sed.sigmaUU,sed.sigmaIQ,sed.sigmaIU,sed.sigmaQU,variances_p,variances_psi,variances_largep,/lac
  290 +sed.sigma_smallp=variances_p
  291 +sed.sigma_largep=variances_largep
  292 +sed.sigma_psi=variances_psi
  293 +ssed=sed[order]
  294 +sed=ssed
  295 +
  296 +;stop
  297 +
  298 +the_end:
  299 +RETURN,sed
  300 +
  301 +END
0 302 \ No newline at end of file
... ...