/
makemom.pro
484 lines (453 loc) · 20 KB
/
makemom.pro
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
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
PRO MAKEMOM, filename, errfile=errfile, rmsest=rmsest, maskfile=maskfile, $
vrange=vrange, xyrange=xyrange, dvref=dvref, baseroot=baseroot, $
thresh=thresh, chmin=chmin, edge=edge, guard=guard, smopar=smopar, $
senmsk=senmsk, replace0=replace0, mask0=mask0, kelvin=kelvin, dorms=dorms, $
pvmom0=pvmom0, gain2err=gain2err, useall=useall
;+
; NAME:
; MAKEMOM
;
; PURPOSE:
; Produce moment maps with masking from PPV image cubes. This program generates
; images of moments 0, 1, 2 and their respective errors, as well as peak brightness
; and peak S/N maps. User should provide an estimate of the noise in the cube.
;
; ARGUMENTS:
; FILENAME -- FITS data cube [no default]. 3rd axis should be velocity (M/S).
; ERRFILE -- FITS error map or cube. Note that any pixels blanked here will
; result in the corresponding data pixels being dropped.
; default: unset => error is assumed constant with position
; RMSEST -- estimate of channel noise, in same units as cube. If ERRFILE is
; given, the error cube is rescaled with this as minimum (e.g. center).
; default: unset => use error map if given, or FITS keyword RMS, or
; get rms from data cube if /DORMS is requested.
; MASKFILE -- external masking cube (value=1 for valid data, 0 otherwise).
; Must match the coordinate grid of the input data cube. This is
; applied AFTER any other requested masking operations.
; default: unset => mask is generated by the program
; VRANGE -- limit the velocity range for moment calculations [km/s, km/s]
; default: unset => whole velocity range is used
; XYRANGE -- limit the xy-region for moment calculations [x0, x1, y0, y1]
; default: unset => whole image region is used
; DVREF -- linewidth in km/s for calculating conservative emom0 map (emom0max).
; default: unset => use velocity range specified by VRANGE
; BASEROOT -- root name for output files
; default: based on FILENAME root
; THRESH -- starting intensity threshold for signal mask, in units of sigma
; default: unset => no masking applied
; CHMIN -- the THRESH mask must everywhere span at least CHMIN channels.
; default: 2
; EDGE -- final threshold for dilated mask, in units of sigma
; default: unset => mask is not dilated
; GUARD -- width of guard band in x, y, v dimensions [pix, pix, pix]
; The mask is extended by the specified no. of pixels in each direction
; if THRESH and/or EDGE are also specified. If only one value is given
; it is used for v; if two values are given they are used for x, y.
; default: unset => mask is not extended in any direction
; SMOPAR -- for mask generation, first degrade the cube angular resolution to
; <SMOPAR[0]> arcsec and/or smooth spectra with a Gaussian function of
; fwhm=<SMOPAR[1]> km/s for signal detection.
; default: [0.,0.] => no pre-smoothing applied
; SENMSK -- regions where the noise level is higher than <SENMSK>*min(ERRCUBE)
; (e.g. edges of field-of-view) will be masked out and excluded from
; conservative error estimates.
; default: unset => no additional masking applied
;
; OPTION KEYWORDS:
; REPLACE0 -- input cube values of 0.0 will be treated as missing data
; MASK0 -- mask pixels where some channels are blanked (e.g. FOV varies w/channel)
; KELVIN -- force conversion from Jy/beam to Kelvin
; DORMS -- make an estimate of the channel noise from the data. This is used to
; scale ERRFILE if given. Cannot be used if RMSEST is given.
; GAIN2ERR -- use this flag if the input cube is primary gain corrected and you
; provide the gain cube as ERRFILE. Noise cube is then proportional to
; 1/ERRFILE. Must be used with /DORMS option unless RMSEST is given
; or the 'RMS' keyword is provided in the FITS header.
; USEALL -- set this to use all channels for rms noise estimation (/DORMS). By
; default only the first 2 and last 2 channels are used.
; PVMOM0 -- produce additional mom0 images by collapsing cube along x & y axes
; (not constant RAs or Decs!)
;
; OUTPUTS:
; <BASEROOT>.*.fits
;
; HISTORY:
;
; 20150527 tw initial version based on Rui Xue's idl_moments code
; 20150601 tw implement guard parameter
; 20150610 tw give estr format I0
; 20150617 tw gain2err and useall parameters
; 20150620 tw use guard[2] in file name
; 20150712 tw output flux table from masked mom-0
; 20160131 tw don't output masked cube; allow 2D input mask; output 2D mask
; 20160201 tw output emom0max image
; 20160215 tw weight mean spectrum by variance
; 20160506 tw gzip compress the output mask cube; better support for fits.gz
; 20160609 tw additional code to avoid DATAMIN = NaN. Use FITS keyword RMS
; for GAIN2ERR if available. Convert to K after errcube finalized.
; 20160917 tw gzip compress the output error cube
; 20161021 tw rename convfac as osamp (it's an oversampling factor)
; 20161112 tw add chmin parameter
;
;-
; DEFAULT SETTINGS AND FILE NAMES
if n_elements(smopar) eq 0 then begin
smopar = [0.,0.]
smostr = ''
endif else if n_elements(smopar) eq 1 then begin
smopar = [smopar[0],0.]
smostr = 'sm'+string(smopar[0],format='(I0)')+'_'
endif else begin
smostr = 'sm'+string(smopar[0],format='(I0)')+'v'+$
string(smopar[1],format='(I0)')+'_'
endelse
if keyword_set(thresh) then begin
tstr = 't'+string(thresh,format='(F3.1)')
endif else begin
thresh = 0.
tstr = ''
endelse
if n_elements(chmin) eq 0 then begin
chmin = 2
endif
if keyword_set(edge) then begin
estr = 'e'+string(edge,format='(I0)')
endif else begin
edge = 0.
estr = ''
endelse
if keyword_set(guard) then begin
if n_elements(guard) eq 1 then guard = [0,0,guard[0]]
if n_elements(guard) eq 2 then guard = [guard[0],guard[1],0]
gstr = 'g'+string(guard[2],format='(I0)')
endif else begin
guard = [0, 0, 0]
gstr = ''
endelse
if strmatch(filename,'*.fits.gz',/f) then galname=file_basename(filename,'.fits.gz',/f)
if strmatch(filename,'*.fits',/f) then galname=file_basename(filename,'.fits',/f)
if ~(strmatch(filename,'*.fits',/f) or strmatch(filename,'*.fits.gz',/f)) then begin
galname=file_basename(filename)
message,'Input file extension not recognized',/info
endif
if keyword_set(baseroot) eq 0 then begin
baseroot = galname + '.' + smostr + tstr + estr + gstr
if strpos(baseroot,'.',/reverse_search) eq strlen(baseroot)-1 then $
baseroot=strmid(baseroot,0,strlen(baseroot)-1)
endif
file_mkdir,file_dirname(baseroot,/m)
; READ IN DATA
data = READFITS(filename, hd, /silent)
RADIOHEAD, hd, s = h
; EXTRACT SUBREGION IF REQUESTED
if n_elements(xyrange) eq 4 then begin
HEXTRACT3D, data, hd, tmp, tmphd, xyrange
SXADDPAR, tmphd, 'DATAMAX', max(tmp,/nan), before='HISTORY'
SXADDPAR, tmphd, 'DATAMIN', min(tmp,/nan), before='HISTORY'
;WRITEFITS, baseroot+'.subreg.fits', float(tmp), tmphd
data=tmp
hd=tmphd
endif
; OUTPUT SOME INFO
print,replicate('-',35)
print,'spectral cube size: ',size(data,/d)
print,'pixel spacing : '+string(abs(h.cdelt[0])*3600.,format='(F5.2)') $
+' arcsec x '+string(abs(h.cdelt[1])*3600.,format='(F5.2)') $
+' arcsec x '+string(h.cdelt[2]/1.e3,format='(F5.1)')+' km/s'
print,'beam FWHM maj, min: '+string(h.bmaj,format='(F6.2)')+' ' $
+string(h.bmin,format='(F6.2)')+' eff: ' $
+string(sqrt(h.bmaj*h.bmin),format='(F6.2)')+' arcsec'
print,replicate('-',35)
; BLANKING
if keyword_set(replace0) then begin
data[where(data eq 0.0,/null)]=!values.f_nan
endif
; READ IN INPUT MASK AND SET VELOCITY WINDOW IF GIVEN
sz = size(data)
if n_elements(maskfile) eq 0 then begin
exmask=data*0.0+1.0
endif else begin
exmask=READFITS(maskfile, exmaskhd, /silent)
if (size(exmask))[0] eq 2 or (size(exmask))[3] eq 1 then begin
exmask0=exmask
exmask=make_array(sz[1],sz[2],sz[3],/float,/nozero)
for i=0,sz[3]-1 do exmask[0,0,i]=exmask0
endif
exmask[where(finite(exmask,/nan),/null)]=0.0
if n_elements(xyrange) eq 4 then begin
HEXTRACT3D, exmask, exmaskhd, tmp, tmphd, xyrange
SXADDPAR, tmphd, 'DATAMAX', max(tmp,/nan), before='HISTORY'
SXADDPAR, tmphd, 'DATAMIN', min(tmp,/nan), before='HISTORY'
exmask=tmp
exmaskhd=tmphd
endif
endelse
if n_elements(vrange) ne 0 then begin
tag_outvrange = where(h.v lt vrange[0] or h.v gt vrange[1])
if tag_outvrange[0] ne -1 then exmask[*,*,[tag_outvrange]]=0.0
endif else begin
vrange = [min(h.v),max(h.v)]
tag_outvrange = -1
endelse
in_vrange = where(h.v ge vrange[0] and h.v le vrange[1])
; GENERATE ERROR CUBE AND 2D MAP IF NOT GIVEN
if n_elements(errfile) eq 0 then begin
if keyword_set(dorms) then begin
ecube=ERR_CUBE(data,useall=useall)
endif else begin
if keyword_set(rmsest) then rms=rmsest else begin
rms = SXPAR( hd, 'RMS', count=ct )
print,'Using error from FITS header: ', rms
if (ct eq 0) then rms = 1.
endelse
ecube=data*0.0 + rms
endelse
; OR ELSE USE PROVIDED ERROR CUBE
endif else begin
ecube = readfits(errfile, ehd, /silent)
if n_elements(xyrange) eq 4 then begin
HEXTRACT3D,ecube,ehd,tmp,tmphd,xyrange
SXADDPAR, tmphd, 'DATAMAX', max(tmp,/nan), before='HISTORY'
SXADDPAR, tmphd, 'DATAMIN', min(tmp,/nan), before='HISTORY'
ecube=tmp
ehd=tmphd
endif
esz = size(ecube)
if esz[0] eq 2 then begin
ecube0=ecube
ecube=make_array(esz[1],esz[2],sz[3],/float,/nozero)
for i=0,sz[3]-1 do ecube[0,0,i]=ecube0
endif
ecube[where(ecube eq 0,/null)]=!values.f_nan
if keyword_set(dorms) then begin
if keyword_set(gain2err) then ecube = 1./ecube
tmp=ERR_CUBE(data,pattern=ecube,useall=useall)
ecube=tmp
endif else begin
if keyword_set(gain2err) then begin
ecube = 1./ecube
if keyword_set(rmsest) then begin
ecube=rmsest*ecube/min(ecube,/nan)
endif else begin
rms = SXPAR( hd, 'RMS', count=ct )
print,'Scaling by error from FITS header: ', rms
if (ct gt 0) then ecube=rms*ecube/min(ecube,/nan)
endelse
endif else begin
if keyword_set(rmsest) then begin
ecube=rmsest*ecube/min(ecube,/nan)
endif
endelse
endelse
data[where(ecube ne ecube,/null)]=!values.f_nan ; ignore data with no errors
endelse
; CONVERT JY/BEAM TO KELVIN, GENERATE 2D ERROR MAP
if keyword_set(kelvin) and strpos(strupcase(sxpar(hd,'BUNIT')),'JY/B') ne -1 then begin
data = temporary(data) * h.jypb2k
SXADDPAR, hd, 'BUNIT', 'K'
ecube=temporary(ecube) * h.jypb2k
SXADDPAR, ehd, 'BUNIT', 'K'
endif
emap = total(ecube, 3, /nan) / (total(ecube eq ecube, 3)>1)
emap[where(emap eq 0.0,/null)] = !values.f_nan
; GENERATE MASKING CUBE (missing data locations are still kept in mask)
mask = GENMASK(data,err=ecube,hd=hd,spar=smopar,sig=thresh,chmin=chmin,$
grow=edge,guard=guard)
mask = mask*exmask
; MASK HIGH NOISE AT EDGES
if n_elements(senmsk) ne 0 then begin
mask1d=float((total(ecube, 3) le senmsk*min(total(ecube, 3),/nan)))
mask1d[where(mask1d eq 0.0,/null)]=!values.f_nan
for i = 0, sz[3]-1 do begin
data[*,*,i] = data[*,*,i] + mask1d - mask1d
ecube[*,*,i] = ecube[*,*,i] + mask1d - mask1d
mask[*,*,i] = mask[*,*,i] + mask1d - mask1d
exmask[*,*,i] = exmask[*,*,i] + mask1d - mask1d
endfor
endif
; CALCULATE MOMENTS
; h.v: velocity in km/s
; note that the vrange input is obeyed for snrpk image
if n_elements(dvref) ne 0 then begin
nchref=ceil(dvref/(abs(h.cdelt[2])/1.0e3))*1.0
endif else begin
nchref=0.
endelse
MASKMOMENT, data[*,*,[in_vrange]], mask[*,*,[in_vrange]], $
ecube[*,*,[in_vrange]], h.v[in_vrange], $
mom0 = mom0, mom1 = mom1, mom2 = mom2, $
emom0 = emom0, emom1 = emom1, emom2 = emom2, $
peak = peak, snrpk=snrpk,$
mask0=mask0, nchref=nchref, emommx = emommx
mhd = hd
bunit = SXPAR(mhd,'BUNIT')
; UPDATE HISTORY IN FITS HEADER
histlabel = 'IDL_MOMMAPS: '
SXADDPAR, mhd, 'HISTORY', histlabel+systime()
SXADDPAR, mhd, 'HISTORY', histlabel+'filename='+filename
if n_elements(errfile) gt 0 then $
SXADDPAR, mhd, 'HISTORY', histlabel+'errfile='+errfile
if n_elements(maskfile) gt 0 then $
SXADDPAR, mhd, 'HISTORY', histlabel+'maskfile='+maskfile
SXADDPAR, mhd, 'HISTORY', histlabel+'smopar=['+strcompress(smopar[0],/r)+','+$
strcompress(smopar[1],/r)+']'
SXADDPAR, mhd, 'HISTORY', histlabel+'thresh='+strcompress(thresh,/r)
SXADDPAR, mhd, 'HISTORY', histlabel+'chmin='+strcompress(chmin,/r)
SXADDPAR, mhd, 'HISTORY', histlabel+'edge='+strcompress(edge,/r)
SXADDPAR, mhd, 'HISTORY', histlabel+'guard=['+strcompress(guard[0],/r)+$
','+strcompress(guard[1],/r)+','+strcompress(guard[2],/r)+']'
if n_elements(dvref) gt 0 then $
SXADDPAR, mhd, 'HISTORY', histlabel+'dvref='+strcompress(dvref,/r)
if n_elements(vrange) gt 0 then $
SXADDPAR, mhd, 'HISTORY', histlabel+'vrange=['+strcompress(vrange[0],/r)+$
','+strcompress(vrange[1],/r)+']'
; OUTPUT MASK CUBE
SXADDPAR,mhd,'DATAMAX',2.0, before='HISTORY'
SXADDPAR,mhd,'DATAMIN',-1.0, before='HISTORY'
nan_tag=where(data ne data,nan_ct)
if nan_ct ne 0 then mask[nan_tag]=!values.f_nan
if total(finite(mask)) ge 1 then $
WRITEFITS,baseroot+'.mask.fits',float(mask),mhd,/compress
; OUTPUT MASKED FLUX SPECTRUM WITH FORMAL AND CONSERVATIVE ERRORS
if strpos(strupcase(sxpar(mhd,'BUNIT')),'JY/B') ne -1 then begin
osamp = h.ppbeam
fluxunit = 'JY'
intfluxunit = 'JY*KM/S'
endif else begin
osamp = 1.
fluxunit = strtrim(bunit,2)+'*PIX'
intfluxunit = strtrim(bunit,2)+'KM/S*PIX'
endelse
; WEIGHTED MEAN SPECTRUM WITH FORMAL ERRORS
wgtdata = total(total(mask*data/ecube^2.,1,/nan),1,/nan)
weights1 = total(total(mask/ecube^2.,1,/nan),1,/nan)
weights1[where(weights1 eq 0.0,/null)] = 1.
meandata = (wgtdata/weights1)
meanerr1 = sqrt(1./weights1)
; CONSERVATIVE ERROR EST.
weights2 = total(total(exmask/ecube^2.,1,/nan),1,/nan)
weights2[where(weights2 eq 0.0,/null)] = 1.
meanerr2 = sqrt(1./weights2)
; INTEGRATED SPECTRUM AND FLUX
fluxdata = total(total( mask*data ,1,/nan),1,/nan) / osamp
fluxerr1 = total(total(( mask*ecube)^2.,1,/nan),1,/nan) * h.ppbeam
fluxerr2 = total(total((exmask*ecube)^2.,1,/nan),1,/nan) * h.ppbeam
intfluxdata = total(fluxdata,/nan) * abs(h.cdelt[2])/1.0e3
intfluxerr1 = sqrt(total(fluxerr1,/nan))/osamp*abs(h.cdelt[2])/1.0e3
intfluxerr2 = sqrt(total(fluxerr2,/nan))/osamp*abs(h.cdelt[2])/1.0e3
fluxerr1 = sqrt(fluxerr1) / osamp
fluxerr2 = sqrt(fluxerr2) / osamp
; OUTPUT WEIGHTED MEAN SPECTRUM
write_csv,baseroot+'.mean.out',h.v,meandata,meanerr1,meanerr2, $
header=['# '+h.ctype[2],'Brightness('+sxpar(mhd,'BUNIT')+')',$
'FormErr','ConsErr']
; OUTPUT FLUX VECTOR
fluxinfo = string(intfluxdata)+' '+intfluxunit+' +/- '+ $
string(intfluxerr1,format='(F0.2)')+' (formal) +/- '+ $
string(intfluxerr2,format='(F0.2)')+' (cons)'
print,'moment 0 flux: ',fluxinfo
write_csv,baseroot+'.flux.out',h.v,fluxdata,fluxerr1,fluxerr2, $
header=['# '+h.ctype[2],'Flux('+fluxunit+')','FormErr','ConsErr'], $
table_header='# '+fluxinfo
; OUTPUT 2D ERROR MAP IF NO ERROR MAP PROVIDED
if n_elements(errfile) eq 0 then begin
SXADDPAR, mhd, 'DATAMAX', max(emap,/nan), before='HISTORY'
SXADDPAR, mhd, 'DATAMIN', min(emap,/nan), before='HISTORY'
WRITEFITS, baseroot+'.rms.fits', float(emap), mhd
; OR OUTPUT RESCALED ERROR CUBE IF /DORMS or /RMSEST WAS SET
endif else if keyword_set(dorms) or keyword_set(rmsest) then begin
if esz[0] eq 2 then begin
emap = ecube[*,*,0]
SXADDPAR, mhd, 'DATAMAX', max(emap,/nan), before='HISTORY'
SXADDPAR, mhd, 'DATAMIN', min(emap,/nan), before='HISTORY'
WRITEFITS, baseroot+'.rms.fits', float(emap), mhd
endif else begin
SXADDPAR, mhd, 'DATAMAX', max(ecube,/nan), before='HISTORY'
SXADDPAR, mhd, 'DATAMIN', min(ecube,/nan), before='HISTORY'
WRITEFITS, baseroot+'.ecube.fits', float(ecube), mhd,/compress
endelse
endif
; DROP 3RD AND 4TH AXIS KEYWORDS FOR REST OF OUTPUTS
SXDELPAR, mhd, ['CTYPE3','CRVAL3','CRPIX3','CDELT3','CUNIT3']
SXDELPAR, mhd, ['CTYPE4','CRVAL4','CRPIX4','CDELT4','CUNIT4']
; PEAK INTENSITY (BUNIT IS SAME)
if total(finite(peak)) ge 1 then begin
SXADDPAR, mhd, 'DATAMAX', max(peak,/nan), before='HISTORY'
SXADDPAR, mhd, 'DATAMIN', min(peak,/nan), before='HISTORY'
WRITEFITS, baseroot+'.peak.fits', float(peak), mhd
endif else begin
print,'Peak intensity is completely blank: not written'
; SXADDPAR, mhd, 'DATAMAX', 0.0, before='HISTORY'
; SXADDPAR, mhd, 'DATAMIN', 0.0, before='HISTORY'
endelse
; PEAK SNR (UNMASKED)
SXADDPAR, mhd, 'BUNIT', 'SNR', before='HISTORY'
SXADDPAR, mhd, 'DATAMAX', max(snrpk,/nan), before='HISTORY'
SXADDPAR, mhd, 'DATAMIN', min(snrpk,/nan), before='HISTORY'
WRITEFITS, baseroot+'.snrpk.fits', float(snrpk), mhd
; MOMENT 0 AND ERROR
SXADDPAR, mhd, 'BUNIT', strtrim(bunit,2)+'.KM/S', before='HISTORY'
mom0gm = mom0 * abs(h.cdelt[2]) / 1.0e3
if total(finite(mom0gm)) ge 1 then begin
SXADDPAR, mhd, 'DATAMAX', max(mom0gm,/nan), before='HISTORY'
SXADDPAR, mhd, 'DATAMIN', min(mom0gm,/nan), before='HISTORY'
WRITEFITS, baseroot+'.mom0.fits', float(mom0gm), mhd
endif else begin
print,'Moment 0 is completely blank: not written'
endelse
emom0gm = emom0 * abs(h.cdelt[2]) / 1.0e3
if total(finite(emom0gm)) ge 1 then begin
SXADDPAR, mhd, 'DATAMAX', max(emom0gm,/nan), before='HISTORY'
SXADDPAR, mhd, 'DATAMIN', min(emom0gm,/nan), before='HISTORY'
WRITEFITS, baseroot+'.emom0.fits', float(emom0gm), mhd
endif
; CONSERVATIVE MOM0 ERROR
emommx = emommx * abs(h.cdelt[2]) / 1.0e3
SXADDPAR, mhd, 'DATAMAX', max(emommx,/nan), before='HISTORY'
SXADDPAR, mhd, 'DATAMIN', min(emommx,/nan), before='HISTORY'
WRITEFITS, baseroot+'.emom0max.fits', float(emommx), mhd
; MOMENT 1 AND ERROR
SXADDPAR, mhd, 'BUNIT', 'KM/S', before='HISTORY'
if total(finite(mom1)) ge 1 then begin
SXADDPAR, mhd, 'DATAMAX', max(mom1,/nan), before='HISTORY'
SXADDPAR, mhd, 'DATAMIN', min(mom1,/nan), before='HISTORY'
WRITEFITS, baseroot+'.mom1.fits', float(mom1), mhd
endif else begin
print,'Moment 1 is completely blank: not written'
endelse
if total(finite(emom1)) ge 1 then begin
SXADDPAR, mhd, 'DATAMAX', max(emom1,/nan), before='HISTORY'
SXADDPAR, mhd, 'DATAMIN', min(emom1,/nan), before='HISTORY'
WRITEFITS, baseroot+'.emom1.fits', float(emom1), mhd
endif
; PLOT MOMENT 0 AND MOMENT 1
PLTMOM, baseroot
; MOMENT 2 AND ERROR
if total(finite(mom2)) ge 1 then begin
SXADDPAR, mhd, 'DATAMAX', max(mom2,/nan), before='HISTORY'
SXADDPAR, mhd, 'DATAMIN', min(mom2,/nan), before='HISTORY'
WRITEFITS, baseroot+'.mom2.fits', float(mom2), mhd
endif else begin
print,'Moment 2 is completely blank: not written'
endelse
if total(finite(emom2)) ge 1 then begin
SXADDPAR, mhd, 'DATAMAX', max(emom2,/nan), before='HISTORY'
SXADDPAR, mhd, 'DATAMIN', min(emom2,/nan), before='HISTORY'
WRITEFITS, baseroot+'.emom2.fits', float(emom2), mhd
endif
; 2D PROJECTION OF THE MASK
if thresh gt 0.0 and total(mask,/nan) gt h.ppbeam then begin
SXADDPAR, mhd, 'BUNIT', ' ', before='HISTORY'
SXADDPAR, mhd, 'DATAMAX', 2.0, before='HISTORY'
SXADDPAR, mhd, 'DATAMIN',-1.0, before='HISTORY'
WRITEFITS,baseroot+'.msk2d.fits',max(float(mask),/nan,dim=3),mhd
endif
; MOM0 in XV and YV (OPTIONAL)
if keyword_set(pvmom0) and total(finite(mask)) ge 1 then begin
MASKMOMENT_PV,data,hd,mask,mom0xv,mom0vy,mom0xvhd=mom0xvhd,mom0vyhd=mom0vyhd,$
vrange=vrange
WRITEFITS, baseroot+'.mom0xv.fits',mom0xv,mom0xvhd
WRITEFITS, baseroot+'.mom0vy.fits',mom0vy,mom0vyhd
PLTMOM_PV, baseroot
endif
END