PRO dustem_cc_verify,instrument $
                     ,nwav=nwav $
                     ,wavrng=wavrng $
                     ,logspace=logspace $
                     ,avg_thresh=avg_thresh $
                     ,max_thresh=max_thresh $
                     ,bb=bb,mbb=mbb,power=power $
                     ,help=help

;+
; NAME:
;    dustem_cc_verify
;
; PURPOSE:
;    Compares color corrections for a given spectrum in a given filter
;    to fiducial values (in DATA/FIDUCIAL_CC)
;    This is a sanity-checking routine.
;    Output should be compared to the color-correction tables published by different
;    instrument teams.
;    e.g. 1  MIPS HandBook, page 98
;    http://irsa.ipac.caltech.edu/data/SPITZER/docs/mips/mipsinstrumenthandbook/MIPS_Instrument_Handbook.pdf
;    http://irsa.ipac.caltech.edu/data/SPITZER/docs/mips/mipsinstrumenthandbook/51/
;    AND
;    e.g. 2  IRAS Explanatory Supplement, Table VI.C.6
;    https://irsa.ipac.caltech.edu/IRASdocs/exp.sup/ch6/tabsupC6.html  
;
; CATEGORY:
;    DustEMWrap, Distributed, Low-Level, User convenience
;
; CALLING SEQUENCE:
;    dustem_cc_verify,instrument[,help=help]
;
; INPUTS:
;    instrument: string, instrument which color corrections will be compared
;                (currently only accepts 1 instrument)
;
; OPTIONAL INPUT PARAMETERS:
;    None
;
; OUTPUTS:
;
; OPTIONAL OUTPUT PARAMETERS:
;
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
;
; COMMON BLOCKS:
;    None
;
; SIDE EFFECTS:
;    None
;
; RESTRICTIONS:
;    The DustEM fortran code must be installed
;    The DustEMWrap IDL code must be installed
;    If restart is not set, this routine should be called after the
;    successful completion of a DustEMWrap run (i.e. without quitting the
;    IDL session) since it will use existing system variables 
;
; PROCEDURES AND SUBROUTINES USED:
;    
; EXAMPLES
;    dustem_cc_verify,'WISE'
;
; MODIFICATION HISTORY:
;    Written by AH Jan 2023
;    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_cc_verify'
  goto,the_end
ENDIF

use_avg_thresh=0.01             ; a mean difference of >=1% between our CCs and the tabulated values will be recorded
use_max_thresh=0.02             ; a max difference of >=2% between our CCs and the tabulated values will be recorded
if keyword_set(avg_thresh) then use_avg_thresh=avg_thresh
if keyword_set(max_thresh) then use_max_thresh=max_thresh
use_lo=25
use_hi=75

;== Initialize DustemWrap
dustem_init
use_Nwav=500000 & use_wavmin=0.01 & use_wavmax=5000.
if keyword_set(Nwav) then use_Nwav=Nwav

if keyword_set(wavrng) then begin
   use_wavmin=wavrng[0]
   use_wavmax=wavrng[1]
endif
wavs=findgen(use_Nwav)/(use_Nwav-1)*(use_wavmax-use_wavmin)+use_wavmin
if keyword_set(logspace) then wavs=10^(alog10(use_wavmin)+alog10(use_wavmax)*(dindgen(use_Nwav+1)/use_Nwav))
;restart=0
;end else begin
;   st=dustem_run()
;   wavs=st.sed.wav
;end

;== Check user gave us an instrument
if n_elements(instrument) ne 1 then begin
   message,'You must provide one instrument name',/info
   stop
end

;== Get associated filters and central wavelengths
filters=dustem_instru2filters(instrument)
Nfilt=n_elements(filters)
wavref=dustem_filter2wav(filters)

;== potential fiducial CC filenames
fidcc_pl_files=file_search(!dustem_wrap_soft_dir+'/Data/FIDUCIAL_CC/'+strlowcase(instrument)+'*_powerlaw.xcat',count=nplfiles)
fidcc_bb_files=file_search(!dustem_wrap_soft_dir+'/Data/FIDUCIAL_CC/'+strlowcase(instrument)+'*_bb.xcat',count=nbbfiles)
; problem in MBB formula?
;fidcc_mbb_files=file_search(!dustem_wrap_soft_dir+'/Data/FIDUCIAL_CC/'+strlowcase(instrument)+'*_mbb.xcat',count=nmbbfiles)

problem_cc=['NONE']

tt=where(tag_names((*!dustem_filters)) eq instrument)
;== Cross-check the power law CCs

pl_ccs:
for i=0,nplfiles-1 do begin
   flag_bad_cc=0
   fidcc_pl_str=read_xcat(fidcc_pl_files[i],/silent)
   pows=fidcc_pl_str.alpha
   Nv=n_elements(pows)
   Kmatrix=fltarr(Nfilt,Nv)
   Kmatrix_fid=fltarr(Nfilt,Nv)
   FOR j=0,Nfilt-1 DO Kmatrix_fid[j,*]=fidcc_pl_str.(j+1) ; alpha is the first tag (j=0)
   FOR k=0L,Nv-1 DO BEGIN
      FOR ff=0L,Nfilt-1 DO BEGIN
         fwavs=*((*!dustem_filters).(tt).filter_wavelengths[ff])
         spec=fwavs^(-1.*pows[k])
         sed=dustem_cc(fwavs,spec,filters[ff],cc=cc,fluxconv=fluxconv,help=help)
         Kmatrix[ff,k]=cc
      ENDFOR
   ENDFOR

   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)

   if ngood gt 0 then begin
      Kmatrix_fid[badidx]=!indef
      Kmatrix[badidx]=!indef
      Kmatrix_fid_red=Kmatrix_fid[goodidx]
      Kmatrix_red=Kmatrix[goodidx]
      diff_red=abs(Kmatrix_fid_red-Kmatrix_red)/Kmatrix_fid_red
      
      cc_diff_plo=percentile(diff_red,use_lo)
      cc_diff_phi=percentile(diff_red,use_hi)
      cc_diff_mean=la_mean(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
      cc_diff_med=la_median(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
      cc_diff_min=la_min(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
      cc_diff_max=la_max(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
      cc_diff_std=la_sigma(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))

      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
         problem_cc=[problem_cc,instrument+file_basename(fidcc_pl_files[i])]
         flag_bad_cc=1
      end
      
      print,'========================================================='
      print,instrument,": Color Corrections for POWER LAWS from file ",file_basename(fidcc_pl_files[i])
      print,"MIN / MAX / MEAN / MEDIAN / P25 / P75 / RMS"
      print,cc_diff_min,cc_diff_max,cc_diff_mean,cc_diff_med,cc_diff_plo,cc_diff_phi,cc_diff_std
      if flag_bad_cc ne 0 then begin
         print,"WARNING: Discrepancy with fiducial CC values larger than avg/max thresholds of ",use_avg_thresh, use_max_thresh
         print,"DustEMWrap:"
         print,(tag_names(fidcc_pl_str))[1:-1]
         FOR k=0L,Nv-1 DO print,Kmatrix[*,k]
         print,"Published powerlaw CCs:"
         print,(tag_names(fidcc_pl_str))[1:-1]
         FOR k=0L,Nv-1 DO print,Kmatrix_fid[*,k]
      end
      if flag_bad_cc eq 0 then print,"CC values agree within avg/max thresholds of ",use_avg_thresh, use_max_thresh
      print,'========================================================='
   end else begin
      message,'No matching filters / power law parameters found?',/info
      stop
   endelse
endfor

;== Cross-check the BB CCs 
bb_ccs:
for i=0,nbbfiles-1 do begin
   flag_bad_cc=0
   fidcc_bb_str=read_xcat(fidcc_bb_files[i],/silent)
   tbb=fidcc_bb_str.tbb
   Nv=n_elements(tbb)
   Kmatrix=fltarr(Nfilt,Nv)
   Kmatrix_fid=fltarr(Nfilt,Nv)
   FOR j=0,Nfilt-1 DO Kmatrix_fid[j,*]=fidcc_bb_str.(j+1) ; tbb is the first tag (j=0)
   FOR k=0L,Nv-1 DO BEGIN
      FOR ff=0L,Nfilt-1 DO BEGIN
         fwavs=*((*!dustem_filters).(tt).filter_wavelengths[ff])
         spec=dustem_planck_function(tbb[k],fwavs)
         sed=dustem_cc(fwavs,spec,filters[ff],cc=cc,fluxconv=fluxconv,help=help)
         Kmatrix[ff,k]=cc
      ENDFOR
   ENDFOR

   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)

   if ngood gt 0 then begin
      Kmatrix_fid[badidx]=!indef
      Kmatrix[badidx]=!indef
      Kmatrix_fid_red=Kmatrix_fid[goodidx]
      Kmatrix_red=Kmatrix[goodidx]
      diff_red=abs(Kmatrix_fid_red-Kmatrix_red)/Kmatrix_fid_red

      cc_diff_plo=percentile(diff_red,use_lo)
      cc_diff_phi=percentile(diff_red,use_hi)
      cc_diff_mean=la_mean(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
      cc_diff_med=la_median(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
      cc_diff_min=la_min(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
      cc_diff_max=la_max(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
      cc_diff_std=la_sigma(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
      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
         problem_cc=[problem_cc,instrument+file_basename(fidcc_pl_files[i])]
         flag_bad_cc=1
      end
      
      print,'========================================================='
      print,instrument,": Color Corrections for BLACKBODIES from file ",file_basename(fidcc_bb_files[i])
      print,"MIN / MAX / MEAN / MEDIAN / P25 / P75 / RMS"
      print,cc_diff_min,cc_diff_max,cc_diff_mean,cc_diff_med,cc_diff_plo,cc_diff_phi,cc_diff_std
      if flag_bad_cc ne 0 then begin
         print,"WARNING: Discrepancy with fiducial CC values larger than avg/max thresholds of ",use_avg_thresh, use_max_thresh
         print,"DustEMWrap:"
         print,(tag_names(fidcc_bb_str))[1:-1]
         FOR k=0L,Nv-1 DO print,Kmatrix[*,k]
         print,"Published Blackbody CCs:"
         print,(tag_names(fidcc_bb_str))[1:-1]
         FOR k=0L,Nv-1 DO print,Kmatrix_fid[*,k]
      end
      if flag_bad_cc eq 0 then print,"CC values agree within avg/max thresholds of ",use_avg_thresh, use_max_thresh
      print,'========================================================='
   end else begin
      message,'No matching filters / blackbody parameters found?',/info
      stop
   endelse
endfor

goto,the_end

;; ;== Cross-check the MBB CCs 
;; ;== currently we have a problem in MBB function/formula?
;; mbb_ccs:
;; for i=0,nmbbfiles-1 do begin
;;    flag_bad_cc=0
;;    fidcc_mbb_str=read_xcat(fidcc_mbb_files[i],/silent)
;;    tmbb=fidcc_mbb_str.tmbb
;;    uniq_tmbb=tmbb[uniq(tmbb(sort(tmbb)))]
;;    alpha=fidcc_mbb_str.alpha
;;    uniq_alpha=alpha[uniq(alpha)]
;;    Nt=n_elements(uniq_tmbb)
;;    Na=n_elements(uniq_alpha)
;;    Kmatrix=fltarr(Nfilt,Nt*Na)
;;    Kmatrix_fid=fltarr(Nfilt,Nt*Na)
;;    FOR j=0,Nfilt-1 DO Kmatrix_fid[j,*]=fidcc_mbb_str.(j+2) ; alpha and tmbb are the first tag (j=0,1)
;;    FOR n=0L,Na-1 DO BEGIN
;;       FOR k=0L,Nt-1 DO BEGIN
;;          lambda0=100.
;;          spec=dustem_mbb_function(tmbb[k],uniq_alpha[n],lambda0,wavs)
;;           sed=dustem_cc(wavs,spec,filters,cc=cc,fluxconv=fluxconv,help=help)
;;          startidx=n*nt+k & endidx=(n+1)*nt+k-1
;;          thiscol=n*Nt+k
;;          print,thiscol
;;          Kmatrix[*,thiscol]=cc
;;       ENDFOR
;;    ENDFOR
   
;;    goodidx=where(Kmatrix_fid ne !indef and Kmatrix ne !indef,ngood,comp=badidx)

;;    if ngood gt 0 then begin
;;       Kmatrix_fid[badidx]=!indef
;;       Kmatrix[badidx]=!indef
;;       cc_diff_mean=la_mean(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
;;       cc_diff_med=la_median(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
;;       cc_diff_min=la_min(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
;;       cc_diff_max=la_max(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
;;       cc_diff_std=la_sigma(la_div((la_abs(la_sub(Kmatrix_fid,Kmatrix))),Kmatrix_fid))
;;       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
;;          problem_cc=[problem_cc,instrument+file_basename(fidcc_pl_files[i])]
;;          flag_bad_cc=1
;;       end
      
;;       print,'========================================================='
;;       print,instrument,": Color Corrections for MODIFIED BLACKBODIES from file ",file_basename(fidcc_bb_files[i])
;;       print,cc_diff_min,cc_diff_max,cc_diff_mean,cc_diff_med,cc_diff_std
;;       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
;;       if flag_bad_cc eq 0 then print,"CC values agree within avg/max thresholds of ",use_avg_thresh, use_max_thresh
;;       print,'========================================================='
;;    end else begin
;;       message,'No matching filters / modified blackbody parameters found?',/info
;;       stop
;;    endelse
;; endfor

stop

;stop

the_end:
END