PRO phangs_smooth_muse_isrf,source_name $
					       ,resolution_filter=resolution_filter $
					       ,save=save $
					       ,show_plots=show_plots $
					       ,nostop=nostop $
					       ,test=test

;+
; NAME:
;    phangs_smooth_muse_isrf
; PURPOSE:
;    Smoothes ISRF to a given resolution
; CATEGORY:
;    PHANGS ISRF
; CALLING SEQUENCE:
;    phangs_smooth_muse_isrf,source_name[,resolution_filter=][,/save][,/test]=
; INPUTS:
;    source_name            : source name
; OPTIONAL INPUT PARAMETERS:
;    resolution_filter      :resolution of the filter to be used for smoothing (defaults=)
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    save      = If set, saves results
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    The following files are written, if they do not already exist:
;    _isrf_min_prediction_*.sav
; RESTRICTIONS:
;    None
; PROCEDURE:
;    
; EXAMPLES
;	phangs_smooth_muse_isrf,'ngc0628',reso_filter='SPIRE3'
; MODIFICATION HISTORY:
;    Written by J.-Ph. Bernard (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,'phangs_smooth_muse_isrf'
  goto,the_end
ENDIF

dustem_define_la_common
dustem_init
obp=[1.1,0,1.15,1]
imrange=[-0.1,1]*1.e-19  ;image range for ISRF images
win=0L

IF keyword_set(test) THEN BEGIN
  nwav_ref=1
ENDIF ELSE BEGIN
	nwav_ref=101
ENDELSE

;========= Restore informations at initial resolution
file=!phangs_data_dir+'/ISRF/WORK/'+source_name+'_ref_header.sav'
st_info=file_info(file)
IF st_info.exists NE 1 THEN BEGIN
	message,'Could not find '+file,/continue
	stop
ENDIF
restore,file,/verb
;% RESTORE: Restored variable: HREF.
href0=href & href=0

file=!phangs_data_dir+'/ISRF/WORK/'+source_name+'_isrf_min_prediction.sav'
st_info=file_info(file)
IF st_info.exists NE 1 THEN BEGIN
	message,'Could not find '+file,/continue
	stop
ENDIF
restore,file,/verb
;% RESTORE: Restored variable: ISRFS.
;% RESTORE: Restored variable: G0S.
;% RESTORE: Restored variable: OBJECT_DISTANCE.
;% RESTORE: Restored variable: OBJECT_THICKNESS.
;% RESTORE: Restored variable: USE_SOURCE_NAME.
ISRFS0=ISRFS & ISRFS=0
G0S0=G0S & G0s=0
;=== remove negative ISRFs values
ind=where(ISRFS0 LT 0 AND ISRFS0 NE la_undef(),count)
IF count NE 0 THEN ISRFS0[ind]=0.

file=!phangs_data_dir+'/ISRF/WORK/'+source_name+'_seds_indices.sav'
st_info=file_info(file)
IF st_info.exists NE 1 THEN BEGIN
	message,'Could not find '+file,/continue
	stop
ENDIF
restore,file,/verb
;% RESTORE: Restored variable: ALL_SEDS_INDICES.
ALL_SEDS_INDICES0=ALL_SEDS_INDICES & ALL_SEDS_INDICES=0

Nwav=(size(ISRFS0))[1]
Nvor=(size(ISRFS0))[2]

vors0=lindgen(Nvor)  ;for plots only

;=== This is to compute the Voronoi bins sizes (not used later)
vor_num=lonarr(Nvor)
FOR i=0L,Nvor-1 DO BEGIN
	IF ptr_valid(all_seds_indices0[i]) THEN BEGIN
		vor_num[i]=n_elements(*all_seds_indices0[i])
		IF vor_num[i] EQ 1 THEN BEGIN
			IF *all_seds_indices0[i] EQ -1 THEN BEGIN
				vor_num[i]=0L
			ENDIF
	  ENDIF
	ENDIF
ENDFOR
vor_sizes=2.*sqrt(1.*vor_num*sxpar(href0,'CDELT2')^2/!pi)   ;FWHM in deg
print,minmax(vor_sizes)*60.^2

;==== Make the ISRF cube
message,'Making the ISRF cube at initial resolution',/info
Nx=sxpar(href0,'NAXIS1')
Ny=sxpar(href0,'NAXIS2')
ISRF_cube0=fltarr(Nx,Ny,Nwav)+la_undef()

im0=ISRF_cube0[*,*,0]
;ISRF0=fltarr(Nwav)+la_undef()
FOR i=0L,Nvor-1 DO BEGIN
	IF i mod 1000 EQ 0 THEN BEGIN
		message,strtrim(1.*i/Nvor*100,2)+' %',/continue
		;stop
	ENDIF
	IF vor_num[i] NE 0 THEN BEGIN
		ij=index2ij(*all_seds_indices0[i],[Nx,Ny])
		FOR j=0L,vor_num[i]-1 DO ISRF_cube0[ij[j,0],ij[j,1],*]=ISRFS0[*,i]
	ENDIF
ENDFOR

IF keyword_set(show_plots) THEN BEGIN
	window,win & win=win+1
	tit='initial ISRF'
	image_cont20,ISRF_CUBE0[*,*,101],href0,/square,/silent,imrange=imrange,image_color_table='jpbloadct',axis_color_table=1,off_bar_pos=obp,tit=tit
ENDIF
IF not keyword_set(nostop) THEN stop

;==== smooth the ISRF
use_reso_filter='SPIRE3'
IF keyword_set(resolution_filter) THEN use_reso_filter=resolution_filter

message,'smoothing the ISRF cube',/info
data_reso=3./60./60.
reso_str='_'+use_reso_filter
final_reso=dustem_filter2reso(use_reso_filter)

;stop
;in fact degrade_res can degrade resolution of cubes, so ...
;=== That's a test
;goto,jump_there
IF keyword_set(test) THEN BEGIN
	ISRF_cube=degrade_res(ISRF_cube0[*,*,100:102],href0,data_reso,final_reso,hout)
	test_str='Test'
	Nwav=3
ENDIF ELSE BEGIN
	ISRF_cube=degrade_res(ISRF_cube0,href0,data_reso,final_reso,hout)
	test_str=''
ENDELSE
header=hout
goto,jump_there2

jump_there:
bidon=degrade_res(ISRF_cube0[*,*,0],href0,data_reso,final_reso,hout,/no_resample,FT_PSF=psf_FT)
;stop
ISRF_cube=fltarr(sxpar(hout,'NAXIS1'),sxpar(hout,'NAXIS2'),Nwav)
IF keyword_set(test) THEN BEGIN
	FOR k=100,102 DO BEGIN
		ISRF_cube[*,*,k-100]=degrade_res(ISRF_cube0[*,*,k],href0,data_reso,final_reso,hout,/no_resample,FT_PSF=psf_FT)
	ENDFOR
	test_str='Test'
	Nwav=3
ENDIF ELSE BEGIN
	FOR k=0L,Nwav-1 DO BEGIN
		message,'Degrading ISRF plane '+strtrim(k,2)+' / '+strtrim(Nwav,2),/info
		res=degrade_res(ISRF_cube0[*,*,k],href0,data_reso,final_reso,hout,/no_resample,FT_PSF=psf_FT)
		ISRF_cube[*,*,k]=res
	ENDFOR
	test_str=''
ENDELSE
header=hout

jump_there2:

IF keyword_set(show_plots) THEN BEGIN
	window,win & win=win+1
	tit='Smoothed ISRF'
	IF keyword_set(test) THEN BEGIN
		image_cont20,ISRF_CUBE[*,*,1],header,/square,/silent,imrange=imrange,image_color_table='jpbloadct',axis_color_table=1,off_bar_pos=obp,tit=tit
	ENDIF ELSE BEGIN
		image_cont20,ISRF_CUBE[*,*,101],header,/square,/silent,imrange=imrange,image_color_table='jpbloadct',axis_color_table=1,off_bar_pos=obp,tit=tit
	ENDELSE
ENDIF

IF not keyword_set(nostop) THEN stop

;=== project ISRFs on reference header for the given resolution
message,'project smoothed ISRF cube on reference header',/info
file=!phangs_data_dir+'/ISRF/WORK/'+source_name+'_ref_header'+reso_str+'.sav'
st_info=file_info(file)
IF st_info.exists NE 1 THEN BEGIN
	message,'Could not find '+file,/continue
	stop
ENDIF
restore,file,/verb
;% RESTORE: Restored variable: HREF.
;stop

NNx=sxpar(href,'NAXIS1')
NNy=sxpar(href,'NAXIS2')
proj_ISRF_cube=fltarr(NNx,NNy,Nwav)
i=0L
proj_ISRF_cube[*,*,i]=project2(header,ISRF_cube[*,*,i],href,/silent)
;FOR i=1L,Nwav-1 DO BEGIN
FOR i=1L,(size(proj_ISRF_CUBE))[3]-1 DO BEGIN
	proj_ISRF_cube[*,*,i]=project2(header,ISRF_CUBE[*,*,i],href,/previous_lut,/silent)
ENDFOR

IF keyword_set(show_plots) THEN BEGIN
	window,win & win=win+1
	tit='Projected ISRF'
	IF keyword_set(test) THEN BEGIN
		image_cont20,proj_ISRF_CUBE[*,*,1],href,/square,/silent,imrange=imrange,image_color_table='jpbloadct',axis_color_table=1,off_bar_pos=obp,tit=tit
	ENDIF ELSE BEGIN
		image_cont20,proj_ISRF_CUBE[*,*,101],href,/square,/silent,imrange=imrange,image_color_table='jpbloadct',axis_color_table=1,off_bar_pos=obp,tit=tit
	ENDELSE
ENDIF

IF not keyword_set(nostop) THEN stop

;==== put back ISRFs on voronoi bins
message,'putting back ISRFs on voronoi bins',/info
file=!phangs_data_dir+'/ISRF/WORK/'+use_source_name+'_muse_data'+reso_str+'.sav'
st_info=file_info(file)
IF st_info.exists NE 1 THEN BEGIN
	message,'Could not find '+file,/continue
	stop
ENDIF
restore,file,/verb
;% RESTORE: Restored variable: ST_TEMPLATES.
;% RESTORE: Restored variable: ST_MUSE_WEIGHTS.
;% RESTORE: Restored variable: VORONOI_ID.
;% RESTORE: Restored variable: AGE_VALUES.
;% RESTORE: Restored variable: METALICITY_VALUES.
;% RESTORE: Restored variable: BINS.
;% RESTORE: Restored variable: HREF.
file=!phangs_data_dir+'/ISRF/WORK/'+use_source_name+'_seds_indices'+reso_str+'.sav'
st_info=file_info(file)
IF st_info.exists NE 1 THEN BEGIN
	message,'Could not find '+file,/continue
	stop
ENDIF
restore,file,/verb
;% RESTORE: Restored variable: ALL_SEDS_INDICES.

IF keyword_set(show_plots) THEN BEGIN
	tit='Voronoi ID'
	window,win & win=win+1
	image_cont20,VORONOI_ID,href,/square,/silent,image_color_table='jpbloadct',axis_color_table=1,off_bar_pos=obp,tit=tit
ENDIF

NNvor=max(voronoi_id)+1
out_ISRFS=fltarr(Nwav,NNvor)
FOR i=0L,NNvor-1 DO BEGIN
	ind=where(voronoi_id EQ i,count)
	IF count NE 1 THEN BEGIN
		message,'count should be 1',/continue
		stop
	ENDIF ELSE BEGIN
		;ij=index2ij(ind,[NNx,NNy])
		ij2=index2ij(*all_seds_indices[i],[NNx,NNy])
    ;IF ij[0,0] NE ij2[0,0] OR ij[0,1] NE ij2[0,1] THEN stop
		out_ISRFS[*,i]=reform(proj_ISRF_cube[ij2[0,0],ij2[0,1],*])
		;ind=where(out_ISRFS[*,i] GT 1.e-14,count)
		;IF count NE 0 THEN stop
	ENDELSE
ENDFOR

vors=lindgen(NNvor)  ;for plots only

IF keyword_set(show_plots) THEN BEGIN
	window,win & win=win+1
	!p.multi=[0,1,2]
	indok=where(ISRFS0[100,*] NE la_undef())
	cgplot,vors0[indok],ISRFS0[101,indok],/ylog
	IF keyword_set(test) THEN BEGIN
		cgplot,vors,out_ISRFS[1,*],/ylog
	ENDIF ELSE BEGIN
		cgplot,vors,out_ISRFS[101,*],/ylog
	ENDELSE
ENDIF

;stop

IF not keyword_set(nostop) THEN stop

;print,la_min(ISRFS0),la_max(ISRFS0)
;  -2.0413261e-21   1.6803689e-16
;print,la_min(ISRF_cube),la_max(ISRF_cube)
; -1.48842e-10  1.52286e-10
;print,la_min(proj_ISRF_cube),la_max(proj_ISRF_cube)
; -1.48842e-10  1.52286e-10
;print,la_min(OUT_ISRFS),la_max(OUT_ISRFS)
; -1.48842e-10  1.52286e-10

;ISRFs=proj_ISRFs
ISRFs=out_ISRFs

IF keyword_set(save) THEN BEGIN
	file_save=!phangs_data_dir+'/ISRF/WORK/'+test_str+source_name+'_isrf_min_prediction'+reso_str+'.sav' 
	save,ISRFs,href,use_reso_filter,file=file_save
	message,'Saved '+file_save,/continue
ENDIF

the_end:
;stop

END