PRO check_phangs_ssps_isrf_prediction,show_map=show_map,from_restore=from_restore,from_classes_restore=from_classes_restore,save=save

;check_phangs_ssps_isrf_prediction
;check_phangs_ssps_isrf_prediction,/from_restore
;check_phangs_ssps_isrf_prediction,/from_classes_restore,/save

;===== predicts the minimum ISRF in voronoi bins

win=0L
;window,win,xsize=900,ysize=1000 & win=win+1

;=== This is where the data is read from and the ISRFs will be stored
;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

IF keyword_set(from_restore) THEN goto,from_restore
IF keyword_set(from_classes_restore) THEN goto,from_classes_restore

source_name='ngc0628'
object_distance=9.84  ;MPc
object_thickness=0.1  ;kpc

dustem_define_la_common
obp=[1.1,0,1.15,1]

;dustem_init,show_plots=show_plots
;needed only for NHCO
restore,data_dir+'ngc0628_jwst_images.sav',/verb
;needed for stellar parameters and voronoi bin info
restore,data_dir+'ngc0628_muse_images.sav',/verb

Nvor=max(voronoi_id)

use_NHmap=NHCO       ;used NH map in 1.e21 (from CO)

;stop
!quiet=1

t1=systime(0,/sec)
first_vid=0L
lambir=dustem_get_wavelengths()
;show_each=10
show_each=100
Nlambir=n_elements(lambir)

ISRFS=fltarr(Nlambir,Nvor)
FOR vid=first_vid,Nvor-1 DO BEGIN
	IF vid mod show_each EQ 0 THEN BEGIN
		t2=systime(0,/sec)
		delta_t_hr=(t2-t1)/60.^2 ;elapsed time [hr]
		frac_perc=vid/Nvor*100
		delta_t_remain=delta_t_hr/frac_perc*100.-delta_t_hr  ;remaining time [hr]
		message,'done '+strtrim(frac_perc)+' % in '+strtrim(delta_t_hr,2)+' hr remaining '+strtrim(delta_t_remain,2)+' hr',/continue
	ENDIF
	;===== get the SSP weights for the given voronoi bins
	weights=phangs_binid2weights(st_muse_weights,vid,st_templates,age_values,metalicity_values,reddening=reddening)
	amplitude=1.
	reddening=0.     ;This sets Muse reddening to 0. Will be contsrained later
	fpd=phangs_stellar_continuum_plugin_weight2params(weights,parameter_values=val,redenning=reddening,/force_include_reddening,amplitude=amplitude,/force_include_amplitude)
   Nparams=n_elements(val)
   key=intarr(Nparams)
   FOR i=0L,Nparams-1 DO BEGIN
	    toto=dustem_parameter_description2type(fpd[i],string_name=string_name,key=one_key)
	    key[i]=one_key
 	ENDFOR
 	ISRFS[*,vid]=dustem_plugin_phangs_stellar_isrf(key=key,val=val,object_distance=object_distance,object_thickness=object_thickness)
ENDFOR

save,ISRFS,object_distance,object_thickness,source_name,file=data_dir+'ngc0628_isrf_min_prediction.sav',/verb

from_restore:

;make_phangs_isrf_classes,/save
make_phangs_isrf_classes

stop

from_classes_restore:

;file=data_dir+'isrf_classes_one_ratio.sav'
restore,data_dir+'isrf_classes_one_ratio.sav',/verbose
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.

Nclasses=n_elements(classes)
refs_wavs=[classes[0].FLUX_RATIO_WAV]
reference_wavelength=classes[0].reference_wavelength
ratios=classes.flux_ratio

lambir=dustem_get_wavelengths()

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

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

;compute flux ratio for each voronoi bin
vor_ratios=fltarr(Nvor)
FOR i=0L,Nvor-1 DO BEGIN
	vor_ratios[i]=interpol(n_ISRFs[*,i],lambir,refs_wavs[0])
ENDFOR
;fill in template ISRF of each class
;method='one_spectrum'
method='average_spectrum'
empty_classes=intarr(Nclasses)
CASE method OF
	'one_spectrum':BEGIN ;This is selection based on a single ISRF
		FOR i=0L,Nclasses-1 DO BEGIN
			ind=where(vor_ratios GT ratios[i],count)
			classes[i].ISRF_template=ptr_new(n_ISRFS[*,ind[0]])
		ENDFOR
	END
	'average_spectrum':BEGIN ;This is based on an average
		delta=1.0 ;This is a bit arbitrary. Chosen so that there is no empty classes
		FOR i=0L,Nclasses-1 DO BEGIN
			;ind=where(vor_ratios GE ratios[i]-ratios[i]*delta and vor_ratios LE ratios[i]+ratios[i]*delta,count)
			ratio_min=classes[i].flux_ratio_min
			ratio_max=classes[i].flux_ratio_max
			;ind=where(vor_ratios GT ratios_min[i] and vor_ratios LE ratios_max[i],count)
			ind=where(vor_ratios GT ratio_min and vor_ratios LE ratio_max,count)
			IF count NE 0 THEN BEGIN
				classes[i].ISRF_template=ptr_new(la_mean(n_ISRFS[*,ind],dim=-1))
			ENDIF ELSE BEGIN
				message,'Found no voronoi bin in class '+strtrim(i+1,2),/continue
				;stop
				empty_classes[i]=1
			ENDELSE
		ENDFOR
	END
ENDCASE
;===== remove empty classes
ind=where(empty_classes EQ 0,count)
classes=classes[ind]
Nclasses=count
;ratios_min=ratios_min[ind]
;ratios_max=ratios_max[ind]
;===== Fill in class number and name
FOR i=0L,Nclasses-1 DO BEGIN
	classes[i].number=i+1       ;because class 0 will be for Mathis ISRF
	str='I('+strtrim(refs_wavs[0],2)+'/I('+strtrim(reference_wavelength,2)+')='+strtrim(ratios[i],2)
	classes[i].name=str
ENDFOR

stop

;classify Voronoi bin according to flux ratio
;This could also be done on chi2, instead of just ratio ...
vor_class=lonarr(Nvor)
FOR i=0L,Nvor-1 DO BEGIN
  ;ratio=interpol(n_ISRFs[*,i],lambir,refs_wavs[0])
  this_ratio=vor_ratios[i]
  ind=where(classes.flux_ratio_min LT this_ratio and classes.flux_ratio_max GE this_ratio,count)
  ;ind=where(ratios_min LT this_ratio and ratios_max GE this_ratio,count)
  xv=ind[0]
  ;xv=round(interpol(xratios,ratios,this_ratio))
  IF ptr_valid(classes[xv].voronoi_members) THEN BEGIN
     classes[xv].voronoi_members=ptr_new([*classes[xv].voronoi_members,i])
  ENDIF ELSE BEGIN
  	  classes[xv].voronoi_members=ptr_new([i])
  ENDELSE
  classes[xv].Nvor=classes[xv].Nvor+1
  this_chi2=total((n_ISRFS[*,i])^2-(*classes[xv].ISRF_template)^2)
  IF ptr_valid(classes[xv].voronoi_chi2) THEN BEGIN
  	 classes[xv].voronoi_chi2=ptr_new([*classes[xv].voronoi_chi2,this_chi2])
  ENDIF ELSE BEGIN
  	 classes[xv].voronoi_chi2=ptr_new([this_chi2])
  ENDELSE
  ;This is a voronoi vector of the belonging class
  vor_class[i]=xv
ENDFOR
;classes[xv].voronoi_chi2=ptr_new((*classes[xv].voronoi_chi2)[1:*])
;classes[xv].voronoi_members=ptr_new((*classes[xv].voronoi_members)[1:*])

;==== plot histogram of voronoi classes
res=histogram(vor_class,min=0,max=Nclasses-1,Nbin=Nclasses,locations=x)

window,win & win=win+1
!p.multi=[0,1,2]
cgplot,x,res,psym=10,xrange=[-1,Nclasses+1],/xstyle,xtit='Class',ytit='Npix'
cgplot,ratios,res,psym=10,xrange=[1.e-5,10],/xstyle,xtit='Ratio',ytit='Npix',/xlog

;==== make a sky map of the classes
window,win & win=win+1
map_classes=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'))+la_undef()
FOR i=0L,Nvor-1 DO BEGIN
    map_classes[*all_seds_indices[i]]=vor_class[i]
ENDFOR
obp=[1.1,0.,1.15,1]
image_cont20,map_classes,href,/square,/silent,image_color_table='jpbloadct',axis_color_table=1,imrange=[-1,Nclasses+1],off_bar=obp,tit='Class map'

;=== plot ISRF templates of all classes
window,win & win=win+1
xtit='Wavelength [mic]'
ytit='Normalized ISRF'
tit='ISRF classes'
loadct,13
!p.multi=0
med_isrf=la_median(n_ISRFS,dim=-1)
cgplot,lambir,med_isrf,color='black',/xlog,/ylog,xrange=[1.e-2,1.e2],/xsty,yrange=[1.e-5,10.],/nodata,xtit=xtit,ytit=ytit,tit=tit
colors=fix(range_gen(Nclasses,[10,250]))
FOR i=0L,Nclasses-1 DO cgoplot,lambir,*classes[i].ISRF_template,color=colors[i]
cgoplot,lambir,med_isrf,color='black',thick=2.
cgoplot,lambir,*class0.isrf_template,color='black',thick=2.
cgoplot,replicate(refs_wavs,2),10^!y.crange,linestyle=2
cgoplot,replicate(reference_wavelength,2),10^!y.crange,linestyle=2


;stop

IF keyword_set(save) THEN BEGIN
	file_save=data_dir+'ngc0682_isrf_classes_one_ratio.sav'
	message,'About to write '+file_save,/continue
	message,'Make sure you are allowed to do that and type .c to proceed.',/continue
	save,classes,vor_class,map_classes,file=file_save,/verb
	message,'Saved '+file_save,/continue
ENDIF

message,'========== Finished computing ISRF classes',/continue

stop


END