PRO make_phangs_isrf_classes,bidon,source_name=source_name,fits=fits,help=help,resolution_filter=resolution_filter

;+
; 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:

;       source_name   = source name. Default = ngc0628
;       fits          = if set, saves the classes into a .fits file
;       help          = if set, print this help

; OUTPUTS:
;	    None
; OPTIONAL INPUT:
;       None
; OPTIONAL OUTPUT:
;       None
; PROCEDURE AND SUBROUTINE USED

;       Uses the following files:
;       *_isrf_min_prediction.sav
;       *_ref_header.sav
;       *_astrosat_seds_muse_pixels.sav

; SIDE EFFECTS:

;       Produces the following files:
;       isrf_classes_one_ratio_on_*.fits
;       *_isrf_classes_map.fits
;       *_isrf_classes_voronoi.sav

; 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

use_source_name='ngc0628'
IF keyword_set(source_name) THEN use_source_name=source_name

reso_str=''
IF keyword_set(resolution_filter) THEN BEGIN
    reso_str='_'+resolution_filter
ENDIF

;==== restore needed data

;CAUTION: Here ISRFs has a different shape wrt voronoi bins when reso_str is set vs not set ...
restore,data_dir+use_source_name+'_isrf_min_prediction'+reso_str+'.sav',/verb

;% RESTORE: Restored variable: ISRFS.
;% RESTORE: Restored variable: OBJECT_DISTANCE.
;% RESTORE: Restored variable: OBJECT_THICKNESS.
;% RESTORE: Restored variable: SOURCE_NAME.

;=== The following is restored only for variables HREF and ALL_SEDS_INDICES
;=== in order to generate the fits class map fo that galaxy

restore,data_dir+use_source_name+'_ref_header'+reso_str+'.sav',/verb

;% RESTORE: Restored variable: HREF.

;=== The following is restored only for variables HREF and ALL_SEDS_INDICES
;=== in order to generate the fits class map fo that galaxy

restore,data_dir+use_source_name+'_seds_indices'+reso_str+'.sav',/verb
;% RESTORE: Restored variable: ALL_SEDS_INDICES.
;restore,data_dir+use_source_name+'_astrosat_seds_muse_pixels.sav',/verb

;% RESTORE: Restored variable: ALL_SEDS.

;% RESTORE: Restored variable: ALL_SEDS_INDICES.

class_map=lonarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))

;stop

;=== Defines ISRF classes, based on a given galaxy

lambir=dustem_get_wavelengths(isrf_wavelengths=isrf_wavelengths)

Nlamb=(size(ISRFS))[1]

Nvor=(size(ISRFS))[2]

;Define normalized ISRFs
n_ISRFS=fltarr(Nlamb,Nvor)     ;ISRFs normalized to 1 micron
n_wavs=fltarr(Nlamb,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

;=== Plot 2D histogram of normalized ISRFs in NGC0628
x=reform(n_wavs,Nlamb*Nvor)
y=reform(n_ISRFS,Nlamb*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)  ;actually not used

;These are the ratios which will be used to classify ISRFs

;==== compute ISRF ration for each Voronoi bin
isrf_ratio_values=fltarr(Nvor)

vor_class=lonarr(Nvor)         ;This is the class of each voronoi bin

FOR i=0L,Nvor-1 DO BEGIN
  isrf_ratio_values[i]=interpol(n_ISRFS[*,i],isrf_wavelengths,refs_wavs[0])
ENDFOR

undef_array=fltarr(Nlamb)
undef_array[*]=la_undef()

one_class={number:0L $                                  ;class number
           ,name:'' $                                   ;class name
           ,ISRF_wavelengths:undef_array $              ;isrf wavelengths
           ,ISRF_template:undef_array $                 ;isrf template, normalized at reference_wavelength
           ,flux_ratio:0. $                             ;flux ratio between flux_ratio_wav and reference_wavelength
           ,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=[1.e-5,10.]    ;reange of flux ratios used in the classification
Nclasses=30L                     ;Number of classes used for classification

ratios=range_gen(Nclasses,class_ratio_range,/log,/get_minmax,min_values=ratios_min,max_values=ratios_max)

classes=replicate(one_class,Nclasses)

;====== fill in class values
Nvor_class=lonarr(Nclasses)

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

median_isrf=fltarr(Nlamb)

FOR i=0L,Nclasses-1 DO BEGIN

  ind=where(isrf_ratio_values GE ratios_min[i] AND isrf_ratio_values LT ratios_max[i],count)
  IF count NE 0 THEN  BEGIN

    vor_class[ind]=i             ;This is the class of each voronoi bin

    FOR j=0L,Nlamb-1 DO BEGIN
      median_isrf[j]=la_median(n_ISRFs[j,ind])
    ENDFOR
  	classes[i].ISRF_template=median_isrf
    classes[i].ISRF_wavelengths=isrf_wavelengths
    Nvor_class[i]=count
    FOR j=0L,count-1 DO BEGIN

      IF ptr_valid(ALL_SEDS_INDICES[ind[j]]) THEN BEGIN
        class_map[*ALL_SEDS_INDICES[ind[j]]]=i
      ENDIF

    ENDFOR
  ENDIF

ENDFOR

;stop

;=== 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)
class0.ISRF_template=n_Mathis_ISRF
class0.ISRF_wavelengths=isrf_wavelengths

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

;=== Do a plot of ISRF classes
window,win & win=win+1
colors=range_gen(Nclasses,[0,255])
yrange=[1.e-5,10]
xrange=[1.e-2,1.e4]
xtit='Wavlenghts [mic]'
ytit='Normalized ISRF (Inu)'
cgplot,class0.ISRF_wavelengths,class0.ISRF_template,/xlog,/ylog,yrange=yrange,xrange=xrange,xtit=xtit,ytit=ytit
FOR i=0L,Nclasses-1 DO BEGIN
  oplot,classes[i].ISRF_wavelengths,classes[i].ISRF_template,color=colors[i]
ENDFOR

;=== plot an histogram of the classes:
window,win & win=win+1
xtit='ISRF class #'
ytit='Number of Voronoi bins'
cgplot,Nvor_class,psym=10,xtit=xtit,ytit=ytit

;=== plot the class map:
cleanplot
window,win & win=win+1
obp=[1.1,0.,1.15,1]
imrange=[0,30]
image_cont20,class_map,href,/square,/silent,image_color_table='jpbloadct',axis_color_table=1,off_bar=obp,imrange=imrange

;stop

;=== save a fits file
IF keyword_set(fits) THEN BEGIN

  file_save=data_dir+'isrf_classes_one_ratio_on_'+use_source_name+reso_str+'.fits'

  ;message,'Will save '+file_save,/continue
  ;message,'.c to proceed ',/continue
  ;stop
  mwrfits,class0,file_save,/create
  mwrfits,classes,file_save

  message,'Saved '+file_save,/continue

  file_save=data_dir+use_source_name+'_isrf_classes_map'+reso_str+'.fits'

  writefits,file_save,class_map,href

  message,'Saved '+file_save,/continue

  file_save=data_dir+use_source_name+'_isrf_classes_voronoi'+reso_str+'.sav'

  save,vor_class,file=file_save

  message,'Saved '+file_save,/continue

ENDIF

the_end:

END