PRO phangs_make_astrosat_images,source_name=source_name,save=save,show_images=show_images,nostop=nostop,resolution_filter=resolution_filter,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')
;    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 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

reso_str=''
;stop
IF keyword_set(resolution_filter) THEN BEGIN
		reso_str='_'+resolution_filter
ENDIF

pdp_define_la_common
dustem_init

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

;data_dir=!phangs_data_dir+'/phangs_drive/PHANGS_ASTROSAT/v1p0/release/'
;data_dir=!phangs_data_dir+'/phangs_drive/Archive/PHANGS_ASTROSAT/v1p0/release/'
data_dir=!phangs_data_dir+'/phangs_drive/Archive/PHANGS_AstroSat/v1p0/release/'
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
file=save_data_dir+use_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 '+st_info,/continue
	stop
ENDIF
restore,file,/verb
;restore,save_data_dir+use_source_name+'_jwst_images.sav',/verb
;% RESTORE: Restored variable: HREF.

astrosat_filters=dustem_filter_names2filters(['F148W','F154W','F169M','F172M','N219M'])
Nastrosat_filter=n_elements(astrosat_filters)
astrosat_images=fltarr(sxpar(href,'NAXIS1'),sxpar(href,'NAXIS2'),2,Nastrosat_filter)
files=data_dir+strupcase(use_source_name)+'_'+['F148','F154','F169','F172','N219']+'_bkg_subtracted_mw_corrected.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(astrosat_filters)
cmic=3.e14   ;mic/sec
facts=wavs_mic/cmic*(wavs_mic*1.e4)/4./!pi*1.e-7*1.e4*1.e20    ;from Flambda in ergs/cm2/s/AA to MJy/sr (note: was divided by 4*pi)
imrange=[-0.5,2]*1e-17    ;erg/s/cm2/AA
imrange=[-0.5,1]/10   ;MJy/sr

win=0L

FOR i=0L,Nastrosat_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 '+file,/continue
		stop
	ENDIF
	d=mrdfits(file,0,h)
	steradians=(sxpar(h,'CDELT2')/!radeg)^2   ;pixel in sr
	fact=facts[i]*4.*!pi/steradians           ;This is to rescale from /pixel to /sr
	ind=where(finite(d) NE 1,count)
	IF count NE 0 THEN d[ind]=la_undef()
	d=la_mul(d,fact)
	window,win & win=win+1
	image_cont20,d,h,/square,imrange=imrange,image_color_table='jpbloadct',axis_color_table=1,/silent,tit=tit,off_bar_pos=obp
	reso_str=''
	;stop
	IF keyword_set(resolution_filter) THEN BEGIN
		reso_str='_'+resolution_filter
    data_reso=sxpar(h,'BMAJ')
		final_reso=dustem_filter2reso(resolution_filter)
  	d=degrade_res(d,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)
	ENDIF
	tit=source_name+' '+astrosat_filters[i]
	window,win & win=win+1
	image_cont20,d,href,/square,imrange=imrange,image_color_table='jpbloadct',axis_color_table=1,/silent,tit=tit,off_bar_pos=obp
	astrosat_images[*,*,0,i]=d
	IF not keyword_set(nostop) THEN stop
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(astrosat_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+'_astrosat_data'+reso_str+'.sav'     ;changed names on April 15 2024
	file_save=save_data_dir+use_source_name+'_astrosat_images'+reso_str+'.sav'
	save,astrosat_images,astrosat_filters,href,file=file_save
	message,'Saved '+file_save,/info
ENDIF

the_end:

END