OV-GSO-DC / dustem-wrapper_idl

Commit 0c48096dd5ec56117117b7065f5f07bd625fb572

Authored by Jean-Philippe Bernard
2 parents 2f404849 e475799b
Exists in master

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

Showing 2 changed files with 315 additions and 12 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_cc_print.pro
  Diff comments   View file @0c48096
1   -PRO dustem_cc_print,filters,restart=restart,help=help
  1 +PRO dustem_cc_print,filters,norestart=norestart,help=help
2 2  
3 3 ;+
4 4 ; NAME:
... ... @@ -34,8 +34,7 @@ PRO dustem_cc_print,filters,restart=restart,help=help
34 34 ;
35 35 ; ACCEPTED KEY-WORDS:
36 36 ; help = If set, print this help
37   -; restart = If set, reinitialise DustEMWrap and use a new
38   -; wavelength vector for integration.
  37 +; norestart = If set, do not reinitialise DustEMWrap
39 38 ;
40 39 ; COMMON BLOCKS:
41 40 ; None
... ... @@ -53,8 +52,8 @@ PRO dustem_cc_print,filters,restart=restart,help=help
53 52 ; PROCEDURES AND SUBROUTINES USED:
54 53 ;
55 54 ; EXAMPLES
56   -; dustem_cc_print,['MIRI1','MIRI6'],/restart
57   -; dustem_cc_print,['IRAS1','IRAC2']
  55 +; dustem_cc_print,['MIRI1','MIRI6']
  56 +; dustem_cc_print,['IRAS1','IRAC2'],/norestart
58 57 ;
59 58 ; MODIFICATION HISTORY:
60 59 ; Written by JPB Apr-2011
... ... @@ -66,8 +65,14 @@ IF keyword_set(help) THEN BEGIN
66 65 goto,the_end
67 66 ENDIF
68 67  
  68 +;== Check if user gave us filters
  69 + if not keyword_set(filters) then begin
  70 + message,'You must provide a list of filters',/info
  71 + stop
  72 + end
  73 +
69 74 ;== Initialize DustemWrap if restart requested
70   - if keyword_set(restart) then begin
  75 + if not keyword_set(norestart) then begin
71 76 dustem_init
72 77 Nwav=50000 & wavmin=0.1 & wavmax=500.
73 78 wavs=findgen(Nwav)/(Nwav-1)*(wavmax-wavmin)+wavmin
... ... @@ -78,12 +83,6 @@ ENDIF
78 83 wavs=st.sed.wav
79 84 end
80 85  
81   -;== Check if user gave us filters
82   - if not keyword_set(filters) then begin
83   - message,'You must provide a list of filters',/info
84   - stop
85   - end
86   -
87 86 Nfilt=n_elements(filters)
88 87 wavref=dustem_filter2wav(filters)
89 88  
... ...
src/idl/dustem_cc_verify.pro 0 → 100644
  Diff comments   View file @0c48096
... ... @@ -0,0 +1,304 @@
  1 +PRO dustem_cc_verify,instrument $
  2 + ,nwav=nwav $
  3 + ,wavrng=wavrng $
  4 + ,logspace=logspace $
  5 + ,avg_thresh=avg_thresh $
  6 + ,max_thresh=max_thresh $
  7 + ,bb=bb,mbb=mbb,power=power $
  8 + ,help=help
  9 +
  10 +;+
  11 +; NAME:
  12 +; dustem_cc_verify
  13 +;
  14 +; PURPOSE:
  15 +; Compares color corrections for a given spectrum in a given filter
  16 +; to fiducial values (in DATA/FIDUCIAL_CC)
  17 +; This is a sanity-checking routine.
  18 +; Output should be compared to the color-correction tables published by different
  19 +; instrument teams.
  20 +; e.g. 1 MIPS HandBook, page 98
  21 +; http://irsa.ipac.caltech.edu/data/SPITZER/docs/mips/mipsinstrumenthandbook/MIPS_Instrument_Handbook.pdf
  22 +; http://irsa.ipac.caltech.edu/data/SPITZER/docs/mips/mipsinstrumenthandbook/51/
  23 +; AND
  24 +; e.g. 2 IRAS Explanatory Supplement, Table VI.C.6
  25 +; https://irsa.ipac.caltech.edu/IRASdocs/exp.sup/ch6/tabsupC6.html
  26 +;
  27 +; CATEGORY:
  28 +; DustEMWrap, Distributed, Low-Level, User convenience
  29 +;
  30 +; CALLING SEQUENCE:
  31 +; dustem_cc_verify,instrument[,help=help]
  32 +;
  33 +; INPUTS:
  34 +; instrument: string, instrument which color corrections will be compared
  35 +; (currently only accepts 1 instrument)
  36 +;
  37 +; OPTIONAL INPUT PARAMETERS:
  38 +; None
  39 +;
  40 +; OUTPUTS:
  41 +;
  42 +; OPTIONAL OUTPUT PARAMETERS:
  43 +;
  44 +; ACCEPTED KEY-WORDS:
  45 +; help = If set, print this help
  46 +;
  47 +; COMMON BLOCKS:
  48 +; None
  49 +;
  50 +; SIDE EFFECTS:
  51 +; None
  52 +;
  53 +; RESTRICTIONS:
  54 +; The DustEM fortran code must be installed
  55 +; The DustEMWrap IDL code must be installed
  56 +; If restart is not set, this routine should be called after the
  57 +; successful completion of a DustEMWrap run (i.e. without quitting the
  58 +; IDL session) since it will use existing system variables
  59 +;
  60 +; PROCEDURES AND SUBROUTINES USED:
  61 +;
  62 +; EXAMPLES
  63 +; dustem_cc_verify,'WISE'
  64 +;
  65 +; MODIFICATION HISTORY:
  66 +; Written by AH Jan 2023
  67 +; Evolution details on the DustEMWrap gitlab.
  68 +; See http://dustemwrap.irap.omp.eu/ for FAQ and help.
  69 +
  70 +IF keyword_set(help) THEN BEGIN
  71 + doc_library,'dustem_cc_verify'
  72 + goto,the_end
  73 +ENDIF
  74 +
  75 +use_avg_thresh=0.01 ; a mean difference of >=1% between our CCs and the tabulated values will be recorded
  76 +use_max_thresh=0.02 ; a max difference of >=2% between our CCs and the tabulated values will be recorded
  77 +if keyword_set(avg_thresh) then use_avg_thresh=avg_thresh
  78 +if keyword_set(max_thresh) then use_max_thresh=max_thresh
  79 +use_lo=25
  80 +use_hi=75
  81 +
  82 +;== Initialize DustemWrap
  83 +dustem_init
  84 +use_Nwav=500000 & use_wavmin=0.01 & use_wavmax=5000.
  85 +if keyword_set(Nwav) then use_Nwav=Nwav
  86 +
  87 +if keyword_set(wavrng) then begin
  88 + use_wavmin=wavrng[0]
  89 + use_wavmax=wavrng[1]
  90 +endif
  91 +wavs=findgen(use_Nwav)/(use_Nwav-1)*(use_wavmax-use_wavmin)+use_wavmin
  92 +if keyword_set(logspace) then wavs=10^(alog10(use_wavmin)+alog10(use_wavmax)*(dindgen(use_Nwav+1)/use_Nwav))
  93 +;restart=0
  94 +;end else begin
  95 +; st=dustem_run()
  96 +; wavs=st.sed.wav
  97 +;end
  98 +
  99 +;== Check user gave us an instrument
  100 +if n_elements(instrument) ne 1 then begin
  101 + message,'You must provide one instrument name',/info
  102 + stop
  103 +end
  104 +
  105 +;== Get associated filters and central wavelengths
  106 +filters=dustem_instru2filters(instrument)
  107 +Nfilt=n_elements(filters)
  108 +wavref=dustem_filter2wav(filters)
  109 +
  110 +;== potential fiducial CC filenames
  111 +fidcc_pl_files=file_search(!dustem_wrap_soft_dir+'/Data/FIDUCIAL_CC/'+strlowcase(instrument)+'*_powerlaw.xcat',count=nplfiles)
  112 +fidcc_bb_files=file_search(!dustem_wrap_soft_dir+'/Data/FIDUCIAL_CC/'+strlowcase(instrument)+'*_bb.xcat',count=nbbfiles)
  113 +; problem in MBB formula?
  114 +;fidcc_mbb_files=file_search(!dustem_wrap_soft_dir+'/Data/FIDUCIAL_CC/'+strlowcase(instrument)+'*_mbb.xcat',count=nmbbfiles)
  115 +
  116 +problem_cc=['NONE']
  117 +
  118 +tt=where(tag_names((*!dustem_filters)) eq instrument)
  119 +;== Cross-check the power law CCs
  120 +
  121 +pl_ccs:
  122 +for i=0,nplfiles-1 do begin
  123 + flag_bad_cc=0
  124 + fidcc_pl_str=read_xcat(fidcc_pl_files[i],/silent)
  125 + pows=fidcc_pl_str.alpha
  126 + Nv=n_elements(pows)
  127 + Kmatrix=fltarr(Nfilt,Nv)
  128 + Kmatrix_fid=fltarr(Nfilt,Nv)
  129 + FOR j=0,Nfilt-1 DO Kmatrix_fid[j,*]=fidcc_pl_str.(j+1) ; alpha is the first tag (j=0)
  130 + FOR k=0L,Nv-1 DO BEGIN
  131 + FOR ff=0L,Nfilt-1 DO BEGIN
  132 + fwavs=*((*!dustem_filters).(tt).filter_wavelengths[ff])
  133 + spec=fwavs^(-1.*pows[k])
  134 + sed=dustem_cc(fwavs,spec,filters[ff],cc=cc,fluxconv=fluxconv,help=help)
  135 + Kmatrix[ff,k]=cc
  136 + ENDFOR
  137 + ENDFOR
  138 +
  139 + goodidx=where(Kmatrix_fid ne !indef and Kmatrix ne !indef and Kmatrix_fid ne !values.f_nan and Kmatrix ne !values.f_nan,ngood,comp=badidx)
  140 +
  141 + if ngood gt 0 then begin
  142 + Kmatrix_fid[badidx]=!indef
  143 + Kmatrix[badidx]=!indef
  144 + Kmatrix_fid_red=Kmatrix_fid[goodidx]
  145 + Kmatrix_red=Kmatrix[goodidx]
  146 + diff_red=abs(Kmatrix_fid_red-Kmatrix_red)/Kmatrix_fid_red
  147 +
  148 + cc_diff_plo=percentile(diff_red,use_lo)
  149 + cc_diff_phi=percentile(diff_red,use_hi)
  150 + cc_diff_mean=la_mean(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
  151 + cc_diff_med=la_median(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
  152 + cc_diff_min=la_min(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
  153 + cc_diff_max=la_max(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
  154 + cc_diff_std=la_sigma(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
  155 +
  156 + if cc_diff_mean gt use_avg_thresh or cc_diff_med gt use_avg_thresh or cc_diff_max gt use_max_thresh then begin
  157 + problem_cc=[problem_cc,instrument+file_basename(fidcc_pl_files[i])]
  158 + flag_bad_cc=1
  159 + end
  160 +
  161 + print,'========================================================='
  162 + print,instrument,": Color Corrections for POWER LAWS from file ",file_basename(fidcc_pl_files[i])
  163 + print,"MIN / MAX / MEAN / MEDIAN / P25 / P75 / RMS"
  164 + print,cc_diff_min,cc_diff_max,cc_diff_mean,cc_diff_med,cc_diff_plo,cc_diff_phi,cc_diff_std
  165 + if flag_bad_cc ne 0 then begin
  166 + print,"WARNING: Discrepancy with fiducial CC values larger than avg/max thresholds of ",use_avg_thresh, use_max_thresh
  167 + print,"DustEMWrap:"
  168 + print,(tag_names(fidcc_pl_str))[1:-1]
  169 + FOR k=0L,Nv-1 DO print,Kmatrix[*,k]
  170 + print,"Published powerlaw CCs:"
  171 + print,(tag_names(fidcc_pl_str))[1:-1]
  172 + FOR k=0L,Nv-1 DO print,Kmatrix_fid[*,k]
  173 + end
  174 + if flag_bad_cc eq 0 then print,"CC values agree within avg/max thresholds of ",use_avg_thresh, use_max_thresh
  175 + print,'========================================================='
  176 + end else begin
  177 + message,'No matching filters / power law parameters found?',/info
  178 + stop
  179 + endelse
  180 +endfor
  181 +
  182 +;== Cross-check the BB CCs
  183 +bb_ccs:
  184 +for i=0,nbbfiles-1 do begin
  185 + flag_bad_cc=0
  186 + fidcc_bb_str=read_xcat(fidcc_bb_files[i],/silent)
  187 + tbb=fidcc_bb_str.tbb
  188 + Nv=n_elements(tbb)
  189 + Kmatrix=fltarr(Nfilt,Nv)
  190 + Kmatrix_fid=fltarr(Nfilt,Nv)
  191 + FOR j=0,Nfilt-1 DO Kmatrix_fid[j,*]=fidcc_bb_str.(j+1) ; tbb is the first tag (j=0)
  192 + FOR k=0L,Nv-1 DO BEGIN
  193 + FOR ff=0L,Nfilt-1 DO BEGIN
  194 + fwavs=*((*!dustem_filters).(tt).filter_wavelengths[ff])
  195 + spec=dustem_planck_function(tbb[k],fwavs)
  196 + sed=dustem_cc(fwavs,spec,filters[ff],cc=cc,fluxconv=fluxconv,help=help)
  197 + Kmatrix[ff,k]=cc
  198 + ENDFOR
  199 + ENDFOR
  200 +
  201 + goodidx=where(Kmatrix_fid ne !indef and Kmatrix ne !indef and Kmatrix_fid ne !values.f_nan and Kmatrix ne !values.f_nan,ngood,comp=badidx)
  202 +
  203 + if ngood gt 0 then begin
  204 + Kmatrix_fid[badidx]=!indef
  205 + Kmatrix[badidx]=!indef
  206 + Kmatrix_fid_red=Kmatrix_fid[goodidx]
  207 + Kmatrix_red=Kmatrix[goodidx]
  208 + diff_red=abs(Kmatrix_fid_red-Kmatrix_red)/Kmatrix_fid_red
  209 +
  210 + cc_diff_plo=percentile(diff_red,use_lo)
  211 + cc_diff_phi=percentile(diff_red,use_hi)
  212 + cc_diff_mean=la_mean(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
  213 + cc_diff_med=la_median(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
  214 + cc_diff_min=la_min(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
  215 + cc_diff_max=la_max(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
  216 + cc_diff_std=la_sigma(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
  217 + if cc_diff_mean gt use_avg_thresh or cc_diff_med gt use_avg_thresh or cc_diff_max gt use_max_thresh then begin
  218 + problem_cc=[problem_cc,instrument+file_basename(fidcc_pl_files[i])]
  219 + flag_bad_cc=1
  220 + end
  221 +
  222 + print,'========================================================='
  223 + print,instrument,": Color Corrections for BLACKBODIES from file ",file_basename(fidcc_bb_files[i])
  224 + print,"MIN / MAX / MEAN / MEDIAN / P25 / P75 / RMS"
  225 + print,cc_diff_min,cc_diff_max,cc_diff_mean,cc_diff_med,cc_diff_plo,cc_diff_phi,cc_diff_std
  226 + if flag_bad_cc ne 0 then begin
  227 + print,"WARNING: Discrepancy with fiducial CC values larger than avg/max thresholds of ",use_avg_thresh, use_max_thresh
  228 + print,"DustEMWrap:"
  229 + print,(tag_names(fidcc_bb_str))[1:-1]
  230 + FOR k=0L,Nv-1 DO print,Kmatrix[*,k]
  231 + print,"Published Blackbody CCs:"
  232 + print,(tag_names(fidcc_bb_str))[1:-1]
  233 + FOR k=0L,Nv-1 DO print,Kmatrix_fid[*,k]
  234 + end
  235 + if flag_bad_cc eq 0 then print,"CC values agree within avg/max thresholds of ",use_avg_thresh, use_max_thresh
  236 + print,'========================================================='
  237 + end else begin
  238 + message,'No matching filters / blackbody parameters found?',/info
  239 + stop
  240 + endelse
  241 +endfor
  242 +
  243 +goto,the_end
  244 +
  245 +;; ;== Cross-check the MBB CCs
  246 +;; ;== currently we have a problem in MBB function/formula?
  247 +;; mbb_ccs:
  248 +;; for i=0,nmbbfiles-1 do begin
  249 +;; flag_bad_cc=0
  250 +;; fidcc_mbb_str=read_xcat(fidcc_mbb_files[i],/silent)
  251 +;; tmbb=fidcc_mbb_str.tmbb
  252 +;; uniq_tmbb=tmbb[uniq(tmbb(sort(tmbb)))]
  253 +;; alpha=fidcc_mbb_str.alpha
  254 +;; uniq_alpha=alpha[uniq(alpha)]
  255 +;; Nt=n_elements(uniq_tmbb)
  256 +;; Na=n_elements(uniq_alpha)
  257 +;; Kmatrix=fltarr(Nfilt,Nt*Na)
  258 +;; Kmatrix_fid=fltarr(Nfilt,Nt*Na)
  259 +;; FOR j=0,Nfilt-1 DO Kmatrix_fid[j,*]=fidcc_mbb_str.(j+2) ; alpha and tmbb are the first tag (j=0,1)
  260 +;; FOR n=0L,Na-1 DO BEGIN
  261 +;; FOR k=0L,Nt-1 DO BEGIN
  262 +;; lambda0=100.
  263 +;; spec=dustem_mbb_function(tmbb[k],uniq_alpha[n],lambda0,wavs)
  264 +;; sed=dustem_cc(wavs,spec,filters,cc=cc,fluxconv=fluxconv,help=help)
  265 +;; startidx=n*nt+k & endidx=(n+1)*nt+k-1
  266 +;; thiscol=n*Nt+k
  267 +;; print,thiscol
  268 +;; Kmatrix[*,thiscol]=cc
  269 +;; ENDFOR
  270 +;; ENDFOR
  271 +
  272 +;; goodidx=where(Kmatrix_fid ne !indef and Kmatrix ne !indef,ngood,comp=badidx)
  273 +
  274 +;; if ngood gt 0 then begin
  275 +;; Kmatrix_fid[badidx]=!indef
  276 +;; Kmatrix[badidx]=!indef
  277 +;; cc_diff_mean=la_mean(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
  278 +;; cc_diff_med=la_median(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
  279 +;; cc_diff_min=la_min(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
  280 +;; cc_diff_max=la_max(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
  281 +;; cc_diff_std=la_sigma(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
  282 +;; if cc_diff_mean gt use_avg_thresh or cc_diff_med gt use_avg_thresh or cc_diff_max gt use_max_thresh then begin
  283 +;; problem_cc=[problem_cc,instrument+file_basename(fidcc_pl_files[i])]
  284 +;; flag_bad_cc=1
  285 +;; end
  286 +
  287 +;; print,'========================================================='
  288 +;; print,instrument,": Color Corrections for MODIFIED BLACKBODIES from file ",file_basename(fidcc_bb_files[i])
  289 +;; print,cc_diff_min,cc_diff_max,cc_diff_mean,cc_diff_med,cc_diff_std
  290 +;; if flag_bad_cc ne 0 then print,"WARNING: Discrepancy with fiducial CC values larger than avg/max thresholds of ",use_avg_thresh, use_max_thresh
  291 +;; if flag_bad_cc eq 0 then print,"CC values agree within avg/max thresholds of ",use_avg_thresh, use_max_thresh
  292 +;; print,'========================================================='
  293 +;; end else begin
  294 +;; message,'No matching filters / modified blackbody parameters found?',/info
  295 +;; stop
  296 +;; endelse
  297 +;; endfor
  298 +
  299 +stop
  300 +
  301 +;stop
  302 +
  303 +the_end:
  304 +END
... ...