phangs_make_herschel_images.pro
7.3 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
PRO phangs_make_herschel_images,source_name=source_name,save=save,show_images=show_images,nostop=nostop,help=help
;+
; 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')
; 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
;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
;JPB : I modified this since we now have a save set that contains only the reference Fits header
;restore,save_data_dir+'ngc0628_jwst_images.sav',/verb
;% RESTORE: Restored variable: JWST_IMAGES.
;% RESTORE: Restored variable: FILTERS.
;% RESTORE: Restored variable: HREF.
;% RESTORE: Restored variable: NHCO.
restore,save_data_dir+'ngc0628_ref_header.sav',/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)
;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'
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
FOR i=0L,Nherschel_filter-1 DO BEGIN
file=files[i]
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
;if (STRPOS('spire250', name_filters[i]) NE -1) then begin
; fact = 1.e-6/omega_250
;endif
;if (STRPOS('spire350', name_filters[i]) NE -1) then begin
; fact = 1.e-6/omega_350
;endif
;if (STRPOS('spire500', name_filters[i]) NE -1) then begin
; fact = 1/e-6/omega_500
;endif
;if (STRPOS(name_filters[i],'spire') EQ -1) then begin
; steradians=(sxpar(h,'CDELT2')/!radeg)^2 ;pixel in sr
; fact=1.e-6/steradians ;This is to rescale from /pixel to /sr
;endif
d=la_mul(d,fact)
;=== JPB: I have re-instored image_cont20. It has to work on your machine, if it does not, we'll have to fix the problem
;=== and you can always not run it by not setting the show_images keyword
IF keyword_set(show_images) THEN BEGIN
image_cont20,d,h,/square,imrange=0,image_color_table='jpbloadct',axis_color_table=1,/silent,tit=tit,off_bar_pos=obp
ENDIF
IF sxpar(h,'NAXIS1') NE sxpar(href,'NAXIS1') OR sxpar(h,'NAXIS2') NE sxpar(href,'NAXIS2') THEN BEGIN
;JPB: here something you could not know:
;project2 can store and restore the projection info. So, if you do the same projection twice, you can do this way:
;and gain a fcator of almost 2 in computing speed
;d=project2(h,d,href,/silent)
;d_err=project2(h,d_err,href,/silent)
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]
;=== JPB: I have re-instored image_cont20. It has to work on your machine, if it does not, we'll have to fix the problem
;=== and you can always not run it by not setting the show_images keyword
IF keyword_set(show_images) THEN BEGIN
image_cont20,d,href,/square,imrange=imrange,image_color_table='jpbloadct',axis_color_table=1,/silent,tit=tit,off_bar_pos=obp
ENDIF
herschel_images[*,*,0,i]=d
;JPB: I believe that d_err is an uncertainty, not a variance, so it has to be squerd to go into this cube for the SED extractor
;herschel_images[*,*,1,i]=d_err
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.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