-
Notifications
You must be signed in to change notification settings - Fork 0
/
visgen_streamers.py
292 lines (253 loc) · 14.7 KB
/
visgen_streamers.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
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
import math
import numpy
import pyfits
import os
from disk import Disk
# TEMPO: NONREDUCED CHI
"""
For the generator to work, have folders named "fits" and "images" in
the directory from which you are calling this script. There should
also be a folder named "hd61005" in which is stored the visibility
file "hd61005_all.vis". From there, just run with command line arguments:
number of pixels across, inclination angle, rotation angle, and name of
the fits file you would like created in the fits directory. Everything
in the image directory is temporary and will be deleted on the next run.
The "images" directory holds the generated visibilities for the model
and the data, both with file extension "uvf".
"""
class VisibilityGenerator:
"""
creates a visibility generator with width pixels,
angle inclination and rotation rotation (in degrees)
"""
def __init__(self, disk, width, angle, rotation, disk_scale, streamer_scale, fits_file):
self.width = width
# convert to radians
self.incline = angle*math.pi/180.0
self.rotation = (rotation+90)*math.pi/180.0
# Flux scaling factor
self.disk_scale = disk_scale
self.streamer_scale = streamer_scale
# name files we need
self.fits_name = fits_file
self.fits_model_image = 'images/' + fits_file
self.mir_model_image = 'images/model.mp'
self.mir_model_vis = 'images/model.vis'
self.fits_model_vis = 'images/model.uvf'
# get the data files we need
self.mir_data_vis = 'hd61005/hd61005_all.ppm.vis'
self.fits_data_vis = 'hd61005/hd61005_all.ppm.uvf'
# get the data from the observed visibility
if not os.path.exists(self.fits_data_vis):
command('fits in=' + self.mir_data_vis + ' op=uvout out=' + self.fits_data_vis)
# store data for computing the chi squared later on
observed = pyfits.getdata(self.fits_data_vis, 'PRIMARY').data[:,0,0,0,0]
self.data_vis_real = observed[:,0]
self.data_vis_imag = observed[:,1]
self.data_vis_weight = observed[:,2]
#Dimensional Information
self.outerRadius = disk.getOuterRadius()
self.innerRadius = disk.getInnerRadius()
self.majoraxis = self.outerRadius*1.496e11
self.minoraxis = self.majoraxis*math.sin(math.pi/2 - self.incline)
# square images
self.center = self.width/2.0
# now find the scaling factor
self.pxwid = 5.0*self.outerRadius/self.width
#Streamer Information
self.streamer_angle = 23.0*math.pi/180.0 #From Buenzli et al, 2010
self.streamer_powerlaw = -4.4 #From Buenzli et al, 2010
self.edge = (self.innerRadius + self.outerRadius)/2 #In AU
self.F_edge = disk.calculatePointFlux(self.edge, 225.5)/math.cos(self.incline)
self.dS_edge = math.sqrt(1+((self.edge*1.496e11)**2+(self.outerRadius*1.496e11)**2/math.tan(self.streamer_angle)**2)/((self.outerRadius*1.496e11)**2-(self.edge*1.496e11)**2))
self.width_edge = 2*math.sqrt((self.outerRadius*1.496e11)**2 - (self.edge*1.496e11)**2)
self.flux_scalar = self.streamer_scale*self.F_edge/(self.width_edge*self.dS_edge*(self.outerRadius*1.496e11)**self.streamer_powerlaw) #Makes the flux at edge equal to F_edge in the disk
"""
generate a FITS file with name fits_file
"""
def generateFITS(self, disk):
# remove all current files which we will be working with
os.system('rm -Rf images/*')
# remove if the fits file already exists
if os.path.exists('fits/' + self.fits_name):
os.system('rm fits/' + self.fits_name)
#print 'fits/' + self.fits_name, 'was overwritten.'
self.total_flux = 0
# x and y are now transformed into new coordinate plane
# self.center is at (0,0) with points sampled at self.center of each pixel
# initialize an empty 2-D "image" of 0's
image = [[0]*self.width for x in range(self.width)]
# exploit the diagonal symmetry for a sub-2x speedup
for j in xrange(self.width/2):
y_face = (j-self.center+0.5)*self.pxwid
for i in xrange(self.width):
x_face = (i-self.center+0.5)*self.pxwid
# first rotate, then transform for inclination angle
x = x_face*math.cos(self.rotation) + y_face*math.sin(self.rotation)
y = (-x_face*math.sin(self.rotation) + y_face*math.cos(self.rotation)) / math.cos(self.incline)
radius = math.sqrt(x**2+y**2)
# have flux in Jy/m^2, now need in Jy/pixel
flux = disk.calculatePointFlux(radius, 225.5)/math.cos(self.incline)*((self.pxwid*1.496e11)**2)*self.disk_scale
if flux:
# only assign if the flux is NOT 0
image[j][i] = flux
image[-j-1][-i-1] = flux
self.total_flux += flux*2
self.disk_flux = self.total_flux
self.streamer_flux = 0
for j in xrange(self.width):
y_face = (j-self.center+0.5)*self.pxwid*1.496e11
for i in xrange(self.width):
x_face = (i-self.center+0.5)*self.pxwid*1.496e11
#Rotate to plane in which cone is upright. Unlike with the disk above, this is immediately in meters.
#Angle is weird only because self.rotation is defined weird.
x = x_face*math.cos(-math.pi - self.rotation) - y_face*math.sin(-math.pi - self.rotation)
y = x_face*math.sin(-math.pi - self.rotation) + y_face*math.cos(-math.pi - self.rotation)
#Calculate flux due to streamers at (x,y)
radius_cone_in = self.innerRadius*1.496e11 - y/math.tan(self.streamer_angle) #Radius from center of cone. Now given a finite width.
radius_cone_out = self.outerRadius*1.496e11 - y/math.tan(self.streamer_angle)
if y**2 >= self.minoraxis**2*(1-x**2/self.majoraxis**2) and y <= 0:
if abs(x) <= radius_cone_in:
width = 2*(math.sqrt(radius_cone_out**2 - x**2) - math.sqrt(radius_cone_in**2 - x**2))
streamer_flux_density = self.flux_scalar*width*(radius_cone_out)**(self.streamer_powerlaw)
flux = streamer_flux_density*math.sqrt(1+(x**2+radius_cone_out**2/math.tan(self.streamer_angle)**2)/(radius_cone_out**2-x**2))*(self.pxwid*1.496e11)**2
image[j][i] += flux
self.streamer_flux += flux
self.total_flux += flux
elif abs(x) > radius_cone_in and abs(x) <= radius_cone_out:
width = 2*math.sqrt(radius_cone_out**2 - x**2)
streamer_flux_density = self.flux_scalar*width*(radius_cone_out)**(self.streamer_powerlaw) #Equals F_edge when y = 0, r = halfway between inner & outer
flux = streamer_flux_density*math.sqrt(1+(x**2+radius_cone_out**2/math.tan(self.streamer_angle)**2)/(radius_cone_out**2-x**2))*(self.pxwid*1.496e11)**2
image[j][i] += flux
self.streamer_flux += flux
self.total_flux += flux
# convert the Python image array into a numpy 4-D array for MIRIAD compatibility
image = numpy.array([[image]])
# calculate data necessary for header
stardist = disk.getStarDistance()
pix_parsec = self.pxwid/stardist/3600.0 # AU to degrees
# create hdu object to encapsulate data with the header
hdu = pyfits.PrimaryHDU(image)
head = hdu.header
head.update('OBJECT', 'HD-61005')
head.update('SIMPLE', 'T')
head.update('NAXIS', 4)
head.update('NAXIS1', self.width)
head.update('NAXIS2', self.width)
head.update('NAXIS3', 1)
head.update('NAXIS4', 1)
head.update('CDELT1', -1.0*pix_parsec)
head.update('CRPIX1', self.center+0.5)
#head.update('CRVAL1', (7 + 35/60.0 + 47.46192/3600.0)*15) #J2000
head.update('CRVAL1', (7 + 35/60.0 + 47.420/3600.0)*15) #Corrected for PM. I'm lying to MIRIAD because the header in the ppm corrected data file is wrong.
head.update('CTYPE1', 'RA---SIN')
head.update('CDELT2', pix_parsec)
head.update('CRPIX2', center+0.5)
#head.update('CRVAL2', -1.0*(32 + 12/60.0 + 14.0431/3600.0)) #J2000
head.update('CRVAL2', -1.0*(32 + 12/60.0 + 13.32/3600.0)) #Corrected for PM
head.update('CTYPE2', 'DEC--SIN')
head.update('EPOCH', 2000)
head.update('BSCALE', 1.0)
head.update('BZERO', 0.0)
head.update('CDELT3', 8.0e11)
head.update('CRVAL3', -1.0)
head.update('CRPIX3', 1.0)
head.update('CTYPE3', 'STOKES')
head.update('CDELT4', 4.0e9)
head.update('CRPIX4', 1.0)
head.update('CRVAL4', 225.538e9)
head.update('CTYPE4', 'FREQ-LSR')
head.update('BUNIT', 'JY/PIXEL')
head.update('BTYPE', 'INTENSITY')
# write to file
hdu.writeto(self.fits_model_image)
"""
print 'AU^2/pixel:', str(round(self.pxwid**2,4))
print 'total flux (Jy):', total_flux
print 'total flux (Jy*MHz):', str(round(total_flux*225.5e9/1e6,2))
print 'actual SED flux at 225.5GHz:', str(round(disk.calculateFlux(3.0e8/225.5e9)/1e6,2))
"""
"""
generate the visibilities of the model using MIRIAD routines
"""
def generateVisibility(self):
# make sure we have the files we need
if not os.path.exists(self.mir_data_vis):
print self.mir_data_vis, "visibility file does not exist."
exit()
# convert FITS model image into a Miriad image
command('fits in=' + self.fits_model_image + ' op=xyin out=' + self.mir_model_image)
# also copy this into the fits directory just for records
command('cp -R ' + self.fits_model_image + ' fits/' + self.fits_name)
# make the Miriad image into a visibility
command('uvmodel model=' + self.mir_model_image + ' vis=' + self.mir_data_vis
+ ' options=replace out=' + self.mir_model_vis)
# convert the Miriad visibilities into FITS visibilities
command('fits in=' + self.mir_model_vis + ' op=uvout out=' + self.fits_model_vis)
def computeChiSquared(self, disk):
# do all the prereqs
self.generateFITS(disk)
self.generateVisibility()
# get the data from created files
# model = pyfits.getdata(self.fits_model_vis, 'PRIMARY').data
model_hdu = pyfits.open(self.fits_model_vis)
model = model_hdu[0].data.field(5)
# compute the chi squared
# = sum of weight*[(model-data)**2] for real and imaginary
numVis = len(model)
print numVis
# use numpy vector operations
real = model[:,0,0,0,0,0]
imag = model[:,0,0,0,0,1]
chi = float((((real-self.data_vis_real)**2)*self.data_vis_weight).sum())
chi += float((((imag-self.data_vis_imag)**2)*self.data_vis_weight).sum())
# degrees of freedom = 2*[number of visibilities] - [number parameters to fit]
# TEMPO: 4 parameters instead of 6 to fit
#dof = 2*numVis-4
#print 'dof =', dof
chi = chi#/dof
model_hdu.close(closed=True)
return chi
def changeParameters(self, disk, disk_scale, streamer_scale):
#Disk Stuff
self.innerRadius = disk.getInnerRadius()
self.outerRadius = disk.getOuterRadius()
self.majoraxis = self.outerRadius*1.496e11
self.minoraxis = self.majoraxis*math.sin(math.pi/2 - self.incline)
#Scaling Stuff
self.disk_scale = disk_scale
self.streamer_scale = streamer_scale
self.edge = (self.innerRadius + self.outerRadius)/2 #In AU
self.F_edge = disk.calculatePointFlux(self.edge, 225.5)/math.cos(self.incline)
self.dS_edge = math.sqrt(1+((self.edge*1.496e11)**2+(self.outerRadius*1.496e11)**2/math.tan(self.streamer_angle)**2)/((self.outerRadius*1.496e11)**2-(self.edge*1.496e11)**2))
self.width_edge = 2*math.sqrt((self.outerRadius*1.496e11)**2 - (self.edge*1.496e11)**2)
self.flux_scalar = self.streamer_scale*self.F_edge/(self.width_edge*self.dS_edge*(self.outerRadius*1.496e11)**self.streamer_powerlaw) #Makes the flux at edge equal to F_edge in the disk
def getStreamerFlux(self):
streamer_flux = 0
for j in xrange(self.width):
y_face = (j-self.center+0.5)*self.pxwid*1.496e11
for i in xrange(self.width):
x_face = (i-self.center+0.5)*self.pxwid*1.496e11
#Rotate to plane in which cone is upright. Unlike with the disk above, this is immediately in meters.
#Angle is weird only because self.rotation is defined weird.
x = x_face*math.cos(-math.pi - self.rotation) - y_face*math.sin(-math.pi - self.rotation)
y = x_face*math.sin(-math.pi - self.rotation) + y_face*math.cos(-math.pi - self.rotation)
#Calculate flux due to streamers at (x,y)
radius_cone_in = self.innerRadius*1.496e11 - y/math.tan(self.streamer_angle) #Radius from center of cone. Now given a finite width.
radius_cone_out = self.outerRadius*1.496e11 - y/math.tan(self.streamer_angle)
if y**2 >= self.minoraxis**2*(1-x**2/self.majoraxis**2) and y <= 0:
if abs(x) <= radius_cone_in:
width = 2*(math.sqrt(radius_cone_out**2 - x**2) - math.sqrt(radius_cone_in**2 - x**2))
streamer_flux_density = self.flux_scalar*width*(radius_cone_out)**(self.streamer_powerlaw)
flux = streamer_flux_density*math.sqrt(1+(x**2+radius_cone_out**2/math.tan(self.streamer_angle)**2)/(radius_cone_out**2-x**2))*(self.pxwid*1.496e11)**2
streamer_flux += flux
elif abs(x) > radius_cone_in and abs(x) <= radius_cone_out:
width = 2*math.sqrt(radius_cone_out**2 - x**2)
streamer_flux_density = self.flux_scalar*width*(radius_cone_out)**(self.streamer_powerlaw) #Equals F_edge when y = 0, r = halfway between inner & outer
flux = streamer_flux_density*math.sqrt(1+(x**2+radius_cone_out**2/math.tan(self.streamer_angle)**2)/(radius_cone_out**2-x**2))*(self.pxwid*1.496e11)**2
streamer_flux += flux
return streamer_flux
def getTotalFlux(self):
return self.total_flux, self.disk_flux
def command(cmd):
os.system(cmd + ' > /dev/null 2>&1')