From f0ac8429e6c266779ab6e754fcf9ed03d33b0a04 Mon Sep 17 00:00:00 2001 From: Jean-Philippe Bernard Date: Wed, 14 Jul 2021 12:00:00 +0200 Subject: [PATCH] forced ordering of filter wavelengths, to avoid gdl message when interpolating --- src/idl/dustem_filter_init.pro | 51 +++++++++++++++++++++++++++++++++++---------------- src/idl/dustem_read_filters.pro | 43 +++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 76 insertions(+), 18 deletions(-) diff --git a/src/idl/dustem_filter_init.pro b/src/idl/dustem_filter_init.pro index f60d452..3fc40a3 100644 --- a/src/idl/dustem_filter_init.pro +++ b/src/idl/dustem_filter_init.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)))) diff --git a/src/idl/dustem_read_filters.pro b/src/idl/dustem_read_filters.pro index a4110a6..5809e35 100644 --- a/src/idl/dustem_read_filters.pro +++ b/src/idl/dustem_read_filters.pro @@ -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+'IRS'+'/' ;====================================================== ;==== 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:'PACS',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:'SPIRE',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:'WMAP',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:'Spitzer-IRAC',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:'Spitzer-MIPS',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:'MSX',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:'BOLOCAM',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:'PRONAOS-SPM',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:'AKARI',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:'IRAS',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:'DIRBE',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:'ARCHEOPS',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:'HFI',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:'LFI',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:'ISOCAM',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:'ISOPHOTP',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:'ISOPHOTC',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:'EFIRAS',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:'WISE',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:'LABOCA',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:'GISMO',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:'SPASS',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:'NIKA2',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:'SCUBA2',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' -- libgit2 0.21.2