-
Notifications
You must be signed in to change notification settings - Fork 0
/
cubespectrum2.py
executable file
·187 lines (144 loc) · 4.2 KB
/
cubespectrum2.py
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
#! /usr/bin/env python
#
# Load a FITS cube , extract the spectrum at a (or reference) pixel
# and operate and plot some and then more....
#
#
# 22-jun-2017 PJT summer project - cloned off cubespectrum.py
# july-2017 Thomas/Peter various improvements
#
# @todo
# - have optional RESTFRQ or RESTFREQ as 3rd argument [done]
# - output the spectrum in a table, much like testCubeSpectrum.tab [done]
# - resample the gauss finer (not 5 points but may be 10x more?)
import os, sys, math
import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.units import Quantity
c = 299792.458 # [km/s] there should be a way to get 'c' from astropy.units ?
na = len(sys.argv)
if na == 7:
# Must be in Km/s
fitsfile = sys.argv[1]
pos = [int(sys.argv[2]),int(sys.argv[3])]
restfreq = float(sys.argv[4])* 1e9
vmin = float(sys.argv[5])
vmax = float(sys.argv[6])
use_vel = True
elif na == 5:
# Must be in GHz
fitsfile = sys.argv[1]
pos = [int(sys.argv[2]),int(sys.argv[3])]
vmin = vmax = None
restfreq = float(sys.argv[4])* 1e9
use_vel = True
elif na == 4:
# Pixel position
fitsfile = sys.argv[1]
pos = [int(sys.argv[2]),int(sys.argv[3])]
restfreq = None
vmin = vmax = None
use_vel = False
elif na == 2:
# Fits file
fitsfile = sys.argv[1]
pos = None
restfreq = None
vmin = vmax = None
use_vel = False
else:
sys.exit(1)
# open the fits file
hdu = fits.open(fitsfile)
print(len(hdu))
# get a reference to the header and data. Data should be 3dim numpy array now
h = hdu[0].header
d = hdu[0].data.squeeze()
print(d.shape)
# grab the restfreq, there are at least two ways how this is done
if restfreq == None:
if 'RESTFRQ' in h:
restfreq=h['RESTFRQ']
elif 'RESTFREQ' in h:
restfreq=h['RESTFREQ']
else:
restfreq= h['CRVAL3']
print("RESTFREQ",restfreq)
if pos == None:
# the FITS reference pixel is always a good backup
xpos = int(h['CRPIX1'])
ypos = int(h['CRPIX2'])
print("No position given, using reference pixel %g %g" % (xpos,ypos))
else:
xpos = pos[0]
ypos = pos[1]
flux = d[:,ypos,xpos]
nchan = d.shape[0]
channeln = np.arange(nchan)
zero = np.zeros(nchan)
cdelt3 = h['CDELT3']
crval3 = h['CRVAL3']
crpix3 = h['CRPIX3']
# to convert the channel to frequency
channelf = (channeln-crpix3+1)*cdelt3 + crval3
# to convert the Frequency to velocity
#channelv = (1.0-channelf/restfreq) * c
#print (channelf)
#print (channelv)
# what we plot
#channel = channelv
#channel = channelf
#channel = channeln
if use_vel:
# to convert the Frequency to velocity
channelv = (1.0-channelf/restfreq) * c
channel = channelv
print (channelv.min())
print (channelv.max())
else:
channel = channelf
print (channelf.min())
print (channelf.max())
ipeak = flux.argmax()
xpeak = channel[ipeak]
ypeak = flux[ipeak]
# moments around the peak
if na == 7:
m = 5
x = channel[ipeak-m:ipeak+m]
y = flux[ipeak-m:ipeak+m]
xmean = (x*y).sum() / y.sum()
xdisp = (x*x*y).sum() / y.sum() - xmean*xmean
if xdisp > 0:
xdisp = math.sqrt(xdisp)
fwhm = 2.355 * xdisp
print("MEAN/DISP/FWHM:",xmean,xdisp,fwhm)
ymodel = ypeak * np.exp(-0.5*(x-xmean)**2/(xdisp*xdisp))
if use_vel == True:
plt.figure()
if vmin != None:
channelv = ma.masked_outside(channelv,vmin,vmax)
plt.xlim([vmin,vmax])
plt.plot(channelv,flux,'o-',markersize=2,label='data')
plt.plot(channelv,zero)
# plt.plot(x,ymodel,label='gauss')
plt.xlabel("Velocity (km/s)")
plt.ylabel("Flux")
plt.title(fitsfile +" @ %g %g" % (xpos,ypos)+ " %g" % (restfreq/1e9)+ 'Ghz')
plt.legend()
plt.show()
else:
plt.figure()
plt.plot(channelf/1e9,flux,'o-',markersize=2,label='data')
plt.plot(channelf/1e9,zero)
plt.xlabel("Frequency (GHz)")
plt.ylabel("Flux")
plt.title(fitsfile + " @ %g %g" % (xpos,ypos))
plt.legend()
plt.show()
#to create a table of the frequency and flux
xtab = channelf /1e9 #to set the freqency to GHz
ytab = flux
np.savetxt('Frequency_Flux.tab',np.c_[xtab,ytab], delimiter=' ',header=("Frequency"" " "Flux"),comments='#',fmt='%.8f')