PRO phangs_make_hi_images,source_name=source_name,save=save,show_images=show_images,nostop=nostop,help=help,resolution_filter=resolution_filter

;+
; NAME:
;    phangs_make_hi_images
; PURPOSE:
;    makes HI images usable by the Phangs ISRF project
; CATEGORY:
;    Dustem Phangs
; CALLING SEQUENCE:
;    phangs_make_hi_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 HI images are asumed to be in units of K*km/s
; EXAMPLES
;	   phangs_make_hi_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_hi_images'
  goto,the_end
ENDIF

win=0L
dustem_define_la_common

dustem_init

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

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

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,/continue
	stop
ENDIF

restore,file,/verb
;% RESTORE: Restored variable: HREF.

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

;data_dir=!phangs_data_dir+'/phangs_drive/PHANGS_HI/'
data_dir=!phangs_data_dir+'/phangs_drive/Archive/PHANGS_HI/'

;file=data_dir+use_source_name+'_21cm_strictmask_mom0.fits'
file=data_dir+strupcase(use_source_name)+'_21cm_strictmask_mom0.fits'

st_info=file_info(file)
IF st_info.exists NE 1 THEN BEGIN
	message,'Could not find '+file,/continue
	stop
ENDIF

d=readfits(file,h)
;stop
;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()

reso_str=''
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 Nx OR sxpar(h,'NAXIS2') NE Ny THEN WHI=project2(h,d,href,/silent) ELSE WHI=d

fact=1.823e18/1.e21

NHI=la_mul(WHI,fact)   ;NH from HI in 1e21 H/cm2
tit=source_name+' '+'NHI [1e21 H/cm2]'
;stop
IF keyword_set(show_images) THEN BEGIN
	window,win & win=win+1
	image_cont20,NHI,href,/square,imrange=[-0.1,5],image_color_table='jpbloadct',/silent,tit=tit,axis_color_table=1
	IF not keyword_set(nostop) THEN stop
ENDIF

IF keyword_set(save) THEN BEGIN

	save_file=save_data_dir+use_source_name+'_HI_images'+reso_str+'.sav'

	save,NHI,href,file=save_file
	message,'Saved '+save_file,/continue
ENDIF

the_end:

END