PRO make_phangs_isrf_classes,save=save,help=help

;+
; NAME:
;       make_phangs_isrf_classes
; CALLING SEQUENCE:
;       make_phangs_isrf_classes[,/save]
; PURPOSE:
;       generates the ISRF classes used in fiting PHANGs data
; INPUTS:
;       None
; OPTIONAL KEYWORDS:
;       save        = if set, save the classes
;       help        = if set, print this help
; OUTPUTS:
;	    None
; OPTIONAL INPUT:
;       None
; OPTIONAL OUTPUT:
;       None
; PROCEDURE AND SUBROUTINE USED
;       None
; SIDE EFFECTS:
;       None
; EXAMPLE:
;       make_phangs_isrf_classes,/save
; MODIFICATION HISTORY:
;       written by Jean-Philippe Bernard
;-

IF keyword_set(help) THEN BEGIN
  doc_library,'make_phangs_isrf_classes'
  goto,the_end
ENDIF

win=0L

;data_dir='/Volumes/PILOT_FLIGHT1/PHANGS-JWST/DR1/'
data_dir=!phangs_data_dir+'/ISRF/WORK/'

use_model='DBP90'    ;Example with default keywords uses the DBP90 model
use_polarization=0   ; initialize Dustemwrap in no polarization mode 

;== INITIALISE DUSTEM
dustem_init,model=use_model,polarization=use_polarization,show_plots=show_plots

;==== restore needed data
restore,data_dir+'ngc0628_isrf_min_prediction.sav',/verb
;% RESTORE: Restored variable: ISRFS.
;% RESTORE: Restored variable: OBJECT_DISTANCE.
;% RESTORE: Restored variable: OBJECT_THICKNESS.
;% RESTORE: Restored variable: SOURCE_NAME.
restore,data_dir+'ngc0628_astrosat_voronoi_prediction_fast.sav',/verb
;% RESTORE: Restored variable: ACTUAL_SEDS.
;% RESTORE: Restored variable: PREDICTED_SEDS.
;% RESTORE: Restored variable: HREF.
;% RESTORE: Restored variable: ALL_SEDS_INDICES.
;% RESTORE: Restored variable: DO_NORMALIZE.
;% RESTORE: Restored variable: ALL_FILTERS.
;% RESTORE: Restored variable: NHVOR.

;stop

;=== attempt to classify ISRFs

;normalized ISRFs
lambir=dustem_get_wavelengths(isrf_wavelengths=isrf_wavelengths)

Nlambir=(size(ISRFS))[1]
Nvor=(size(ISRFS))[2]

n_ISRFS=fltarr(Nlambir,Nvor)
n_wavs=fltarr(Nlambir,Nvor)
ind=where(isrf_wavelengths GT 1.)
indd=ind[0]
FOR vid=0L,Nvor-1 DO BEGIN
  n_ISRFS[*,vid]=ISRFS[*,vid]/ISRFS[indd,vid]
  n_wavs[*,vid]=isrf_wavelengths
ENDFOR

;=== 2D histogram of normalized ISRFs in NGC0628
x=reform(n_wavs,Nlambir*Nvor)
y=reform(n_ISRFS,Nlambir*Nvor)

;jpbloadct
loadct,13
Nx=200 & Ny=200
Nx=100 & Ny=100
xmin=-1.5 & xmax=2.
ymin=-5 & ymax=1.
xtit='log10(wav [mic])'
ytit='log10(ISRF/ISRF(1 mic))'
range=[-1,5]

window,win & win=win+1
plot_bin_correl,alog10(x),alog10(y),Nx,Ny,xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,xtit=xtit,ytit=ytit,range=range;,/xlog ;,/ylog

;compute the median ISRFS
med_isrf=la_median(n_ISRFS,dim=-1)

;refs_wavs=[0.11]      ;This is set to be just before the Lymann break
reference_wavelength=1.  ;all ISRFs are normalized to 1 mic
refs_wavs=[0.15]      ;This is set to be just before the Lymann break
ref_ratios=interpol(med_isrf,isrf_wavelengths,refs_wavs)
;These are the ratios which will be used to classify ISRFs

one_class={number:0L $;class number
           ,name:'' $    ;class name
           ,ISRF_template:ptr_new() $  ;isrf template
           ,flux_ratio:0. $ ;flux ratio
           ,flux_ratio_wav:0. $ ;flux ratio wavelength
           ,reference_wavelength:reference_wavelength $ ;ISRF normalisation wavelength
           ,flux_ratio_min:0. $ ;minimum flux ratio to belong to this class
           ,flux_ratio_max:0. $ ;maximum flux ratio to belong to this class           
           ,flux_ref:0. $ ;
           ,Nvor:0L  $   ;number of Voronoi bins in the class
           ,voronoi_members:ptr_new() $ ;voronoi pixel members of that class
           ,voronoi_chi2:ptr_new()}    ;chi2 between ref ISRF and voronoi member's ISRF

;class_ratio_range=[8.e-5,3.]
;class_ratio_range=[5.e-5,3.]
;Nclasses=20L
class_ratio_range=[1.e-5,10.]
Nclasses=30L
;Nclasses=200L  ;for test
ratios=range_gen(Nclasses,class_ratio_range,/log,/get_minmax,min_values=ratios_min,max_values=ratios_max)
xratios=indgen(Nclasses)
classes=replicate(one_class,Nclasses)

;fill in class values
classes.number=lindgen(Nclasses)
classes.flux_ratio=ratios
classes.flux_ratio_wav=refs_wavs[0]
classes.flux_ratio_min=ratios_min
classes.flux_ratio_max=ratios_max
FOR i=0L,Nclasses-1 DO BEGIN
	classes[i].ISRF_template=ptr_new(n_ISRFs[*,i])
ENDFOR

;=== Do a class (class 0) with a Mathis field
class0=one_class
class0.number=0L
class0.name='Mathis'
class0.reference_wavelength=1.
class0.flux_ratio_wav=refs_wavs[0]
file=!dustem_soft_dir+'/data/ISRF_MATHIS.DAT'
;file='/Users/jpb/Soft_Libraries/dustem_fortran/data/ISRF.DAT'
readcol,file,Mathis_wavs,Mathis_ISRF
Mathis_ISRF=interpol(Mathis_ISRF,Mathis_wavs,isrf_wavelengths)
n_Mathis_ISRF=Mathis_ISRF/interpol(Mathis_ISRF,isrf_wavelengths,class0.reference_wavelength)
class0.ISRF_template=ptr_new(n_Mathis_ISRF)
ratio=interpol(n_Mathis_ISRF,isrf_wavelengths,class0.flux_ratio_wav)
class0.flux_ratio=ratio
class0.flux_ratio_min=ratio*(1.-0.1)  ;arbitrary +-10% around ratio
class0.flux_ratio_max=ratio*(1.+0.1)  ;arbitrary +-10% around ratio
Mathis_1mic=interpol(Mathis_ISRF,isrf_wavelengths,1.0)
message,'Mathis (G0=1) has 4*pi*Inu='+strtrim(Mathis_1mic,2)+' ergs/s/cm2/Hz',/continue
message,'Mathis (G0=1) has Inu(0.15 mic)/Inu(1 mic)='+strtrim(ratio,2)+' ',/continue

;stop

;save empty classes defined above
IF keyword_set(save) THEN BEGIN
	file_save=data_dir+'isrf_classes_one_ratio.sav'
	message,'Will save '+file_save,/continue
	message,'.c to proceed ',/continue
	stop
	save,classes,class0,file=file_save,/verb
	message,'Saved '+file_save,/continue
ENDIF

the_end:

END