-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathss_apo_gaussian.py
196 lines (170 loc) · 6.6 KB
/
ss_apo_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
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
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Nov 2 18:19:29 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
from astropy.io import ascii
try:
fd = open("data/arms_apogee_data.pkl", "rb")
arms_apogee_data = cp.load(fd)
fd.close()
print("ids, vlos, and ucac data loaded from file")
except:
candidate_ids = ['2M06150356-0021091', '2M04221969+4439373', '2M06434070+0051170']
vlos_dict = [[candidate_ids[0], 11.864100456237793], [candidate_ids[1], 46.98500061035156], [candidate_ids[2], 26.171199798583984]]
vlos_dict = dict(vlos_dict)
#grab the ages
#solar_age = 4.6*(10**9) #try and grab an astropy constant here
#age_table = open('data/apo_ages.dat')
#age_and_mass = ascii.read(age_table.readlines())
##here we have
#ids_with_age = "col2"
#ln_mass = "col7"
#ln_age = "col8"
#ln_mass_e = "col15"
#ln_age_e = "col16"
#inspect them if we want
#for ID in candidate_ids:
# for i in age_and_mass:
# if ID == i[1]:
# print(ID, ":", i[8])
#grab the distances for these candidates
hdudist = fits.open('data/allStar+-v603.150209.fits')
distances = hdudist[1].data.field('DISO').T[1].T
ap_id = hdudist[1].data.field('APOGEE_ID')
dist_table = np.vstack([ap_id, distances]).T
candidate_dist = []
for ID in candidate_ids:
for i in dist_table:
if ID == i[0]:
candidate_dist.append(i)
candidate_dist = dict(candidate_dist)
#grab the UCACA4 values
ucac = np.loadtxt('data/ss_ucac4.csv', dtype = str, delimiter=',')
#clean it up
for i in range(len(ucac)):
for e in range(len(ucac[i])):
ucac[i][e] = ucac[i][e][2:-1]
#make a dcit with apo_id as key and ucac values as a list
ucac_dict = []
for i in ucac:
ucac_dict.append([i[0], i[1:]])
ucac_dict = dict(ucac_dict)
arms_apogee_data = [candidate_ids, vlos_dict, ucac_dict, candidate_dist]
fd = open("data/arms_apogee_data.pkl", "wb")
cp.dump(arms_apogee_data, fd)
fd.close()
print("saved some shit")
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
#arm_apogee_data = [candidate_ids, vnolos_dict, ucac_dict]
candidate_ids = arms_apogee_data[0]
vlos_dict = arms_apogee_data[1]
ucac_dict = arms_apogee_data[2]
candidate_dist = arms_apogee_data[3]
"""
ucac4 values work as follows with apo id as the key
ucac_dict = [ucac4_id, ra, ra_err, dec, dec_err, pmra, pmra_err, pmdec, pmdec_err]
"""
print("initialising 3D orbits for solar siblings...")
#get orbits for each of them with schoenrich
def initialise_orbit(ID):
ra_mean = float(ucac_dict[ID][1])
ra_err = float(ucac_dict[ID][2])
dec_mean = float(ucac_dict[ID][3])
dec_err = float(ucac_dict[ID][4])
pmRA_mean = float(ucac_dict[ID][5])
pmRA_err = float(ucac_dict[ID][6])
pmDec_mean = float(ucac_dict[ID][7])
pmDec_err = float(ucac_dict[ID][8])
#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')
pmRA = generate_rand(pmRA_mean, pmRA_err)
pmDec = generate_rand(pmDec_mean, pmDec_err)
#no error for distance in this set
dist = float(candidate_dist[ID])
#VLOS is reliable and stays as is
Vlos = float(vlos_dict[ID])
orbit = Orbit(vxvv=[ra,dec,dist,pmRA,pmDec,Vlos],radec=True,ro=8.,vo=220., solarmotion = "schoenrich")
return orbit
#sun's orbit for comparison:
print("Preparing data for 2D analysis")
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')
print("Flattening orbits and potential")
mwp_2D = [i.toPlanar() for i in mwp]
sun_2D = sun.toPlanar()
sun_2D.integrate(ts,mwp_2D,method='odeint')
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 = -sp_fr_0 #to get from ramirez to galpy divide by 220 below so its in natural units
sp_component = SteadyLogSpiralPotential(amp = 1, omegas = sp_spv/220, A = sp_A, alpha = sp_a, gamma = sp_gamma)
mwp_2D.append(sp_component) # not yet
run = 25
counts = []
for ID in candidate_ids[0:]:
print("finally calculating distance from sun for 10 iterations " + ID)
count = 0
arms = 2
for i in range(run):
# add_spiral(arms)
o = initialise_orbit(ID)
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)
break
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
plt.semilogy(sun_2D.time(ts), delta)
break
counts.append(count)
break
print(run, arms, counts)
#plt.title("Solar Sibling Proximity")
#plt.xlabel("Time (Gyr)")
#plt.ylabel("Distance (pc)")
#plt.rcParams.update({'font.size': 16})