This repository has been archived by the owner on Mar 5, 2020. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathbindata_keplerian.py
113 lines (93 loc) · 4.09 KB
/
bindata_keplerian.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
#bindata_keplerian.py
import csv
import numpy as np
import scipy as sp
import astropy as ap
from astropy import units as u
from astropy.coordinates import Angle
def calculate_delay(r_angle, r_phase, r_rad, timescale):
"""Delay relative to continuum observed at angle for emission at radius"""
# Calculate delay compared to continuum photons
# Draw plane at r_rad_min out. Find x projection of disk position.
# Calculate distance travelled to that plane from the current disk position
# Delay relative to continuum is thus (distance from centre to plane)
# + distance from centre to point
vr_disk = np.array([r_rad.value*np.cos(r_phase), 0.0]) * u.m
vr_normal = np.array([np.cos(r_angle), np.sin(r_angle)])
vr_plane = r_rad * vr_normal
return (np.dot((vr_plane - vr_disk), vr_normal) * u.m / ap.constants.c).to(timescale)
def keplerian_velocity(r_mass, r_rad):
"""Calculates Keplerian velocity at radius"""
return np.sqrt(ap.constants.G * r_mass / r_rad)
def path_to_delay(r_path):
"""Converts path to time delay"""
return r_path / ap.constants.C
def doppler_shift(r_wave, r_vel):
"""Doppler shifts the passed wavelength"""
return r_wave * (1 + (r_vel / ap.constants.c))
# INPUTS
#r_line_inp = 1550.7
r_line_inp = 6562.8
r_rad_min_inp = 50
r_rad_max_inp = 100
#r_mass_bh_inp = 1e9
r_mass_bh_inp = 1e7
r_angle_inp = 40
i_steps = 100
# PROCESSING
print "=== STARTING ==="
print "--- PARSE INPUT ---"
r_line = r_line_inp * u.angstrom
r_angle = Angle(r_angle_inp, u.deg)
r_mass_bh = r_mass_bh_inp * ap.constants.M_sun
r_rad_grav = (6 * ap.constants.G * r_mass_bh / np.power(ap.constants.c, 2))
r_rad_min = r_rad_min_inp * r_rad_grav
r_rad_max = r_rad_max_inp * r_rad_grav
print "--- SETUP ARRAYS ---"
ar_wave = np.zeros(i_steps) * u.angstrom
ar_delay = np.zeros(i_steps) * u.s
ar_phase = np.linspace(0, np.pi*2, i_steps)
ar_rad = np.linspace(r_rad_min,r_rad_max*20.0,i_steps)
print "--- SETUP OUTPUT ---"
csvfile = open('bindata_keplerian.csv', 'w')
csvout = csv.writer(csvfile, delimiter=' ', quotechar='"', quoting=csv.QUOTE_MINIMAL)
as_row_titles = ["Wavelength", "Delay"]
csvout.writerow(as_row_titles)
# ITERATE OVER INNER EDGE
r_vel_max = keplerian_velocity( r_mass_bh, r_rad_min )
print "Wave: %.2e, Vel: %.2e" % (doppler_shift(r_line, r_vel_max).value, r_vel_max.value)
print "Wave: %.2e, Vel: %.2e" % (doppler_shift(r_line, r_vel_max*np.sin(r_angle)).value, r_vel_max.value*np.sin(r_angle))
for r_phase, r_wave, r_delay in np.nditer([ar_phase, ar_wave, ar_delay], op_flags=['readwrite']):
r_vel = r_vel_max * np.sin(r_phase) * np.sin(r_angle)
r_wave[...] = doppler_shift( r_line, r_vel )
r_delay[...] = calculate_delay( r_angle, r_phase, r_rad_min, u.day)
for r_wave, r_delay in np.nditer([ar_wave, ar_delay]):
csvout.writerow((r_wave, r_delay))
csvout.writerow([])
# # ITERATE OVER OUTER EDGE
# r_vel_max = keplerian_velocity( r_mass_bh, r_rad_max )
# for r_phase, r_wave, r_delay in np.nditer([ar_phase, ar_wave, ar_delay], op_flags=['readwrite']):
# r_vel = r_vel_max * np.sin(r_phase) * np.sin(r_angle)
# r_wave[...] = doppler_shift( r_line, r_vel )
# r_delay[...] = calculate_delay( r_angle, r_phase, r_rad_max, u.day)
# for r_wave, r_delay in np.nditer([ar_wave, ar_delay]):
# csvout.writerow((r_wave, r_delay))
# csvout.writerow([])
# ITERATE OVER BLUE BOUND
for r_rad, r_wave, r_delay in np.nditer([ar_rad, ar_wave, ar_delay], op_flags=['readwrite']):
r_rad = r_rad * u.m
r_vel = keplerian_velocity( r_mass_bh, r_rad ) * np.sin(r_angle)
r_wave[...] = doppler_shift( r_line, r_vel )
r_delay[...] = calculate_delay( r_angle, np.pi/2, r_rad, u.day )
for r_wave, r_delay in np.nditer([ar_wave, ar_delay]):
csvout.writerow((r_wave, r_delay))
csvout.writerow([])
# ITERATE OVER RED BOUND
for r_rad, r_wave, r_delay in np.nditer([ar_rad, ar_wave, ar_delay], op_flags=['readwrite']):
r_rad = r_rad * u.m
r_vel = -keplerian_velocity( r_mass_bh, r_rad ) * np.sin(r_angle)
r_wave[...] = doppler_shift( r_line, r_vel )
r_delay[...] = calculate_delay( r_angle, np.pi/2, r_rad, u.day )
for r_wave, r_delay in np.nditer([ar_wave, ar_delay]):
csvout.writerow((r_wave, r_delay))
csvout.writerow([])