diff --git a/src/idl/dustem_extract_sed_example.pro b/src/idl/dustem_extract_sed_example.pro index c5f74ca..3503dc4 100644 --- a/src/idl/dustem_extract_sed_example.pro +++ b/src/idl/dustem_extract_sed_example.pro @@ -1,18 +1,73 @@ -PRO test_dustem_sed_extractor +PRO dustem_extract_sed_example + +;+ +; NAME: +; dustem_extract_sed_example +; +; PURPOSE: +; This routine contains a simple example of how to extract an SED to use with DustEMWrap. +; +; CATEGORY: +; DustEMWrap, Distributed, High-Level, User Example +; +; CALLING SEQUENCE: +; dustem_extract_sed_example +; +; INPUTS: +; None +; +; OPTIONAL INPUT PARAMETERS: +; None +; +; OUTPUTS: +; None +; +; OPTIONAL OUTPUT PARAMETERS: +; +; ACCEPTED KEY-WORDS: +; None +; +; COMMON BLOCKS: +; None +; +; SIDE EFFECTS: +; None +; +; RESTRICTIONS: +; The DustEM fortran code must be installed +; The DustEMWrap IDL code must be installed +; +; PROCEDURES AND SUBROUTINES USED: +; +; EXAMPLES +; +; MODIFICATION HISTORY: +; Written by JPB Jan 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_extract_sed_example' + goto,the_end +END define_la_common -dustem_init +use_model='DBP90' +use_polarization=0 +dustem_init,model=use_model,polarization=use_polarization -Nmaps=2 ;This is total intensity and variances +Nmaptypes=2 ; number of map types in the structure. Here 2 (Stokes I, variance in I) filters=['IRAS3','IRAS4','PACS1','PACS2','PACS3','PILOT1','SPIRE1','SPIRE2','SPIRE3','HFI3','HFI4','HFI5','HFI6'] Ndata_sets=n_elements(filters) -maps=fltarr(100,100,Nmaps,Ndata_sets) +maps=fltarr(100,100,Nmaptypes,Ndata_sets) maps[*]=la_undef() +;; Define the region of interest (ON) and OFF region mask=maps[*,*,0,0] -mask[5:35,5:35]=1 -mask[10:30,10:30]=mask[10:30,10:30]+1 +mask[5:35,5:35]=1 ; OFF will be a guard band of 5 pixel around the ON defined below +mask[10:30,10:30]=mask[10:30,10:30]+1 ; ON source index_on=where(mask EQ 2,count_on) index_off=where(mask EQ 1,count_off) @@ -21,34 +76,82 @@ indI=0 indII=1 maps_order=[indI,indII] -temp_on=40. ;K +;; Loop through the maps corresponding to different filters +;; Generate some fake data using a BB function +;; ON pixels will be set to blackbody emission at 40K +;; OFF pixels will be set to blackbody emission at 20K +;; Loop through the maps corresponding to different filters +temp_on=40. ;Kelvin temp_off=20. FOR i=0L,Ndata_sets-1 DO BEGIN - ;fill in intensity data wave=dustem_filter2wav(filters[i]) this_map0=maps[*,*,indI,i] this_map=this_map0 this_map[index_on]=dustem_planck_function(temp_on,wave) this_map[index_off]=dustem_planck_function(temp_off,wave) maps[*,*,indI,i]=this_map - ;fill in intensity uncertainties +;; assign a variance value to each valid pixel maps[*,*,indII,i]=la_power(la_mul(maps[*,*,indI,i],0.3),2.) ENDFOR +;; extract the default type of Stokes I only SED for the ON pixels (method = 'MEAN') sed=dustem_sed_extractor(maps,index_on,filters,maps_order=maps_order,/total_intensity_only) -nsigma=10 +;; show SED with 3-sigma error bars +nsigma=3 cgplot,sed.wave,sed.stokesI,/xlog,/ylog,xtit='Wavelength [mic]',ytit='SED [MJy/sr]' cgoplot,sed.wave,sed.stokesI,psym='Filled Circle',color='red' cgerrplot,sed.wave,sed.stokesI-nsigma*sqrt(sed.sigmaII),sed.stokesI+nsigma*sqrt(sed.sigmaII) +;; extract a MEAN_ONOFF Stokes I only SED for the ON pixels (method = 'MEAN_ONOFF') sed2=dustem_sed_extractor(maps,index_on,filters,maps_order=maps_order,/total_intensity_only,index_off=index_off,method='MEAN_ONOFF') cgoplot,sed2.wave,sed2.stokesI cgoplot,sed2.wave,sed2.stokesI,psym='Filled Circle',color='blue' -stop +;; Try a simple BB fit of a dust model to SED + +;=== INFORMATION TO RUN THE FITS +tol=1.e-16 ;fit tolerence +use_Nitermax=10 + +;=== INFORMATION TO MAKE THE PLOTS +yr=[1.00e1,8.00E8] ; y-axis limits +xr=[5.00E0,2.00e4] ; x-axis limits +tit='DUSTEM SED EXTRACT EXAMPLE' ; plot title +ytit=textoidl('I_\nu (MJy/sr)') ; y-axis title +xtit=textoidl('\lambda (\mum)') ; x-axis title + +pd = [ $ + '(*!dustem_params).G0', $ ;G0 + 'dustem_plugin_continuum_1', $ ;Temperature of a BB + 'dustem_plugin_continuum_2'] ;Peak amplitude of a BB + +iv = [1., 10., 1.e6] ; initial guess + +Npar=n_elements(pd) +ulimed=replicate(0,Npar) +llimed=replicate(1,Npar) +llims=replicate(1.e-15,Npar) + + +fpd = [ $ + '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH + '(*!dustem_params).grains(1).mdust_o_mh',$ ;VSG + '(*!dustem_params).grains(2).mdust_o_mh'] ;BG + +fiv = [1.e-12,1.e-12,1.e-12] + +ind=where(sed.sigmaII LT (0.2*sed.StokesI)^2,count) +IF count NE 0 THEN sed[ind].sigmaII=(0.2*sed[ind].StokesI)^2 +dustem_set_data,m_fit=sed,m_show=sed +dustem_init_params,use_model,pd,iv,fpd=fpd,fiv=fiv,ulimed=ulimed,llimed=llimed,ulims=ulims,llims=llims + +res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol $ + ,/xlog,/ylog,xr=xr,yr=yr,xtit=xtit,ytit=ytit,title=tit $ + ,legend_xpos=legend_xpos,legend_ypos=legend_ypos $ + ,errors=errors,chi2=chi2,rchi2=rchi2,show_plot=show_plot) the_end: -END \ No newline at end of file +END diff --git a/src/idl/dustem_sed_extractor.pro b/src/idl/dustem_sed_extractor.pro index 388d8f7..063bae4 100644 --- a/src/idl/dustem_sed_extractor.pro +++ b/src/idl/dustem_sed_extractor.pro @@ -1,26 +1,33 @@ -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 +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 +; 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] +; 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 maps [Nmaps] +; 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 @@ -28,27 +35,40 @@ FUNCTION dustem_sed_extractor,maps,index_roi,filters, $ ; 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 writtena s comments into the SED. +; 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 +; 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. @@ -126,73 +146,73 @@ 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 + 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 @@ -298,4 +318,4 @@ sed=ssed the_end: RETURN,sed -END \ No newline at end of file +END -- libgit2 0.21.2