forced ordering of filter wavelengths, to avoid gdl message when interpolating

Jean-Philippe Bernard
1 parent 3dd007df
Showing 2 changed files with 76 additions and 18 deletions Show diff stats
src/idl/dustem_filter_init.pro
src/idl/dustem_read_filters.pro
@@ -46,35 +46,52 @@ IF keyword_set(help) THEN BEGIN
   goto,the_end
 ENDIF
+order_by_wavelength=1   ;to order all transmissions by increasing wavelength
+
 ;=== Read all filters
 st=dustem_read_filters()
 !dustem_filters=ptr_new(st)
+;======================================================
+;==== order by increasing wavelength if needed
+;======================================================
+IF order_by_wavelength EQ 1 THEN BEGIN
+  filter_names=tag_names((*!dustem_filters))
+  FOR i=0L,n_tags(*!dustem_filters)-1 DO BEGIN
+    FOR j=0L,(*!dustem_filters).(i).Nbands-1 DO BEGIN
+      ;message,'Ordering '+filter_names[i]+strtrim(j+1,2),/continue
+      transmissions=*((*!dustem_filters).(i).filter_transmissions(j))
+      wavelengths=*((*!dustem_filters).(i).filter_wavelengths(j))
+      frequencies=*((*!dustem_filters).(i).filter_frequencies(j))
+      order=sort(wavelengths)
+      wavelengths=wavelengths[order]
+      transmissions=transmissions[order]
+      frequencies=frequencies[order]
+      un=uniq(wavelengths)
+      wavelengths=wavelengths[un]
+      transmissions=transmissions[un]
+      frequencies=frequencies[un]
+      ((*!dustem_filters).(i).filter_transmissions(j))=ptr_new(transmissions)
+      ((*!dustem_filters).(i).filter_wavelengths(j))=ptr_new(wavelengths)
+      ((*!dustem_filters).(i).filter_frequencies(j))=ptr_new(frequencies)
+    ENDFOR
+  ENDFOR
+ENDIF
+
 ;== run standard model to set the wavelengths used by Dustem
 ;== This is needed to integrate the spectrum accurately in the filters
 IF not keyword_set(dir) THEN BEGIN
   CASE !dustem_which OF
-;  CASE getenv('DUSTEM_WHICH') OF
     'DESERT':dir_in=!dustem_soft_dir+'DESERT_POST_ISO_MORE/'
     'COMPIEGNE': dir_in=!dustem_soft_dir+'MC_DAT/'
     'VERSTRAETE': dir_in=!dustem_soft_dir+'d_3.5/'
     'WEB3p8': dir_in=!dustem_soft_dir
-;    'WEB3p8': dir_in=!dustem_soft_dir+'/src/dustem3.8_web/'
-;    'DESERT':dir_in=getenv('DUSTEM_SOFT_DIR')+'DESERT_POST_ISO_MORE/'
-;    'COMPIEGNE': dir_in=getenv('DUSTEM_SOFT_DIR')+'MC_DAT/'
-;    'VERSTRAETE': dir_in=getenv('DUSTEM_SOFT_DIR')+'d_3.5/'
-;    'WEB3p8': dir_in=getenv('DUSTEM_SOFT_DIR')+'/src/dustem3.8_web/'
     ELSE: dir_in=!dustem_soft_dir+'/Data/les_DAT/'
-;    ELSE: dir_in=getenv('DUSTEM_SOFT_DIR')+'/Data/les_DAT/'
   ENDCASE
 ENDIF ELSE BEGIN
   dir_in=!dustem_soft_dir
-;  dir_in=!dustem_soft_dir+'/src/dustem3.8_web/'
-;  dir_in=getenv('DUSTEM_SOFT_DIR')+'/src/dustem3.8_web/'
-;  dir_in=getenv('DUSTEM_SOFT_DIR')+'/Data/'+dir
 ENDELSE
-;CASE getenv('DUSTEM_WHICH') OF
 CASE !dustem_which OF
   'WEB3p8':BEGIN
      dir_in_qabs=dir_in+'/oprop/'
@@ -84,7 +101,6 @@ CASE !dustem_which OF
      dustem_waves=st_lambda.lambda
   END
   ELSE: BEGIN
-;    dir_out=getenv('DUSTEM_DAT')
     dir_out=!dustem_dat
     st_model=dustem_read_all(dir_in,/silent)
     dustem_write_all,st_model,dir_out
@@ -95,10 +111,11 @@ ENDCASE
 cmic=2.9979246D14
-;== Fill in the !dustem_filters structure
+;== interpolating transmission on use_waavelengths
+filter_names=tag_names((*!dustem_filters))
 FOR i=0L,n_tags(*!dustem_filters)-1 DO BEGIN
   FOR j=0L,(*!dustem_filters).(i).Nbands-1 DO BEGIN
-;    all_waves=[*((*!dustem_filters).(i).filter_wavelengths(j)),st.sed.wav]
+    ;message,'Doing '+filter_names[i]+strtrim(j+1,2),/continue
     all_waves=[*((*!dustem_filters).(i).filter_wavelengths(j)),dustem_waves]
     order=sort(all_waves)
     all_waves=all_waves(order)
@@ -109,9 +126,11 @@ FOR i=0L,n_tags(*!dustem_filters)-1 DO BEGIN
     all_waves=all_waves(ind)
     ind=where(all_waves GE min(*((*!dustem_filters).(i).filter_wavelengths(j))) AND $
               all_waves LE max(*((*!dustem_filters).(i).filter_wavelengths(j))),count)
-    all_waves=all_waves(ind)
+    all_waves=all_waves[ind]
+    all_frequencies=cmic/all_waves
+    ;==== fill up structure and interpol in wavelengths
     (*!dustem_filters).(i).use_wavelengths(j)=ptr_new(all_waves)
-    (*!dustem_filters).(i).use_frequencies(j)=ptr_new(cmic/all_waves)
+    (*!dustem_filters).(i).use_frequencies(j)=ptr_new(all_frequencies)
     (*!dustem_filters).(i).use_transmissions(j)=ptr_new(interpol(*((*!dustem_filters).(i).filter_transmissions(j)), $
                                                                  *((*!dustem_filters).(i).filter_wavelengths(j)), $
                                                                  *((*!dustem_filters).(i).use_wavelengths(j))))
@@ -50,6 +50,7 @@ ccm=2.9979246D10
 ;which_hfi_filters='GL'
 which_hfi_filters='OFFICIAL'
+print_messages=0      ;to print filter names when reading filters
 ;======================================================
 ;==== directories containing filter information
@@ -85,7 +86,7 @@ dir_irs=filter_dir+&#39;IRS&#39;+&#39;/&#39;
 ;======================================================
 ;==== Herschel PACS
 ;======================================================
-
+IF print_messages THEN message,'Reading PACS filters',/continue
 ;Read Herschel PACS filter Transmission
 ;wfilt_pacs=[75,110,170.] ;microns
 filt_names=['PACS1','PACS2','PACS3']
@@ -124,6 +125,7 @@ pacs_struct={Name:&#39;PACS&#39;,Nbands:Nband,filter_names:filt_names,central_wavelength
 ;==== Herschel SPIRE
 ;======================================================
 ;stop
+IF print_messages THEN message,'Reading SPIRE filters',/continue
 filt_names=['SPIRE1','SPIRE2','SPIRE3']
 wfilt_spire=dustem_filter2wav(filt_names)
 Nband=n_elements(filt_names)
@@ -162,6 +164,7 @@ spire_struct={Name:&#39;SPIRE&#39;,Nbands:Nband,filter_names:filt_names,central_waveleng
 ;======================================================
 ;==== WMAP
 ;======================================================
+IF print_messages THEN message,'Reading WMAP filters',/continue
 ;Read WMAP filter Transmission
 filt_names=['WMAP1','WMAP2','WMAP3','WMAP4','WMAP5']
@@ -194,6 +197,7 @@ wmap_struct={Name:&#39;WMAP&#39;,Nbands:Nband,filter_names:filt_names,central_wavelength
 ;======================================================
 ;==== Spitzer IRAC
 ;======================================================
+IF print_messages THEN message,'Reading IRAC filters',/continue
 ;Read Spitzer IRAC filter Transmission
 ;Strongly inspired from N. Flagey's routine get_irac_flux.pro
@@ -217,6 +221,7 @@ irac_struct={Name:&#39;Spitzer-IRAC&#39;,Nbands:Nband,filter_names:filt_names,central_wa
 ;======================================================
 ;==== Spitzer MIPS
 ;======================================================
+IF print_messages THEN message,'Reading MIPS filters',/continue
 ;Read Spitzer MIPS filter Transmission
 ;Strongly inspired from N. Flagey's routine get_mips_flux.pro
@@ -239,6 +244,7 @@ mips_struct={Name:&#39;Spitzer-MIPS&#39;,Nbands:Nband,filter_names:filt_names,central_wa
 ;======================================================
 ;==== MSX
 ;======================================================
+IF print_messages THEN message,'Reading MSX filters',/continue
 ;Read MSX filter Transmission
 ;Strongly inspired from N. Flagey's routine get_msx_flux.pro
@@ -262,6 +268,8 @@ msx_struct={Name:&#39;MSX&#39;,Nbands:Nband,filter_names:filt_names,central_wavelengths:
 ;======================================================
 ;==== BOLOCAM
 ;======================================================
+IF print_messages THEN message,'Reading BOLOCAM filters',/continue
+
 filt_names=['BOLOCAM1']
 wfilt=dustem_filter2wav(filt_names)
 Nband=n_elements(filt_names)
@@ -277,6 +285,7 @@ bolocam_struct={Name:&#39;BOLOCAM&#39;,Nbands:Nband,filter_names:filt_names,central_wave
 ;======================================================
 ;==== Pronaos-SPM
 ;======================================================
+IF print_messages THEN message,'Reading SPM filters',/continue
 ;Read SPM filter Transmission
 filt_names=['SPM1','SPM2','SPM3','SPM4']
@@ -307,6 +316,7 @@ spm_struct={Name:&#39;PRONAOS-SPM&#39;,Nbands:Nband,filter_names:filt_names,central_wave
 ;======================================================
 ;==== Akari
 ;======================================================
+IF print_messages THEN message,'Reading Akari filters',/continue
 ;written by Ryou Ohsawa
 ;Bugs corrected by Deborah Paradis
@@ -338,6 +348,8 @@ akari_struct={Name:&#39;AKARI&#39;,Nbands:Nband,filter_names:filt_names,central_waveleng
 ;======================================================
 ;==== IRAS
 ;======================================================
+IF print_messages THEN message,'Reading IRAS filters',/continue
+
 Nband=4
 filt_names=['IRAS1','IRAS2','IRAS3','IRAS4']
 wfilt_iras=dustem_filter2wav(filt_names)
@@ -361,6 +373,8 @@ iras_struct={Name:&#39;IRAS&#39;,Nbands:Nband,filter_names:filt_names,central_wavelength
 ;======================================================
 ;==== DIRBE
 ;======================================================
+IF print_messages THEN message,'Reading DIRBE filters',/continue
+
 file=dir_dirbe+'dirbe_response.asc'
 readcol,file,wavemic,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,/silent
 ind=where(wavemic GT 1 AND wavemic LT 500,count)
@@ -391,6 +405,8 @@ dirbe_struct={Name:&#39;DIRBE&#39;,Nbands:Nband,filter_names:filt_names,central_waveleng
 ;======================================================
 ;==== ARCHEOPS
 ;======================================================
+IF print_messages THEN message,'Reading ARCHEOPS filters',/continue
+
 Nband=4
 filt_names=['ARCHEOPS4','ARCHEOPS3','ARCHEOPS2','ARCHEOPS1']
 wfilt_archeops=dustem_filter2wav(filt_names)
@@ -421,6 +437,7 @@ archeops_struct={Name:&#39;ARCHEOPS&#39;,Nbands:Nband,filter_names:filt_names,central_wa
 ;Filters for each detectors are just averaged simply
 ;PSBs are taken the same way as other bolometers, which
 ;is probably wrong for Intensity (should devide them by 2) ...
+IF print_messages THEN message,'Reading HFI filters',/continue
 Nband=6
 filt_names=['HFI6','HFI5','HFI4','HFI3','HFI2','HFI1']
@@ -504,6 +521,7 @@ hfi_struct={Name:&#39;HFI&#39;,Nbands:Nband,filter_names:filt_names,central_wavelengths:
 ;======================================================
 ;Obtained from Marco Bersaneli, via Cyril Roset
 ;Filters for each detectors are just averaged simply
+IF print_messages THEN message,'Reading LFI filters',/continue
 Nband=3
 filt_names=['LFI3','LFI2','LFI1']
@@ -540,6 +558,8 @@ lfi_struct={Name:&#39;LFI&#39;,Nbands:Nband,filter_names:filt_names,central_wavelengths:
 ;======================================================
 ;==== ISOCAM
 ;======================================================
+IF print_messages THEN message,'Reading ISOCAM filters',/continue
+
 file=dir_iso+'cam_filter_transmission.save'
 restore,file
 Nband_sw=11
@@ -582,6 +602,8 @@ isocam_struct={Name:&#39;ISOCAM&#39;,Nbands:Nband,filter_names:filt_names, $
 ;======================================================
 ;==== ISOphotP
 ;======================================================
+IF print_messages THEN message,'Reading ISOPHOTP filters',/continue
+
 ;Obtained from R. Laureijs email to JPB on March 10 2011
 filt_names='ISOPHOTP'+['1',   '2',  '3',   '4',  '5',  '6', '7',   '8',   '9',   '10','11','12','13','14']
 band_names='P*'+      ['3_29','3_6','4_85','7_3','7_7','10','11_3','11_5','12_8','16','20','25','60','100']+'FM_R.DAT'
@@ -606,6 +628,8 @@ isophotP_struct={Name:&#39;ISOPHOTP&#39;,Nbands:Nband,filter_names:filt_names,central_wa
 ;======================================================
 ;==== ISOphotC
 ;======================================================
+IF print_messages THEN message,'Reading ISOPHOTC filters',/continue
+
 ;Obtained from R. Laureijs email to JPB on March 10 2011
 filt_names='ISOPHOTC'+['1', '2', '3', '4', '5',  '6',  '7',  '8',  '9',  '10', '11']
 band_names='C_*'+     ['50','60','70','90','100','105','120','135','160','180','200']+'FM_R.DAT'
@@ -630,6 +654,8 @@ isophotC_struct={Name:&#39;ISOPHOTC&#39;,Nbands:Nband,filter_names:filt_names,central_wa
 ;======================================================
 ;==== FIRAS
 ;======================================================
+IF print_messages THEN message,'Reading equivalent FIRAS filters',/continue
+
 ;Equivalent FIRAS filters (used in Boulanger et al 90)
 ;Defined as top-hat filters with 30% bandwidth (although not given in
 ;Boulanger et al).
@@ -661,6 +687,8 @@ efiras_struct={Name:&#39;EFIRAS&#39;,Nbands:Nband,filter_names:filt_names,central_wavele
 ;======================================================
 ;==== WISE
 ;======================================================
+IF print_messages THEN message,'Reading WISE filters',/continue
+
 Nband=4
 filt_names=['WISE1','WISE2','WISE3','WISE4']
 wfilt_wise=dustem_filter2wav(filt_names)
@@ -684,6 +712,8 @@ wise_struct={Name:&#39;WISE&#39;,Nbands:Nband,filter_names:filt_names,central_wavelength
 ;======================================================
 ;==== LABOCA
 ;======================================================
+IF print_messages THEN message,'Reading LABOCA filters',/continue
+
 Nband=1
 filt_names=['LABOCA1']
 wfilt_laboca=dustem_filter2wav(filt_names)
@@ -702,6 +732,8 @@ laboca_struct={Name:&#39;LABOCA&#39;,Nbands:Nband,filter_names:filt_names,central_wavele
 ;======================================================
 ;==== GISMO
 ;======================================================
+IF print_messages THEN message,'Reading GISMO filters',/continue
+
 Nband=1
 filt_names=['GISMO1']
 wfilt_gismo=dustem_filter2wav(filt_names)
@@ -720,6 +752,8 @@ gismo_struct={Name:&#39;GISMO&#39;,Nbands:Nband,filter_names:filt_names,central_waveleng
 ;======================================================
 ;==== SPASS
 ;======================================================
+IF print_messages THEN message,'Reading SPASS filters',/continue
+
 Nband=1
 filt_names=['SPASS1']
 wfilt_spass=dustem_filter2wav(filt_names)
@@ -741,6 +775,8 @@ spass_struct={Name:&#39;SPASS&#39;,Nbands:Nband,filter_names:filt_names,central_waveleng
 ;======================================================
 ;==== NIKA2
 ;======================================================
+IF print_messages THEN message,'Reading NIKA2 filters',/continue
+
 ;# NIKA1200
 ;# photon
 ;# Average of array 1 and array3 of NIKA2, 1.2mm arrays
@@ -772,6 +808,8 @@ nika2_struct={Name:&#39;NIKA2&#39;,Nbands:Nband,filter_names:filt_names,central_waveleng
 ;======================================================
 ;==== SCUBA2
 ;======================================================
+IF print_messages THEN message,'Reading SCUBA2 filters',/continue
+
 Nband=1
 filt_names=['SCUBA21']
 wfilt_spass=dustem_filter2wav(filt_names)
@@ -792,6 +830,8 @@ scuba2_struct={Name:&#39;SCUBA2&#39;,Nbands:Nband,filter_names:filt_names,central_wavele
 ;======================================================
 ;==== IRS
 ;======================================================
+IF print_messages THEN message,'Reading IRS filters',/continue
+
 ;# IRS16
 ;# energy
 ;# Spitzer IRS 16 micron
@@ -837,7 +877,6 @@ dm_filter_struct={IRAC:irac_struct,MIPS:mips_struct,MSX:msx_struct,IRAS:iras_str
 ;======================================================
 ;==== plot transmissions
 ;======================================================
-
 IF keyword_set(plot_it) THEN BEGIN
   xtit=textoidl('\lambda (\mum)')
   ytit='Filter transmission'