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

;+
; 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')
; 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

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

;JPB : I modified this since we now have a save set that contains only the reference Fits header
;restore,save_data_dir+'ngc0628_jwst_images.sav',/verb

;% RESTORE: Restored variable: JWST_IMAGES.
;% RESTORE: Restored variable: FILTERS.
;% RESTORE: Restored variable: HREF.
;% RESTORE: Restored variable: NHCO.

restore,save_data_dir+'ngc0628_ref_header.sav',/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'

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]
	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

 	;if (STRPOS('spire250', name_filters[i]) NE -1) then begin
  ;   	fact = 1.e-6/omega_250
  ;endif
  ;if (STRPOS('spire350', name_filters[i]) NE -1) then begin
  ;     	fact = 1.e-6/omega_350
  ;endif
  ;if (STRPOS('spire500', name_filters[i]) NE -1) then begin
  ;     	fact = 1/e-6/omega_500
  ;endif
  ;if (STRPOS(name_filters[i],'spire') EQ -1) then begin
 	; 	steradians=(sxpar(h,'CDELT2')/!radeg)^2   ;pixel in sr
  ;  fact=1.e-6/steradians           ;This is to rescale from /pixel to /sr   	
  ;endif
              	

 	d=la_mul(d,fact)

	;=== JPB: I have re-instored image_cont20. It has to work on your machine, if it does not, we'll have to fix the problem
	;=== and you can always not run it by not setting the show_images keyword
	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

 	IF sxpar(h,'NAXIS1') NE sxpar(href,'NAXIS1') OR sxpar(h,'NAXIS2') NE sxpar(href,'NAXIS2') THEN BEGIN

      ;JPB: here something you could not know:
      ;project2 can store and restore the projection info. So, if you do the same projection twice, you can do this way:
      ;and gain a fcator of almost 2 in computing speed
     	;d=project2(h,d,href,/silent)
     	;d_err=project2(h,d_err,href,/silent)
     	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]

	;=== JPB: I have re-instored image_cont20. It has to work on your machine, if it does not, we'll have to fix the problem
	;=== and you can always not run it by not setting the show_images keyword
	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
	ENDIF

	herschel_images[*,*,0,i]=d

	;JPB: I believe that d_err is an uncertainty, not a variance, so it has to be squerd to go into this cube for the SED extractor
	;herschel_images[*,*,1,i]=d_err
	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.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