-
Notifications
You must be signed in to change notification settings - Fork 0
/
gk_ss_gaussian.py
136 lines (112 loc) · 4.25 KB
/
gk_ss_gaussian.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
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Nov 12 13:08:06 2016
@author: john
"""
import numpy as np
import matplotlib.pyplot as plt
from math import tan, cos
from astropy.io import fits
from galpy.orbit import Orbit
from galpy.potential import MWPotential2014
from galpy.actionAngle import actionAngleAdiabatic
from galpy.actionAngle import estimateDeltaStaeckel
from galpy.actionAngle import actionAngleStaeckel
from galpy.potential import KeplerPotential, SteadyLogSpiralPotential
from galpy.util import bovy_conversion
from astropy import units
import pickle as cp
fd = open("data/dbscans/gk_ss_cands.pkl", "rb")
gk_candidates = cp.load(fd)
fd.close()
gk_candidates = np.array(gk_candidates)
def generate_rand(measurement, error, measurement_type = None):
#the type is the important part for converting the measurement
#the reason for this function is to handle measurement confusion
if measurement_type == 'mas':
error = error/3600000
value = np.random.normal(loc = measurement, scale = error)
return value
def initialise_orbit(ID):
cand = gk_candidates[(gk_candidates.T[-1] == int(ID))]
cand = cand[0]
ra_mean = cand[0]
ra_err = 0 #unavailable
dec_mean = cand[1]
dec_err = 0
pmRA_mean = cand[3]
pmRA_err = cand[6]
pmDec_mean = cand[4]
pmDec_err = cand[7]
#calc a random value from gaussian
# ra = generate_rand(ra_mean, ra_err, measurement_type = 'mas')
# dec = generate_rand(dec_mean, dec_err, measurement_type = 'mas')
ra = ra_mean
dec = dec_mean
pmRA = generate_rand(pmRA_mean, pmRA_err)
if pmDec_mean == pmDec_err: #some issue with things being wrong just ignore it
pmDec = pmDec_mean
else:
pmDec = generate_rand(pmDec_mean, pmDec_err)
#no error for distance in this set
dist = cand[2]
#VLOS is reliable and stays as is
Vlos = cand[5]
orbit = Orbit(vxvv=[ra,dec,dist,pmRA,pmDec,Vlos],radec=True,ro=8.,vo=220., solarmotion = "schoenrich")
return orbit
def add_spiral(arms):
#spiral parameters and defaults then generate a value for it
sp_m = arms #4
sp_spv = (22.5, 7.5) #20 15-30
sp_spv = generate_rand(sp_spv[0], sp_spv[1])
if sp_m == 4:
sp_i = (-12*(3.1415/180), 2*(3.1415/180)) #12, err = 2 6, err = 1 (2 arms)
sp_i = generate_rand(sp_i[0], sp_i[1])
if sp_m == 2:
sp_i = (-6*(3.1415/180), 1*(3.1415/180)) #12, err = 2 6, err = 1 (2 arms)
sp_i = generate_rand(sp_i[0], sp_i[1])
sp_x_0 = (-120*(3.1415/180), 10*(3.1415/180)) #-120 err = 10
sp_x_0 = generate_rand(sp_x_0[0], sp_x_0[1])
sp_a = -sp_m/tan(sp_i)
sp_gamma = sp_x_0/sp_m
sp_fr_0 = (0.05, 0.01)
sp_fr_0 = generate_rand(sp_fr_0[0], sp_fr_0[1])
# ro = 8
# vo = 220
sp_A = -(1)*sp_fr_0 #*((ro*vo)**2) #not sure on this part but adding these makes it too big so leave them
sp_component = SteadyLogSpiralPotential(amp = 1, omegas = sp_spv/220, A = sp_A, alpha = sp_a, gamma = sp_gamma)
mwp_2D.append(sp_component)
ts = np.linspace(0,-150,10000)
mwp = MWPotential2014
sun = Orbit(vxvv=[0, 0, 0, 0, 0, 0],radec=True,ro=8.,vo=220., solarmotion = "schoenrich")
sun.integrate(ts,mwp,method='odeint')
mwp_2D = [i.toPlanar() for i in mwp]
sun_2D = sun.toPlanar()
sun_2D.integrate(ts,mwp_2D,method='odeint')
runs = 25
counts = []
for sobject in gk_candidates.T[-1][4:]:
print("finally calculating distance from sun for 50 iterations " + str(int(sobject)))
count = 0
# plots = 0
arms = 2
for i in range(runs):
add_spiral(arms)
o = initialise_orbit(sobject)
o_2D = o.toPlanar()
o_2D.integrate(ts,mwp_2D,method='odeint')
delta_x = o_2D.x(ts) - sun_2D.x(ts)
delta_y = o_2D.y(ts) - sun_2D.y(ts)
delta = (delta_x**2 + delta_y**2)**0.5
delta = delta*1000
# plt.semilogy(sun_2D.time(ts), delta*1000)
delta_before_4gyr = delta[(sun_2D.time(ts) <= -4)][(sun_2D.time(ts)[(sun_2D.time(ts) <= -4)] >= -5.2)]
if delta_before_4gyr.min() <= 100:
count += 1
# if plots <= 3:
# plt.semilogy(sun.time(ts), delta)
# plots += 1
counts.append(count)
break
print(runs, arms, counts)