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

;+
; NAME:
;    phangs_make_jwst_images
; PURPOSE:
;    makes JWST images usable by the Phangs ISRF project
; CATEGORY:
;    Dustem Phangs
; CALLING SEQUENCE:
;    phangs_make_jwst_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 JWST images are asumed to be in units of MJy/sr
; EXAMPLES
;	 phangs_make_jwst_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_jwst_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='/Volumes/PILOT_FLIGHT1/PHANGS-JWST/DR1/'
data_dir=!phangs_data_dir+'/phangs_drive/PHANGS-JWST/DR1/'
save_data_dir=!phangs_data_dir+'/ISRF/WORK/'

;========== get JWST reference header
file=data_dir+use_source_name+'_nircam_lv3_f300m_i2d_anchor_atgauss1.fits' 
d=mrdfits(file,1,h)
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
print,sxpar(h,'CDELT1')
href=cd2astro_header(h)
sxaddpar,href,'EQUINOX',2000.
IF keyword_set(show_images) THEN BEGIN
	window,win & win=win+1
	image_cont20,d,href,/square,imrange=[-0.1,5],image_color_table='jpbloadct',/silent,tit=tit
ENDIF

;=== save reference header
IF keyword_set(save) THEN BEGIN
	save_file=save_data_dir+use_source_name+'_ref_header.sav'
	save,href,file=save_file
	message,'Saved '+save_file,/continue
ENDIF

stop


Nx=sxpar(href,'NAXIS1')
Ny=sxpar(href,'NAXIS2')

filters_names=['F200W','F300M','F335M','F360M','F0770W','F1000W','F1130W','F2100W']
filters=dustem_filter_names2filters(filters_names)
Nfilters=n_elements(filters)
jwst_images=fltarr(Nx,Ny,2,Nfilters)

IF not keyword_set(nostop) THEN stop

;========== JWST stuff
i=0L
file=data_dir+use_source_name+'_nircam_lv3_f200w_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
help,d
;D               DOUBLE    = Array[7289, 9476]
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
IF keyword_set(show_images) THEN BEGIN
	window,win & win=win+1
	image_cont20,d,href,/square,imrange=[-0.1,5],image_color_table='jpbloadct',/silent,tit=tit
	IF not keyword_set(nostop) THEN stop
ENDIF

;stop

file=data_dir+use_source_name+'_nircam_lv3_f300m_i2d_anchor_atgauss1.fits' 
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
help,d
;D               DOUBLE    = Array[3515, 4576]
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
IF keyword_set(show_images) THEN BEGIN
	window,win & win=win+1
	image_cont20,d,href,/square,imrange=[-0.1,5],image_color_table='jpbloadct',/silent,tit=tit
	IF not keyword_set(nostop) THEN stop
ENDIF

file=data_dir+use_source_name+'_nircam_lv3_f335m_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
IF keyword_set(show_images) THEN BEGIN
	window,win & win=win+1
	image_cont20,d,href,/square,imrange=[-0.1,5],image_color_table='jpbloadct',/silent,tit=tit
	IF not keyword_set(nostop) THEN stop
ENDIF

file=data_dir+use_source_name+'_nircam_lv3_f360m_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
IF keyword_set(show_images) THEN BEGIN
	window,win & win=win+1
	image_cont20,d,href,/square,imrange=[-0.2,5],image_color_table='jpbloadct',/silent,tit=tit
	IF not keyword_set(nostop) THEN stop
ENDIF

file=data_dir+use_source_name+'_miri_lv3_f770w_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
IF keyword_set(show_images) THEN BEGIN
	window,win & win=win+1
	image_cont20,d,href,/square,imrange=[-0.2,10],image_color_table='jpbloadct',/silent,tit=tit,off_bar=obp
	IF not keyword_set(nostop) THEN stop
ENDIF

file=data_dir+use_source_name+'_miri_lv3_f1000w_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
IF keyword_set(show_images) THEN BEGIN
	window,win & win=win+1
	image_cont20,d,href,/square,imrange=[-0.2,5],image_color_table='jpbloadct',/silent,tit=tit
	IF not keyword_set(nostop) THEN stop
ENDIF

file=data_dir+use_source_name+'_miri_lv3_f1130w_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
IF keyword_set(show_images) THEN BEGIN
	window,win & win=win+1
	image_cont20,d,href,/square,imrange=[-0.3,20],image_color_table='jpbloadct',/silent,tit=tit
	IF not keyword_set(nostop) THEN stop
ENDIF
IF not keyword_set(nostop) THEN stop

file=data_dir+use_source_name+'_miri_lv3_f2100w_i2d_anchor_atgauss1.fits'
d=mrdfits(file,0,h0)
d=mrdfits(file,1,h)
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN d=project2(h,d,href,/silent)
jwst_images[*,*,0,i]=d & i=i+1
print,sxpar(h,'EXTNAME')
filter_name=sxpar(h0,'FILTER')
print,filter_name
tit=source_name+' '+filter_name
IF keyword_set(show_images) THEN BEGIN
	window,win & win=win+1
	image_cont20,d,href,/square,imrange=[-0.5,10],image_color_table='jpbloadct',/silent,tit=tit
	IF not keyword_set(nostop) THEN stop
ENDIF
IF not keyword_set(nostop) THEN stop

;Invent variances (will have to do better)
;jwst_images[*,*,1,*]=(jwst_images[*,*,0,*]*5./100.)^2   ;assumed intensity variance
perc_error=5./100.
jwst_images[*,*,1,*]=la_power(la_mul(jwst_images[*,*,0,*],perc_error),2)   ;assumed intensity variance
;=== check for null variances
ind=where(jwst_images[*,*,1,*] EQ 0.,count)
IF count NE 0 THEN BEGIN
	message,'There are null variances ...',/continue
	stop
ENDIF

;WCO map:
file=NH_data_dir+use_source_name+'_12m+7m+tp_co21_broad_mom0.fits'
;file='/Volumes/PILOT_FLIGHT1/PHANGS-JWST/'+use_source_name+'_12m+7m+tp_co21_broad_mom0.fits'
d=readfits(file,h)
;sxaddpar,h,'CTYPE1','RA---TAN'
;sxaddpar,h,'CTYPE2','RA---TAN'
sxaddpar,h,'EQUINOX',2000.
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
IF sxpar(h,'NAXIS1') NE Nx OR sxpar(h,'NAXIS2') NE Ny THEN WCO=project2(h,d,href,/silent) ELSE WCO=d
fact=4.e20/1.e21
NHCO=la_mul(WCO,fact)   ;NH from CO in 1e21 H/cm2
tit=source_name+' '+'NHCO [1e21 H/cm2]'
IF keyword_set(show_images) THEN BEGIN
	window,win & win=win+1
	image_cont20,NHCO,href,/square,imrange=[-0.5,10],image_color_table='jpbloadct',/silent,tit=tit
	IF not keyword_set(nostop) THEN stop
ENDIF

IF keyword_set(save) THEN BEGIN
	save_file=save_data_dir+use_source_name+'_jwst_images.sav'
	save,jwst_images,filters,href,NHCO,file=save_file
	message,'Saved '+save_file,/continue
ENDIF

the_end:

END