Merge branch 'master' of https://gitlab.irap.omp.eu/OV-GSO-DC/dustem-wrapper_idl

Jean-Philippe Bernard
2 parents 66d0acb9 a2335cf9
Showing 16 changed files with 82 additions and 62 deletions Show diff stats
src/idl/dustem_cc.pro
src/idl/dustem_compute_stokes.pro
src/idl/dustem_filter2instru.pro
src/idl/dustem_fit_sed_qsed_used_readme.pro
src/idl/dustem_fit_sed_readme.pro
src/idl/dustem_init.pro
src/idl/dustem_instru2filters.pro
src/idl/dustem_read_align.pro
src/idl/dustem_read_all_web3p8.pro
src/idl/dustem_read_filters.pro
src/idl/dustem_read_qabs_lv.pro
src/idl/dustem_read_spin.pro
src/idl/dustem_write_align.pro
src/idl/dustem_write_all_web3p8.pro
src/idl/dustem_write_grain_web3p8.pro
src/idl/dustem_write_qabs.pro
@@ -132,7 +132,7 @@ FOR i=0L,Nbands-1 DO BEGIN
       'NUINU=CSTE': BEGIN
         spec_int=interpol(spec2,wavein2,*(st.use_wavelengths(ii)))
         num=integral(*(st.use_wavelengths(ii)),spec_int*(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii)))^2.,st.use_wmin(ii),st.use_wmax(ii),/double)
-        den=integral(*(st.use_wavelengths(ii)),(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii))),st.use_wmin(ii),st.use_wmax(ii),/double)
+        den=integral(*(st.use_wavelengths(ii)),(*(st.use_transmissions(ii)))/(*(st.use_wavelengths(ii))),st.use_wmin(ii),st.use_wmax(ii),/double)            
         cc(i)=num/den/spec0*(st.central_wavelengths(ii))
       END
       'MIPS': BEGIN
@@ -94,10 +94,10 @@ FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
 ENDFOR
+;stop
 ;-----------------------------------------
-
 ;INITIALIZING THE STOKES SEDS
 dustem_qsed = (*!dustem_data.sed).values * 0.
 dustem_used = dustem_qsed
@@ -105,7 +105,6 @@ dustem_used = dustem_qsed
 if not isarray(stI) THEN stop 
-
 ;Performing color corrections on the dustem spectra and generating the SEDs at the given filters 
 IF !dustem_do_cc NE 0 AND !dustem_never_do_cc EQ 0 THEN BEGIN
@@ -128,14 +127,14 @@ IF count_polsed NE 0 THEN BEGIN ; INCLUDE COLOR CORRECTIONS OF (Q,U) for dust &amp; 
     dustem_used(ind_polsed)=sused
   endif
   ;--------------------------------------------------------------------------
-  
 ENDIF
+
+
 ;For spectrum data points, interpolate in log-log
 ;Linear interpolation leads to wrong values, in particular where few
 ;wavelengths points exist in the model (long wavelengths).
-
 ind_spec=where((*!dustem_data.sed).filt_names EQ 'SPECTRUM',count_spec)
 IF count_spec NE 0 THEN BEGIN
   if isa(dustem_qsed) then dustem_qsed(ind_spec)=interpol(Q_spec,st.polsed.wav,(((*!dustem_data.sed).wav)(ind_spec)))  
@@ -143,7 +142,6 @@ IF count_spec NE 0 THEN BEGIN
 ENDIF
 out_st=st
-
 ;clean pointers
 heap_gc
@@ -52,7 +52,6 @@ ENDFOR
 instru=!dustem_instrument_description(ind1).instru
 the_end:
-
 RETURN,instru
 END
@@ -52,7 +52,7 @@ dustem_define_la_common
 ;=== initialise dustem
 ;dustem_init,mode='COMPIEGNE_ETAL2010',/pol
-use_mode='G17_MODELA' ; try with model D 
+use_mode='THEMIS';'G17_MODELA' ; try with model D 
 dustem_init,mode=use_mode,/pol
@@ -67,7 +67,7 @@ help,!run_pol
 ;filters=['DIRBE1','DIRBE2','DIRBE3','DIRBE4','IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3']
-filters=['IRAS1','IRAS2','IRAS3','IRAS4','PACS3','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3']
+filters=['IRAS1','IRAS2','IRAS3','IRAS4','PACS3','PILOT1','SPIRE1','SPIRE2','SPIRE3','HFI2','HFI3','HFI4','HFI5','HFI6','LFI1','LFI2','LFI3']
 Nfilt=n_elements(filters)
 sed=dustem_initialize_sed(Nfilt)
 sed.filter=filters
@@ -78,11 +78,11 @@ sed.instru=dustem_filter2instru(filters)
 ;INITIALIZING IQU SEDS and associated errors
 sed.StokesI=sed.wave*0
-;sed.StokesQ=sed.wave*0
-;sed.StokesU=sed.wave*0
+sed.StokesQ=sed.wave*0
+sed.StokesU=sed.wave*0
 sed.sigmaII=sed.wave*0
-;sed.sigmaQQ=sed.wave*0
-;sed.sigmaUU=sed.wave*0
+sed.sigmaQQ=sed.wave*0
+sed.sigmaUU=sed.wave*0
@@ -97,17 +97,23 @@ pd = [ $
              '(*!dustem_params).grains(0).mdust_o_mh',$  
              '(*!dustem_params).grains(1).mdust_o_mh',$     
              ;'(*!dustem_params).grains(1).mdust_o_mh',$    
-             '(*!dustem_params).grains(2).mdust_o_mh'$    
+             '(*!dustem_params).grains(2).mdust_o_mh',$    
              ;'dustem_plugin_synchrotron_2', $               ;Synchrotron amplitude
              ;'dustem_plugin_synchrotron_4', $               ;Synchrotron polarization angle
-             ;'dustem_plugin_modify_dust_polarization_2', $ ;Polarization angle applied to the Fortran dust model via a core plugin
-              ;'dustem_plugin_modify_spinning_polarization_1', $
-             ;'dustem_plugin_modify_spinning_polarization_2' $
+             'dustem_plugin_modify_dust_polarization_2',$ ;Polarization angle applied to the Fortran dust model via a core plugin
+             'dustem_plugin_modify_spinning_polarization_1', $
+             'dustem_plugin_modify_spinning_polarization_2' $
               ]                   
-;p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04,5e-3,20.,68.,0.01,30]    
+;p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04,5e-3,20.,68.,0.03,30.]    
+;p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04,68.,0.03,30.]
+;p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04];,0.01,30.]
+
+
+
+p_truth=[1.,0.17E-03,0.63E-03,0.51E-03,68.,0.03,30.]
+
-p_truth=[1.,7.8000E-04,7.8000E-04,1.6500E-04];,0.01,30.]
 ;p_truth=[68.]
@@ -131,8 +137,10 @@ fixed=[1,1,1,1,0];0,1,1]
 ;dustem_init_fixed_params,fpd,fiv
+iv=p_truth+[0.,5.000E-04,5.000E-04,5.000E-03,-10.,1.00e-2,10.]      ;shifted from solution ,10.
-iv=p_truth+[0.,5.000E-04,5.000E-04,5.000E-04];,5.e-4,1.]      ;shifted from solution ,10.
+
+;iv=p_truth+[0.,5.000E-04,5.000E-04,5.000E-04];,5.e-4,1.]      ;shifted from solution ,10.
 Npar=n_elements(pd)  
@@ -154,19 +162,19 @@ dustem_init_plugins,pd;, fpd=fpd
 sed.StokesI=dustem_compute_sed(p_truth,ste)   ;sst is not set on purpose, for dustem_compute_sed to compute the sed using fortran code
-sed.sigmaII=(sed.StokesI*0.3)^2
+sed.sigmaII=double((sed.StokesI*0.3)^2)
-; toto=dustem_compute_stokes(p_truth,ste,dustem_qsed,dustem_used) ;this procedure also allows for the extraction of the spectra
+toto=dustem_compute_stokes(p_truth,stee,dustem_qsed,dustem_used) ;this procedure also allows for the extraction of the spectra
-; sed.StokesQ=dustem_qsed
-; sed.sigmaQQ=(double((sed.StokesQ*0.03)^2)) > 1d-50
-; ;tes=where(finite(sed.sigmaQQ) eq 0,coun)
-; ;if coun ne 0 then sed.sigmaQQ(tes)=0.
+sed.StokesQ=dustem_qsed
+sed.sigmaQQ=(double((sed.StokesQ*0.03)^2)) > 1d-40
+;tes=where(finite(sed.sigmaQQ) eq 0,coun)
+;if coun ne 0 then sed.sigmaQQ(tes)=0.
-; sed.StokesU=dustem_used
-; sed.sigmaUU=(double((sed.StokesU*0.03)^2)) > 1d-50
-; ;tes=where(finite(sed.sigmaUU) eq 0,coun)
-; ;if coun ne 0 then sed.sigmaUU(tes)=0.
+sed.StokesU=dustem_used
+sed.sigmaUU=(double((sed.StokesU*0.03)^2)) > 1d-40
+;tes=where(finite(sed.sigmaUU) eq 0,coun)
+;if coun ne 0 then sed.sigmaUU(tes)=0.
 ;======================================================================================================================
@@ -71,7 +71,7 @@ ENDIF
 IF keyword_set(model) THEN BEGIN
   use_model=strupcase(model)
 ENDIF ELSE BEGIN
-  use_model='G17_MODELA';'COMPIEGNE_ETAL2010'    ;Default is last dustem model
+  use_model='THEMIS';G17_MODELA';'COMPIEGNE_ETAL2010'    ;Default is last dustem model
 ENDELSE
 IF keyword_set(png) THEN BEGIN
@@ -347,13 +347,11 @@ CASE use_model OF
       'THEMIS':BEGIN
          pd = [ $
              '(*!dustem_params).G0', $      ;G0
-             '(*!dustem_params).grains(0).mdust_o_mh',$  ;PAH0 mass fraction
-             '(*!dustem_params).grains(1).mdust_o_mh',$     ;PAH1 mass fraction
-             '(*!dustem_params).grains(2).mdust_o_mh', $    ;Gra
-             '(*!dustem_params).grains(3).mdust_o_mh', $    ;Gra
-             '(*!dustem_params).grains(4).mdust_o_mh', $    ;aSil
+             '(*!dustem_params).grains(0).mdust_o_mh',$  
+             '(*!dustem_params).grains(1).mdust_o_mh',$     
+             '(*!dustem_params).grains(2).mdust_o_mh', $  
              'dustem_plugin_continuum_2']                ;Intensity of NIR continuum
-        iv =   [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001]
+        iv =   [1.0,5.4e-4,2.33e-3,8.27e-3,0.001]
         Npar=n_elements(pd)
         ulimed=replicate(0,Npar)
         llimed=replicate(1,Npar)
@@ -169,8 +169,6 @@ defsysv, &#39;!dustem_params&#39;, ptr_new() ;Contains the values of all Desert Model pa
 ;ENDIF
-
-
 defsysv,'!dustem_filters',ptr_new()     ;filter information (filled by dustem_filter_init.pro)
 ;; IDL> help,*!dustem_filters,/str
 ;; ** Structure <1b95608>, 17 tags, length=8232, data length=8118, refs=1:
@@ -241,7 +239,7 @@ defsysv, &#39;!dustem_model&#39;, model
     END
     'THEMIS':BEGIN
       (*!dustem_inputs).grain='GRAIN_THEMIS.DAT'
-      (*!dustem_inputs).align='ALIGN_G17_MODELD.DAT'
+      (*!dustem_inputs).align='ALIGN_G17_MODELA.DAT'
     END
     ELSE:BEGIN
       message,'model '+model+' unknown',/continue
@@ -19,5 +19,6 @@ ENDIF ELSE BEGIN
 ENDELSE
 RETURN,filters
+stop
 END
@@ -41,16 +41,13 @@ readf,unit,str,format=&#39;(A100)&#39;
 xx=strsplit(str,' ',/regex,/extr)
 anisG0 = float(xx(0))
 ncurrent+=+1
-nalig=nlines-ncurrent
-
+nalig=nlines-ncurrent  ; why do we need this? 
 ;stop
 one_st={aligned:0.,law:'',athresh:0.,plev:0.,pstiff:0.,TdsTg:0.,a0:0.,rvcut:0.}
 ;Ngrains=st_grains.Ngrains
 Ngrains=n_elements(st_grains)
-
 st=replicate(one_st,Ngrains)
-
 ;st.aligned = stregex(st_grains.grains.type_keywords, 'pol', /bool) 
 ;count=0
@@ -92,8 +89,9 @@ FOR i=0L,Ngrains-1 DO BEGIN
     IF !run_univ eq 1 THEN BEGIN
       st = replicate(st(i),Ngrains)
       ;st.aligned = stregex(st_grains.grains.type_keywords, 'pol', /bool) 
-      st.aligned = stregex(st_grains(i).type_keywords, 'pol', /bool) 
+      st(i).aligned = stregex(st_grains(i).type_keywords, 'pol', /bool)
       break
+      ;stop  
     ENDIF
   ENDIF
 ENDFOR
@@ -150,7 +148,6 @@ close,unit
 free_lun,unit
 full_st={keywords:key_str,anisG0:anisG0,gamma:0,grains:st}
-
 RETURN,full_st
 END
@@ -42,9 +42,9 @@ IF keyword_set(help) THEN BEGIN
   goto,the_end
 ENDIF
-dir_in_dat=dir_in+'/data/'
-dir_in_qabs=dir_in+'/oprop/'
-dir_in_capa=dir_in+'/hcap/'
+dir_in_dat=dir_in+'data/'
+dir_in_qabs=dir_in+'oprop/'
+dir_in_capa=dir_in+'hcap/'
 file=dir_in_dat+(*!dustem_inputs).grain
 ;st_grains=dustem_read_grain(file,silent=silent)
@@ -64,7 +64,6 @@ st_gas=dustem_read_gas(file,silent=silent)
 st_qabs=ptrarr(Ngrains)
 FOR i=0L,Ngrains-1 DO BEGIN
-  ;stop
   read_densities=0
   IF st_grains[i].rho LE 0 THEN read_densities=1
   Qabs_file=dir_in_qabs+'Q_'+st_grains(i).grain_type+'.DAT'
@@ -78,6 +78,7 @@ dir_lfi=filter_dir+&#39;LFI&#39;+&#39;/&#39;
 dir_wmap=filter_dir+'WMAP'+'/'
 dir_spire=filter_dir+'SPIRE'+'/'
 dir_pacs=filter_dir+'PACS'+'/'
+dir_pilot=filter_dir+'PILOT'+'/'
 dir_akari=filter_dir+'AKARI'+'/'
 dir_spass=filter_dir+'SPASS'+'/'
 dir_nika2=filter_dir+'NIKA2'+'/'
@@ -125,6 +126,29 @@ pacs_struct={Name:&#39;PACS&#39;,Nbands:Nband,filter_names:filt_names,central_wavelength
             use_frequencies:use_nu,use_transmissions:use_T,use_wavelengths:use_w,use_wmin:dblarr(Nband),use_wmax:dblarr(Nband), $
             cc_method:'dustem_cc_pacs'}
+
+
+
+;======================================================
+;==== PILOT
+;======================================================
+IF print_messages THEN message,'Reading PILOT filter',/continue
+
+;Read WMAP filter Transmission
+filt_names=['PILOT1']
+wfilt_pilot=dustem_filter2wav(filt_names)
+Nband=n_elements(filt_names)
+nu_ptr=ptrarr(Nband) & T_ptr=ptrarr(Nband) & wav_ptr=ptrarr(Nband)
+use_nu=ptrarr(Nband) & use_T=ptrarr(Nband) & use_w=ptrarr(Nband)
+stP=read_xcat(dir_pilot+'pilot_transmission.xcat',/silent)
+Trans=stP.transmission & lamb=stP.wavelength;cmic/(stP.frequency*1e9)  ;microns
+T_ptr(0)=ptr_new(Trans/max(Trans)) & nu_ptr(0)=ptr_new(cmic/lamb) & wav_ptr(0)=ptr_new(lamb)
+
+pilot_struct={Name:'PILOT',Nbands:Nband,filter_names:filt_names,central_wavelengths:wfilt_pilot,central_frequencies:cmic/wfilt_pilot, $
+            filter_wavelengths:wav_ptr,filter_frequencies:nu_ptr,filter_transmissions:T_ptr, $
+            use_frequencies:use_nu,use_transmissions:use_T,use_wavelengths:use_w,use_wmin:dblarr(Nband),use_wmax:dblarr(Nband), $
+            cc_method:'dustem_cc_pilot'}
+            
 ;======================================================
 ;==== Herschel SPIRE
 ;======================================================
@@ -937,6 +961,7 @@ dm_filter_struct={IRAC:irac_struct,MIPS:mips_struct,MSX:msx_struct,IRAS:iras_str
                   HFI:hfi_struct,LFI:lfi_struct, $
                   WMAP:wmap_struct, $
                   SPIRE:spire_struct,PACS:pacs_struct, $
+                  PILOT:pilot_struct, $
                   AKARI:akari_struct, $
                   BOLOCAM:bolocam_struct, $
                   WISE:wise_struct, $
@@ -24,7 +24,7 @@ Ncol=Nsizes
 IF keyword_set(read_densities) THEN BEGIN
   readf,unit,st
-  stop
+  ;stop
   IF st EQ '#### Qabs ####' THEN BEGIN
     message,'Problem with size dependent densities in Qabs.DAT',/continue
     ;stop
@@ -77,6 +77,7 @@ close,unit
 free_lun,unit
 qsca_str=sts(0:nlines-1)
+;stop
 Nlines=n_elements(qsca_str)
@@ -30,11 +30,11 @@ WHILE NOT EOF(unit) DO BEGIN
   tmp = [tmp, str]
 ENDWHILE
 ;==save and return structure
-full_st={file:file,fspin:tmp,comments:comments}
-  
+full_st={file:file,fspin:tmp,comments:comments}  
+
 close,unit
 free_lun,unit
-
+;stop
 RETURN,full_st
 END
@@ -17,9 +17,7 @@ c(10)=&#39;# alignmentlaw (idg, rat, par), parameters (par1 par2 par3)&#39;
 c(11)='#'
 openw,unit,file,/get_lun
-
 FOR ii=0,Ncomments-1 DO printf,unit,c(ii)
-
 printf,unit,st_align.keywords,format='(A-40)'
 printf,unit,st_align.anisG0,format='(E10.2)'
@@ -34,6 +32,5 @@ endelse
 close,unit
 free_lun,unit
-
 END
 PRO dustem_write_all_web3p8,st,dir_out
-dir_out_dat=dir_out+'data/'
-dir_out_qabs=dir_out+'oprop/'
-dir_out_capa=dir_out+'hcap/'
+dir_out_dat=dir_out+'/data/'
+dir_out_qabs=dir_out+'/oprop/'
+dir_out_capa=dir_out+'/hcap/'
 ;== ISRF
 file_out=dir_out_dat+'ISRF.DAT'
@@ -70,7 +70,6 @@ ENDFOR
 close,unit
 free_lun,unit
-;stop
 END
@@ -8,7 +8,7 @@ Qstring=&#39;QABS_&#39;
 Estring='.DAT'
 Nfiles=n_elements(st)
 frmt="(50E15.5)"
-
+stop
 FOR i=0L,Nfiles-1 DO BEGIN
   filename=dir+Qstring+(*st(i)).material+Estring
   Nsizes=n_elements((*st(i)).sizes)