FUNCTION dustem_sed_extractor,maps $ ,index_roi $ ,filters, $ comments=comments, $ maps_order=maps_order, $ method=method, $ index_off=index_off, $ reference_filter=reference_filter, $ mask_I_undef=mask_I_undef, $ total_intensity_only=total_intensity_only, $ _extra=_extra ;+ ; NAME: ; sed_extractor ; ; PURPOSE: ; Extract SEDs from a set of maps ; ; CATEGORY: ; DustEMWrap, Distributed, High-Level, Analysis Helper ; ; CALLING SEQUENCE: ; sed=sed_extractor(maps,index,filters[,comments=][,maps_order=][,method=][,mask_I_undef=][,/total_intensity_only]) ; ; INPUTS: ; maps : a set of Nmaps maps of same dimension [Nx,Ny,Nmaps,Ndata_sets] ; index_roi : spatial index in Nx,Ny where the SED is to be extracted from (ROI=region of interest) ; filters : DustEMWrap filters corresponding to each dataset [Ndata_sets] ; ; OPTIONAL INPUT PARAMETERS: ; comments = a vector of strings to be attached to the SED. Also an output. ; 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 ; (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]. ; method = extraction method (default = 'MEAN')/ Possible values : 'MEAN','MEAN_ONOFF','QU_CORREL' ; index_off = same as index for the OFF region (needed only for method='MEAN_ONOFF') ; reference_filter = reference filter to be correlated with each map (needed only for method='QU_CORREL') ; ; OUTPUTS: ; SED : a dustemwrap SED structure containing te extracted SED ; ; OPTIONAL OUTPUT PARAMETERS: ; mask_I_undef = a map of saturated pixels found in Stokes I ; comments = a vector of strings to be written as comments into the SED. ; ; ACCEPTED KEY-WORDS: ; total_intensity_only = if set, data is assumed to contain only total intensity (no polarization) ; comments = comments to be added to the sed. sed_extractor adds comments to this ; ; COMMON BLOCKS: ; None ; ; SIDE EFFECTS: ; None ; ; RESTRICTIONS: ; None ; ; PROCEDURE: ; None ; ; EXAMPLES ; ; MODIFICATION HISTORY: ; Written by JPB 2022 ; Evolution details on the DustEMWrap gitlab. ; See http://dustemwrap.irap.omp.eu/ for FAQ and help. ;- IF keyword_set(help) THEN BEGIN doc_library,'dustem_sed_extractor' sed=0. goto,the_end ENDIF Nfilt=n_elements(filters) sed=dustem_initialize_sed(Nfilt,comments=comments) sed.filter=filters sed.instru=dustem_filter2instru(sed.filter) sed.wave=dustem_filter2wav(filters) comments=[comments,'\Computed using sed_extractor on '+systime(0)] order=sort(sed.wave) ;wavelength order count_roi=n_elements(index_roi) use_method='MEAN' ;This is the default extraction method IF keyword_set(method) THEN BEGIN use_method=method ENDIF comments=[comments,'SED extracted using method '+use_method] comments=[comments,'ROI has '+strtrim(count_roi,2)+' map pixels'] N1=(size(maps))[1] N2=(size(maps))[2] Nplanes=(size(maps))[3] IF keyword_set(total_intensity_only) THEN BEGIN indI=0L & indII=1L IF keyword_set(maps_order) THEN BEGIN indI=maps_order[0] indII=maps_order[1] ENDIF ENDIF ELSE BEGIN indI=0L & indQ=1L & indU=2L & indII=3L & indQQ=4L & indUU=5L &indIQ=6L & indIU=7L &indQU=8L IF keyword_set(maps_order) THEN BEGIN indI=maps_order[0] indQ=maps_order[1] indU=maps_order[2] indII=maps_order[3] indQQ=maps_order[4] indUU=maps_order[5] indIQ=maps_order[6] indIU=maps_order[7] indQU=maps_order[8] ENDIF ENDELSE ;===== do a map of saturated stokesI pixels mask_I_undef=maps[*,*,0,0] mask_I_undef[*]=0 FOR j=0L,Nfilt-1 DO BEGIN this_map=maps[*,*,indI,j] ind=where(this_map EQ la_undef(),count) IF count NE 0 AND count NE N1*N2 THEN BEGIN message,'Found '+strtrim(count,2)+' undefined values for map of '+filters[j],/continue ;stop mask_I_undef[ind]=1 ENDIF ENDFOR used_maps=maps ;set saturated stokesI pixels to la_undef in indbad=where(mask_I_undef EQ 1,count) IF count NE 0 THEN BEGIN message,'Setting Stokes I saturated values to undef in the whole cube',/continue FOR j=0L,Nfilt-1 DO BEGIN FOR k=0L,Nplanes-1 DO BEGIN this_map=used_maps[*,*,k,j] this_map[indbad]=la_undef() used_maps[*,*,k,j]=this_map ENDFOR ENDFOR ENDIF CASE use_method OF 'MEAN':BEGIN FOR j=0L,Nfilt-1 DO BEGIN this_map=used_maps[*,*,indI,j] sed[j].stokesI=la_mean(this_map[index_roi]) this_map=used_maps[*,*,indII,j] sed[j].sigmaII=la_div(la_mean(this_map[index_roi]),count_roi) IF not keyword_set(total_intensity_only) THEN BEGIN this_map=used_maps[*,*,indQ,j] sed[j].stokesQ=la_mean(this_map[index_roi]) this_map=used_maps[*,*,indU,j] sed[j].stokesU=la_mean(this_map[index_roi]) this_map=used_maps[*,*,indQQ,j] sed[j].sigmaQQ=la_div(la_mean(this_map[index_roi]),count_roi) this_map=used_maps[*,*,indUU,j] sed[j].sigmaUU=la_div(la_mean(this_map[index_roi]),count_roi) this_map=used_maps[*,*,indIQ,j] sed[j].sigmaIQ=la_div(la_mean(this_map[index_roi]),count_roi) this_map=used_maps[*,*,indIU,j] sed[j].sigmaIU=la_div(la_mean(this_map[index_roi]),count_roi) this_map=used_maps[*,*,indQU,j] sed[j].sigmaQU=la_div(la_mean(this_map[index_roi]),count_roi) ENDIF ENDFOR END 'MEAN_ONOFF': BEGIN ;This is mean ON - mean MEAN_ONOFF IF not keyword_set(index_off) THEN BEGIN message,'index_off keyword must be set with method '+use_method,/continue stop ENDIF count_off=n_elements(index_off) comments=[comments,'OFF has '+strtrim(count_off,2)+' map pixels'] FOR j=0L,Nfilt-1 DO BEGIN this_map=used_maps[*,*,indI,j] sed[j].stokesI=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off])) this_map=used_maps[*,*,indII,j] sig_on=la_div(la_mean(this_map[index_roi]),count_roi) sig_off=la_div(la_mean(this_map[index_off]),count_off) sed[j].sigmaII=la_add(sig_on,sig_off) IF not keyword_set(total_intensity_only) THEN BEGIN this_map=used_maps[*,*,indQ,j] sed[j].stokesQ=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off])) this_map=used_maps[*,*,indU,j] sed[j].stokesU=la_sub(la_mean(this_map[index_roi]),la_mean(this_map[index_off])) ;variances this_map=used_maps[*,*,indQQ,j] sig_on=la_div(la_mean(this_map[index_roi]),count_roi) sig_off=la_div(la_mean(this_map[index_off]),count_off) sed[j].sigmaQQ=la_add(sig_on,sig_off) this_map=used_maps[*,*,indUU,j] sig_on=la_div(la_mean(this_map[index_roi]),count_roi) sig_off=la_div(la_mean(this_map[index_off]),count_off) sed[j].sigmaUU=la_add(sig_on,sig_off) ;co-variances this_map=used_maps[*,*,indIQ,j] sig_on=la_div(la_mean(this_map[index_roi]),count_roi) sig_off=la_div(la_mean(this_map[index_off]),count_off) sed[j].sigmaIQ=la_add(sig_on,sig_off) this_map=used_maps[*,*,indIU,j] sig_on=la_div(la_mean(this_map[index_roi]),count_roi) sig_off=la_div(la_mean(this_map[index_off]),count_off) sed[j].sigmaIU=la_add(sig_on,sig_off) this_map=used_maps[*,*,indQU,j] sig_on=la_div(la_mean(this_map[index_roi]),count_roi) sig_off=la_div(la_mean(this_map[index_off]),count_off) sed[j].sigmaQU=la_add(sig_on,sig_off) ENDIF ENDFOR ;stop END 'QU_CORREL': BEGIN ;This is to compute SED of polarization correlated to a reference IF keyword_set(total_intensity_only) THEN BEGIN message,'mode QU_CORREL can only be used with polarization data' ENDIF IF not keyword_set(reference_filter) THEN BEGIN message,'reference_filter keyword must be set with method '+use_method,/continue stop ENDIF indf=where(filters EQ reference_filter,countf) IF countf EQ 0 THEN BEGIN message,'reference filter '+reference_filter+' not found in provided filter list',/continue stop ENDIF mapI_ref=used_maps[*,*,indI,indf] mapQ_ref=used_maps[*,*,indQ,indf] mapU_ref=used_maps[*,*,indU,indf] mapII_ref=used_maps[*,*,indII,indf] mapQQ_ref=used_maps[*,*,indQQ,indf] mapUU_ref=used_maps[*,*,indUU,indf] Iref=la_mean(mapI_ref[index_roi]) Qref=la_mean(mapQ_ref[index_roi]) Uref=la_mean(mapU_ref[index_roi]) FOR j=0L,Nfilt-1 DO BEGIN ;==== do the correlation in Stokes I x=mapI_ref[index_roi] y=(used_maps[*,*,indI,j])[index_roi] sx=la_power(mapII_ref[index_roi],0.5) sy=la_power((used_maps[*,*,indI,j])[index_roi],0.5) ind_good=where(x NE la_undef() AND y NE la_undef() AND sx NE la_undef() and sy NE la_undef(),count_good) comments=[comments,'I corelation has '+strtrim(count_good,2)+' map pixels'] IF count_good NE 0 THEN BEGIN x=x[ind_good] y=y[ind_good] sx=sx[ind_good] sy=sy[ind_good] start=[1.,1.] try=linear_mpfit(x,y,sx,sy,start,status=status $ ,perror=perror,bestnorm=bestnorm,_extra=_extra) ;stop ENDIF ELSE BEGIN status=0 try=[la_undef(),la_undef()] ENDELSE IF status NE 0 AND status NE -16 THEN BEGIN sed[j].stokesI=try[0]*Iref sed[j].sigmaII=la_power(la_mul(perror[0],Iref),2.) ENDIF ;Do the correlation between Q maps x=[mapQ_ref[index_roi],mapU_ref[index_roi]] y=[(used_maps[*,*,indQ,j])[index_roi],(used_maps[*,*,indU,j])[index_roi]] sx=la_power([mapQQ_ref[index_roi],mapUU_ref[index_roi]],0.5) sy=la_power([(used_maps[*,*,indQ,j])[index_roi],(used_maps[*,*,indU,j])[index_roi]],0.5) ind_good=where(x NE la_undef() AND y NE la_undef() AND sx NE la_undef() and sy NE la_undef(),count_good) comments=[comments,'QU corelation has '+strtrim(count_good,2)+' map pixels'] IF count_good NE 0 THEN BEGIN x=x[ind_good] y=y[ind_good] sx=sx[ind_good] sy=sy[ind_good] start=[1.,1.] try=linear_mpfit(x,y,sx,sy,start,status=status $ ,perror=perror,bestnorm=bestnorm,_extra=_extra) ;stop ENDIF ELSE BEGIN status=0 try=[la_undef(),la_undef()] ENDELSE IF status NE 0 AND status NE -16 THEN BEGIN sed[j].stokesQ=la_mul(try[0],Qref) sed[j].stokesU=la_mul(try[0],Uref) sed[j].sigmaQQ=la_power(la_mul(perror[0],Qref),2.) sed[j].sigmaUU=la_power(la_mul(perror[0],Uref),2.) ;for now assume no co-variances. But we could aclculate this, right ? sed[j].sigmaIQ=0. sed[j].sigmaIU=0. sed[j].sigmaQU=0. ;IF filters[j] EQ reference_filter THEN stop ENDIF ENDFOR comments=[comments,'Reference filter for QU correlation '+reference_filter] END ELSE: BEGIN message,'Method '+use_method+' not implemented',/continue sed=0L goto,the_end END ENDCASE polar_iqu2ippsi,sed.stokesI,sed.stokesQ,sed.stokesU,p_values,psi_values,/lac,largep=largep sed.smallp=p_values sed.psi=psi_values sed.largep=largep 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 sed.sigma_smallp=variances_p sed.sigma_largep=variances_largep sed.sigma_psi=variances_psi ssed=sed[order] sed=ssed ;stop the_end: RETURN,sed END