PRO phangs_make_herschel_images,source_name=source_name $
															 ,save=save $
															 ,show_images=show_images $
															 ,nostop=nostop $
															 ,help=help $
															 ,resolution_filter=resolution_filter

;+
; NAME:
;    phangs_make_herschel_images
; PURPOSE:
;    makes astrosat images usable by the Phangs ISRF project
; CATEGORY:
;    Dustem Phangs
; CALLING SEQUENCE:
;    phangs_make_herschel_images[,source_name=][,/save][,/show_images][,/nostop]
; INPUTS:
;    None
; OPTIONAL INPUT PARAMETERS:
;    source_name            : source name (default = 'ngc0628')
;    resolution_filter      : if set makes images at this resolution (default none)
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
;    save      = If set, save result
;    show_images= if set, show images
;	 nostop    = if set, does not stop
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    A file is written
; RESTRICTIONS:
;    None
; PROCEDURE:
;    input Herschel images are assumed to be in Jy/beam and Jy/pixel for SPIRE and PACS respectively
;    
; EXAMPLES
;	phangs_make_herschel_images,source_name='ngc0628',/show,/nostop
;	phangs_make_herschel_images,source_name='ngc0628',/save,/show,/nostop
; 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_make_herschel_images'
  goto,the_end
ENDIF

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

dustem_define_la_common
;JPB: I'm runing dustem_init there to get the filter informations into memory
dustem_init

;window,0
obp=[1.1,0,1.15,1]   ;JPB: see with image_cont20 how this is used. This is an offset position for where the color bar should be drawn.
win=0L

;JPB: modified this as this is the arborescence where the files are located in alma1.
;     make sure you put them in the same arborescence in your PC
;data_dir=!phangs_data_dir+'/phangs_drive/PHANGS_Herschel_JWST_cycle1/'
data_dir=!phangs_data_dir+'/phangs_drive/Archive/Ancillary/herschel/PHANGS_Herschel_v1/'

save_data_dir=!phangs_data_dir+'/ISRF/WORK/'

;This is to get the reference header
file=save_data_dir+'ngc0628_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.

;== I added PMW and PLW, since the data is available on alma1
;name_loc_filters = ['PACS_BLUE','PACS_GREEN','PACS_RED','PSW']
name_loc_filters = ['PACS_BLUE','PACS_GREEN','PACS_RED','PSW','PMW','PLW']
herschel_filters=dustem_filter_names2filters(name_loc_filters)
Nherschel_filter=n_elements(herschel_filters)
herschel_images=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),2,Nherschel_filter)

;JPB: added missing SPIRE filters
name_filters = ['pacs70','pacs100','pacs160','spire250','spire350','spire500']

;=== JPB: added .gz, as this is how files are nammed in alma1 and mrdfits can actually read .gz files (!)
;files=data_dir+use_source_name+'_scanamorphos_v25_'+name_filters+'_0.fits.gz'
files=data_dir+use_source_name+'_scanamorphos_v25_'+name_filters+'_0.fits'
wavs_mic=dustem_filter2wav(herschel_filters)
cmic=3.e14   ;mic/sec

;JPB; as I want arrays which all have the same number of elements, I added dummy values (la_undef()) for PACS filters
;spire_beams_as = [450., 795., 1665.] ; in arcsec^2
herschel_beams_as = [la_undef(),la_undef(),la_undef(),450., 795., 1665.] ; in arcsec^2
;JPB: This below
;omega_250 = spire_beams_as[0]*(((1.0/3600.0)*!pi/180.0)^2) ;omega_beam in sr
;omega_350 = spire_beams_as[1]*(((1.0/3600.0)*!pi/180.0)^2) ;omega_beam in sr
;omega_500 = spire_beams_as[2]*(((1.0/3600.0)*!pi/180.0)^2) ;omega_beam in sr
;can be done like this in one line:
;omegas=spire_beams_as*(((1.0/3600.0)*!pi/180.0)^2) ;omega_beam in sr
;JPB: This is the same as above but preserves undefined values
omegas=la_mul(herschel_beams_as,((1.0/3600.0)*!pi/180.0)^2) ;omega_beam in sr
;JPB: This array will tell you which instrument the filter is from (quite useful for your loop later)
instrus=dustem_filter2instru(herschel_filters)

imrange=[-10,100]   ;MJy/sr

FOR i=0L,Nherschel_filter-1 DO BEGIN
	file=files[i]
	st_info=file_info(file)
	IF st_info.exists NE 1 THEN BEGIN
		message,'Could not find '+file,/continue
		stop
	ENDIF
	d_narray=mrdfits(file,0,h)
	d = d_narray[*,*,0]
	d_err = d_narray[*,*,1]
 	ind=where(finite(d) NE 1,count)
	IF count NE 0 THEN d[ind]=la_undef()

  ;JPB: look, I've reduced 4 IF statements into 1 statement, and no need to use clumpsy STRPOS
  IF instrus[i] EQ 'SPIRE' THEN BEGIN
  	fact=1.e-6/omegas[i]
  ENDIF ELSE BEGIN
  	steradians=(sxpar(h,'CDELT2')/!radeg)^2   ;pixel in sr
    fact=1.e-6/steradians            ;This is to rescale from /pixel to /sr
  ENDELSE
              	
 	d=la_mul(d,fact)

	IF keyword_set(show_images) THEN BEGIN
		image_cont20,d,h,/square,imrange=0,image_color_table='jpbloadct',axis_color_table=1,/silent,tit=tit,off_bar_pos=obp
	ENDIF

	reso_str=''
	IF keyword_set(resolution_filter) THEN BEGIN
    ;stop
		reso_str='_'+resolution_filter
		data_reso=dustem_filter2reso(herschel_filters[i])      ;This is undefined in instrument_descritpion.xcat
		final_reso=dustem_filter2reso(resolution_filter)
  	d=degrade_res(d,h,data_reso,final_reso,hout)
  	d_err=degrade_res(d_err,h,data_reso,final_reso,hout)
    h=hout
	ENDIF

 	IF sxpar(h,'NAXIS1') NE sxpar(href,'NAXIS1') OR sxpar(h,'NAXIS2') NE sxpar(href,'NAXIS2') THEN BEGIN
     	d=project2(h,d,href,/silent,save_it='/tmp/saved_proj_info.sav')
     	d_err=project2(h,d_err,href,/silent,restore_it='/tmp/saved_proj_info.sav')
  ENDIF
	tit=use_source_name+' '+herschel_filters[i]
	IF keyword_set(show_images) THEN BEGIN
		image_cont20,d,href,/square,imrange=imrange,image_color_table='jpbloadct',axis_color_table=1,/silent,tit=tit,off_bar_pos=obp
		IF not keyword_set(nostop) THEN stop
	ENDIF
	herschel_images[*,*,0,i]=d
	herschel_images[*,*,1,i]=la_power(d_err,2.)
ENDFOR

;=== check for null variances
ind=where(herschel_images[*,*,1,*] EQ 0.,count)
IF count NE 0 THEN BEGIN
	message,'There are null variances ...',/continue
	stop
ENDIF

IF keyword_set(save) THEN BEGIN
	;=== JPB: I changed the file name into images.sav, for compatibility with other files for other missions
	;file_save=save_data_dir+use_source_name+'_herschel_data.sav'
	file_save=save_data_dir+use_source_name+'_herschel_images'+reso_str+'.sav'
	save,herschel_images,herschel_filters,href,file=file_save
	message,'Saved '+file_save,/info
ENDIF
; 
; mwrfits,herschel_images,save_data_dir+use_source_name+'_herschel_data.fits',href

the_end:

END