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

;+
; NAME:
;    phangs_make_astrosat_images
; PURPOSE:
;    makes astrosat images usable by the Phangs ISRF project
; CATEGORY:
;    Dustem Phangs
; CALLING SEQUENCE:
;    phangs_make_astrosat_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 astrosat images are asumed to be in Angstrom-1 cm-2 erg s-1 PER PIXEL (even though BUNIT   = 'Angstrom-1 cm-2 erg s-1' in header)
;    
; EXAMPLES
;	phangs_make_astrosat_images,source_name='ngc0628',/show,/nostop
;	phangs_make_astrosat_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_astrosat_images'
  goto,the_end
ENDIF

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

;pdp_define_la_common
;window,0
obp=[1.1,0,1.15,1]
win=0L

data_dir=!phangs_data_dir+'/phangs_drive/PHANGS_Herschel_JWST_cycle1/'
save_data_dir=!phangs_data_dir+'/ISRF/WORK/'

astrosat_data:
;BUNIT   = 'Angstrom-1 cm-2 erg s-1'
;in fact Angstrom-1 cm-2 erg s-1 PER PIXEL ?!

;This is to get the reference 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.

name_loc_filters = ['PACS_BLUE','PACS_GREEN','PACS_RED','PSW']
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)

name_filters = ['pacs70','pacs100','pacs160','spire250']

files=data_dir+use_source_name+'_scanamorphos_v25_'+name_filters+'_0.fits'
;NGC0628_F148_bkg_subtracted_mw_corrected.fits	NGC0628_F169_bkg_subtracted_mw_corrected.fits	NGC0628_N219_bkg_subtracted_mw_corrected.fits
;NGC0628_F154_bkg_subtracted_mw_corrected.fits	NGC0628_F172_bkg_subtracted_mw_corrected.fits
wavs_mic=dustem_filter2wav(herschel_filters)
cmic=3.e14   ;mic/sec

spire_beams_as = [450., 795., 1665.] ; in arcsec^2


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

imrange=[-0.5,2]*1e-17    ;erg/s/cm2/AA
imrange=[-0.5,1]   ;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 (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
              	
	IF count NE 0 THEN d[ind]=la_undef()
 	d=la_mul(d,fact)
 	
	;image_cont20,d,h,/square,imrange=imrange,image_color_table='jpbloadct',/silent,tit=tit
 	IF sxpar(h,'NAXIS1') NE sxpar(href,'NAXIS1') OR sxpar(h,'NAXIS2') NE sxpar(href,'NAXIS2') THEN BEGIN
     	d=project2(h,d,href,/silent)
     	d_err=project2(h,d_err,href,/silent)
    ENDIF
	tit=use_source_name+' '+herschel_filters[i]
	;image_cont20,d,href,/square,imrange=imrange,image_color_table='jpbloadct',/silent,tit=tit
	herschel_images[*,*,0,i]=d
	herschel_images[*,*,1,i]=d_err
	
ENDFOR

;Invent variances (will have to do better)
; perc_error=5./100.
; astrosat_images[*,*,1,*]=la_power(la_mul(astrosat_images[*,*,0,*],perc_error),2)   ;assumed intensity variance
;=== 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
	file_save=save_data_dir+use_source_name+'_herschel_data.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:
print,'bonjour'

END