AST1100-SolarSystem/simContinueFineTuning.py at master · halvarsu/AST1100-SolarSystem · 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
"""
File: simContinueFineTuning.py
Author: Halvard Sutterud
Email: halvard.sutterud@gmail.com
Github: https://github.com/halvarsu
Description: Used for finding the leg from gravity assist planet to target
planet, with different values for boost
"""

from MySateliteSim import MySateliteSim
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d


def rot_matrix(angle):
    c = np.cos(angle)
    s = np.sin(angle)
    return np.array(((c,-s),(s,c)))

def boost(velocity, boost, boost_angle):
    direction = velocity/np.linalg.norm(velocity)
    boost_direction = rot_matrix(boost_angle).dot(direction)

    new_vel = vCl + boost_direction*boost/system.au*system.year
    return new_vel


system = MySateliteSim(87464)
au = system.au
launch_v = 9594

p = np.load('positionsHomePlanet.npy')
t = np.load('times.npy')
st = np.load('data/times%d.npy' %launch_v)
s = np.load('data/pos%d.npy' % launch_v)

posFunc = interp1d(t[:-1], p)
planetPos = posFunc(st)

satFunc = interp1d(st, s.T)
def vel(t,dt=0.0000001):
    return (satFunc(t+dt)-satFunc(t-dt))/(2*dt)

planet4 = planetPos[:,4]
relPos = s.T-planet4
relDist = np.linalg.norm(relPos, axis =0)
closest_dist = np.amin(relDist[:-10])
closest_i = np.argmin(relDist[:-10])
relposCl = relPos[:,closest_i]
hCl = closest_dist - system.radius[4]*1000./au
tCl = st[closest_i]
pCl = s[closest_i]
aCl = np.arctan(relposCl[1]/relposCl[0])#angle
vCl = vel(tCl)

#if hCl < 0:
    #print "Wont bother with values inside planet"
    #continue

print "time of closest approach: ", tCl
print "satpos at closest approach: ", pCl
print "distance from surface at closest approach: ", hCl
print "relpos at closest approach", relposCl
print "angle at closest approach", aCl

boost_v = 1480

boost_a = 1.081*np.pi/4
vel0 = boost(vCl, boost_v, boost_a)
pos0 = pCl
sim_length = 6.2
close_speed_factor = 2

# 0.336138702386  b1000 a1.06pi/4 lv 9594
# 0.0596081065397 b2000
# 0.00503893068436 a1.07pi/4
# 0.0268986034289 a1.065pi/4 v0.629758507124
# 0.00576089852027 b1500 a1.075pi/4 ?
# 0.0026280081882 b1500 a1.079pi/4
# 0.000108329130473 b1480 a1.081pi/4

#Trying with larger dt, want to find better boost angle

dt_far = 1/(50000.)
print "---------------------------------------------"
print "continuing the launch ", launch_v
psim, vsim, asim, tsim = system.satelite_sim(t0 = tCl,
        tN = tCl+sim_length, start_planet=4, target_planet =5,
        vel0 = vel0, pos0 = pos0, first_boost = False,
        dt_far=dt_far, sec_per_dt_close = close_speed_factor)
save = True
if close_speed_factor > 1:
    if raw_input("Sure you want to save?(y/N)") == 'y':
        save = True
    else:
        save = False
if save:
    filenames = ['%sv%db%d' %(name,launch_v,boost_v) for name in ('pos','vel','times')]
    system.saveData(fname = filenames, folder = 'postGAdata')
system.plot()

#print close_list
#rang = [1990 - 10*x for x in range(5)]
#plt.scatter(close_list, rang)
#plt.show()