phangs_make_herschel_images.pro
7.25 KB
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
64
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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
PRO phangs_make_herschel_images,source_name=source_name $
,save=save $
,show_images=show_images $
,nostop=nostop $
,help=help $
,resolution_filter=resolution_filter $
,resolution_value=resolution_value
;+
; NAME:
; phangs_make_herschel_images
; PURPOSE:
; makes astrosat images usable by the Phangs ISRF project
; CATEGORY:
; Dustem Phangs
; CALLING SEQUENCE:
; phangs_make_herschel_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 Herschel images are assumed to be in Jy/beam and Jy/pixel for SPIRE and PACS respectively
;
; EXAMPLES
; phangs_make_herschel_images,source_name='ngc0628',/show,/nostop
; phangs_make_herschel_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_herschel_images'
goto,the_end
ENDIF
use_source_name='ngc0628'
IF keyword_set(source_name) THEN use_source_name=source_name
dustem_define_la_common
;JPB: I'm runing dustem_init there to get the filter informations into memory
dustem_init
reso_str=''
IF keyword_set(resolution_filter) OR keyword_set(resolution_value) THEN BEGIN
IF keyword_set(resolution_filter) THEN BEGIN
reso_str='_'+resolution_filter
final_reso=dustem_filter2reso(resolution_filter)
ENDIF ELSE BEGIN
reso_str='_'+strtrim(resolution_value,2)
final_reso=resolution_value
ENDELSE
ENDIF
;window,0
obp=[1.1,0,1.15,1] ;JPB: see with image_cont20 how this is used. This is an offset position for where the color bar should be drawn.
win=0L
;JPB: modified this as this is the arborescence where the files are located in alma1.
; make sure you put them in the same arborescence in your PC
;data_dir=!phangs_data_dir+'/phangs_drive/PHANGS_Herschel_JWST_cycle1/'
data_dir=!phangs_data_dir+'/phangs_drive/Archive/Ancillary/herschel/PHANGS_Herschel_v1/'
save_data_dir=!phangs_data_dir+'/ISRF/WORK/'
;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,/continue
stop
ENDIF
restore,file,/verb
;% RESTORE: Restored variable: HREF.
;== I added PMW and PLW, since the data is available on alma1
;name_loc_filters = ['PACS_BLUE','PACS_GREEN','PACS_RED','PSW']
name_loc_filters = ['PACS_BLUE','PACS_GREEN','PACS_RED','PSW','PMW','PLW']
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)
;stop
;JPB: added missing SPIRE filters
name_filters = ['pacs70','pacs100','pacs160','spire250','spire350','spire500']
;=== JPB: added .gz, as this is how files are nammed in alma1 and mrdfits can actually read .gz files (!)
;files=data_dir+use_source_name+'_scanamorphos_v25_'+name_filters+'_0.fits.gz'
files=data_dir+use_source_name+'_scanamorphos_v25_'+name_filters+'_0.fits'
wavs_mic=dustem_filter2wav(herschel_filters)
cmic=3.e14 ;mic/sec
;JPB; as I want arrays which all have the same number of elements, I added dummy values (la_undef()) for PACS filters
;spire_beams_as = [450., 795., 1665.] ; in arcsec^2
herschel_beams_as = [la_undef(),la_undef(),la_undef(),450., 795., 1665.] ; in arcsec^2
;JPB: This below
;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
;can be done like this in one line:
;omegas=spire_beams_as*(((1.0/3600.0)*!pi/180.0)^2) ;omega_beam in sr
;JPB: This is the same as above but preserves undefined values
omegas=la_mul(herschel_beams_as,((1.0/3600.0)*!pi/180.0)^2) ;omega_beam in sr
;JPB: This array will tell you which instrument the filter is from (quite useful for your loop later)
instrus=dustem_filter2instru(herschel_filters)
imrange=[-10,100] ;MJy/sr
win=0L
FOR i=0L,Nherschel_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,/continue
stop
ENDIF
d_narray=mrdfits(file,0,h)
d = d_narray[*,*,0]
d_err = d_narray[*,*,1]
ind=where(finite(d) NE 1,count)
IF count NE 0 THEN d[ind]=la_undef()
;JPB: look, I've reduced 4 IF statements into 1 statement, and no need to use clumpsy STRPOS
IF instrus[i] EQ 'SPIRE' THEN BEGIN
fact=1.e-6/omegas[i]
ENDIF ELSE BEGIN
steradians=(sxpar(h,'CDELT2')/!radeg)^2 ;pixel in sr
fact=1.e-6/steradians ;This is to rescale from /pixel to /sr
ENDELSE
d=la_mul(d,fact)
IF keyword_set(show_images) THEN BEGIN
window,win & win=win+1
image_cont20,d,h,/square,imrange=0,image_color_table='jpbloadct',axis_color_table=1,/silent,tit=tit,off_bar_pos=obp
ENDIF
IF keyword_set(resolution_filter) or keyword_set(resolution_value) THEN BEGIN
data_reso=dustem_filter2reso(herschel_filters[i]) ;This is undefined in instrument_descritpion.xcat
IF final_reso GE data_reso THEN BEGIN
d=degrade_res(d,h,data_reso,final_reso,hout)
d_err=degrade_res(d_err,h,data_reso,final_reso,hout)
h=hout
ENDIF ELSE BEGIN
message,'Not convolving image as data_reso='+strtrim(data_reso,2)+' > finale_reso ='+strtrim(final_reso,2),/continue
stop
;left blank intentionnally
ENDELSE
ENDIF
print,la_min(d),la_max(d)
;stop
IF sxpar(h,'NAXIS1') NE sxpar(href,'NAXIS1') OR sxpar(h,'NAXIS2') NE sxpar(href,'NAXIS2') THEN BEGIN
d=project2(h,d,href,/silent,save_it='/tmp/saved_proj_info.sav')
d_err=project2(h,d_err,href,/silent,restore_it='/tmp/saved_proj_info.sav')
ENDIF
tit=use_source_name+' '+herschel_filters[i]
IF keyword_set(show_images) THEN BEGIN
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
IF not keyword_set(nostop) THEN stop
ENDIF
herschel_images[*,*,0,i]=d
herschel_images[*,*,1,i]=la_power(d_err,2.)
ENDFOR
;=== 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
;=== JPB: I changed the file name into images.sav, for compatibility with other files for other missions
;file_save=save_data_dir+use_source_name+'_herschel_data.sav'
file_save=save_data_dir+use_source_name+'_herschel_images'+reso_str+'.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:
END