Blame view

LabTools/IRAP/JPB/phangs_make_astrosat_images.pro 3.98 KB
9beef1e9   Jean-Philippe Bernard   first commit
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
PRO phangs_make_astrosat_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_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
daa89ea5   Jean-Philippe Bernard   improved
64
restore,save_data_dir+use_source_name+'_jwst_images.sav',/verb
9beef1e9   Jean-Philippe Bernard   first commit
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
;% RESTORE: Restored variable: JWST_IMAGES.
;% RESTORE: Restored variable: FILTERS.
;% RESTORE: Restored variable: HREF.
;% RESTORE: Restored variable: NHCO.

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]   ;MJy/sr

FOR i=0L,Nastrosat_filter-1 DO BEGIN
	file=files[i]
	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)
	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 d=project2(h,d,href,/silent)
	tit=source_name+' '+astrosat_filters[i]
	image_cont20,d,href,/square,imrange=imrange,image_color_table='jpbloadct',/silent,tit=tit
	astrosat_images[*,*,0,i]=d
	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.sav'
	save,astrosat_images,astrosat_filters,href,file=file_save
	message,'Saved '+file_save,/info
ENDIF

the_end:

END