Ast1100/orbitalcalculator.py at master · dulte/Ast1100 · GitHub

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
import numpy as np
import matplotlib.pyplot as plt
from AST1100SolarSystem import AST1100SolarSystem
import scipy.interpolate as inter
from numpy.linalg import norm
import time as tid
import seaborn
import sys


def calc_r(pos,time,posfunc):
    return norm(pos-posfunc(time)[:,1])

def calc_r_vec(pos,time,posfunc):
    return (pos-posfunc(time)[:,1])

def  calc_correction_burn(pos,tim,posfunc):
    r = (pos-posfunc(time+0.001)[:,1])
    r_norm = 1.0/(norm((pos-posfunc(time)[:,1])))*r
    return 100*r

def calc_planet_vel(time,func):
    h = 1.0/300000
    return (func(time+h)[:,1] - func(time)[:,1])/(h)

def calc_injection_burn(time,func,pos,planet_vel,sat_vel,p_mass):
    G = 4*np.pi**2

    r = pos - func(time)[:,1]

    #print norm(r)
    theta = np.arctan2(r[1],r[0])
    orb_vel = np.sqrt(G*p_mass/norm(r))
    #print "orb_vel: ",orb_vel
    #print "Vel diff",planet_vel-sat_vel

    return np.array([-orb_vel*np.sin(theta),orb_vel*np.cos(theta)]) - sat_vel + planet_vel

def print_theta(func,pos,time):
    theta_p = np.arctan2(func(time)[1,1],func(time)[0,1])
    theta_s = np.arctan2(pos[1],pos[0])

    print "Theta: ",theta_p - theta_s


seed = 75041
sys = AST1100SolarSystem(seed)

km_to_au = 6.685e-9

with open("planetPositions.npy", "rb") as npy:
    [planetPositions,t] = np.load(npy)

steps =  planetPositions[0,0,:].size
print steps
#t = np.linspace(0,20,steps)

planetPosFunction = inter.interp1d(t, planetPositions)

# pos = np.array([0.80253929065 , 4.16457524661])
# time = 8.81981666666
#
# print calc_r(pos,time,planetPosFunction)
# print calc_r_vec(pos,time,planetPosFunction)
# print calc_correction_burn(pos,time,planetPosFunction)
# print norm(calc_correction_burn(pos,time,planetPosFunction))
#
#
#
# pos2 = np.array([0.792907737234 , 4.16748703398])
# time2 = 8.82231666666
#
# print "-----"
# print calc_r_vec(pos2,time2,planetPosFunction)
# print "When burning: ",calc_r(pos2,time2,planetPosFunction)
# print "-----"
#
# vel_inj = np.array([-3.77504600841 , 1.09819863799])
# pos_inj = np.array([0.792907737234 , 4.16748703398])
# time_inj = 8.82231666666
#
# planet_vel = calc_planet_vel(time,planetPosFunction)
# p_mass = sys.mass[1]
#
#
# print "injection burn: " ,calc_injection_burn(time_inj,planetPosFunction,pos_inj,planet_vel,vel_inj,p_mass)
# #print norm(calc_injection_burn(time_inj,planetPosFunction,pos_inj,planet_vel,vel_inj,p_mass))
# print "-----"
#
# time = 10.82231666666
# pos = np.array([-4.1747684284 ,  -1.02819079768])
#
# print "After 2 years: ",calc_r(pos,time,planetPosFunction)
# print_theta(planetPosFunction,pos,time)


pos = np.array([0.850583057699 , 4.15384934247])
time = 8.80881666666

print calc_r(pos,time,planetPosFunction)
print calc_correction_burn(pos,time,planetPosFunction)
print norm(calc_correction_burn(pos,time,planetPosFunction))
print "--------"

pos = np.array([0.833365798449 , 4.15916091718])
time = 8.81281666666

print calc_r(pos,time,planetPosFunction)
print "------"


vel_inj = np.array([-4.33874216505 , 1.33664064517])
pos_inj = np.array([0.833361459735 , 4.15916225383])
time_inj = 8.81281666666

planet_vel = calc_planet_vel(time,planetPosFunction)
p_mass = sys.mass[1]


print "injection burn: " ,calc_injection_burn(time_inj,planetPosFunction,pos_inj,planet_vel,vel_inj,p_mass)
print norm(calc_injection_burn(time_inj,planetPosFunction,pos_inj,planet_vel,vel_inj,p_mass))
print "-----"


time = 10.81281666666
pos = np.array([-4.18412274402 , -0.989071568075])

print "After 2 years: ",calc_r(pos,time,planetPosFunction)
print_theta(planetPosFunction,pos,time)