-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgetxyuvtwoparab_apper.py
153 lines (144 loc) · 4.8 KB
/
getxyuvtwoparab_apper.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
import numpy as np
import math
def getxyuveqtwoparabolas(initdat,isper):
phi,orbitangle,orbitalradius,eccentricity, masses=initdat
print("initdat", orbitalradius, phi, eccentricity, np.cos(phi), np.sin(phi))
metersperAU=1
Gconstant=1
#fix x0 y0 at one star, disregard initial data, use orbital radius as separation between stars
#this is consistant with choice in previous part
cosphi=np.cos(phi)
sinphi=np.sin(phi)
coordsep=orbitalradius #/2.
print("coordsep",coordsep)
#x0=orbitalradius/2.*np.cos(phi)*metersperAU
#y0=orbitalradius/2.*np.sin(phi)*metersperAU
count=0
for phi0 in phi:
if phi0==0:
print("zero")
#x0[count]=orbitalradius[count]/2.
#y0[count]=0
cosphi[count]=1.0
sinphi[count]=0.0
if phi0==math.pi:
cosphi[count]=-1.0
sinphi[count]=0.0
print("pi")
#x0[count]=-orbitalradius[count]/2.
#y0[count]=0
if phi0==math.pi/2.:
cosphi[count]=0.0
sinphi[count]=1.0
#x0[count]=0
#y0[count]=orbitalradius[count]/2.
if phi0==3*math.pi/2.:
cosphi[count]=0.0
sinphi[count]=-1.0
#x0[count]=0
#y0[count]=-orbitalradius[count]/2.
count+=1
x0=coordsep*cosphi
y0=coordsep*sinphi
#x0[1]=0.0
#y0[1]=0.0
z0=np.zeros(2)
print(x0)
print(y0)
v=np.zeros(2)
a=np.zeros(2)
ux0=0.
uy0=0.
uz0=0.
ax0=0.
ay0=0.
az0=0.
vx0=0.
vy0=0.
vz0=0.
if eccentricity==1:
#elliptical
orbitalr=orbitalradius #NO REDUCED MASS
coordsep=2*orbitalr #start at aphelion
isper=True
if isper:
coordsep=orbitalr
x0=(coordsep)*cosphi
y0=(coordsep)*sinphi
starsep=np.sqrt((x0[0]-x0[1])**2+(y0[0]-y0[1])**2) #two stars, at opposite ends of the orbit
Fapastron=masses[1]*masses[0]/starsep**2
v=np.zeros(2)
#v= np.sqrt(masses[1]*masses[0]/masses*(2./starsep-1./(2*orbitalr)))
vapsq=masses[1]*masses[0]/masses/orbitalr/2.
if isper:
vapsq=masses[1]*masses[0]/masses/orbitalr/2.
#vapsq=masses[1]*masses[0]/masses*(1./(coordsep)-1./orbitalr)/2.
v=np.sqrt(vapsq)
#v= np.sqrt(masses[1]*masses[0]/masses*(1./orbitalr-2./starsep))
ux0=-v*sinphi
uy0=v*cosphi #initial data in y only
uz0=np.zeros(2)
a=np.zeros(2)
a=Fapastron/masses
ax0=-a*cosphi
ay0=-a*sinphi
az0=np.zeros(2)
elif eccentricity<1 and eccentricity>0:
#elliptical
orbitalr=orbitalradius #NO REDUCED MASS
coordsep=orbitalr*(1+eccentricity) #start at aphelion
if isper:
coordsep=orbitalr*(1-eccentricity)
x0=(coordsep)*cosphi
y0=(coordsep)*sinphi
starsep=np.sqrt((x0[0]-x0[1])**2+(y0[0]-y0[1])**2) #two stars, at opposite ends of the orbit
Fapastron=masses[1]*masses[0]/starsep**2
v=np.zeros(2)
#v= np.sqrt(masses[1]*masses[0]/masses*(2./starsep-1./(2*orbitalr)))
vapsq=masses[1]*masses[0]/masses/orbitalr*(1-eccentricity)/(1+eccentricity)/4.
if isper:
vapsq=masses[1]*masses[0]/masses/orbitalr*(1+eccentricity)/(1-eccentricity)/4.
#vapsq=masses[1]*masses[0]/masses*(1./(coordsep)-1./orbitalr)/2.
v=np.sqrt(vapsq)
#v= np.sqrt(masses[1]*masses[0]/masses*(1./orbitalr-2./starsep))
ux0=-v*sinphi
uy0=v*cosphi #initial data in y only
uz0=np.zeros(2)
a=np.zeros(2)
a=Fapastron/masses
ax0=-a*cosphi
ay0=-a*sinphi
az0=np.zeros(2)
elif eccentricity[0]==0.0: #circular
#start at perihelion for both (eliptical, doesn't generalize to three body)
#actually start with circular orbit
ux0=np.zeros(2) #*149597870700
#centrepital force balances gravitational force
metersperAU=1 #natural units
#G=1
Gconstant=1
r0=2.*orbitalradius #Mystery factor of 2
print("r0", r0)
v=np.zeros(2)
for i in np.arange(2):
v[i]=np.sqrt(Gconstant*masses[(i+1)%2]/np.abs(r0[i]))
print(v)
#r0=orbitalradius #np.sqrt(x0**2+y0**2)
ux0=-v[0]*sinphi
uy0=v[1]*cosphi #initial data in y only
uz0=np.zeros(2)
a=np.zeros(2)
a=Fapastron/masses
ax0=-a[0]*cosphi
ay0=-a[1]*sinphi
az0=np.zeros(2)
statevec=[]
avec=[]
for i in np.arange(len(x0)):
stateveci=np.array([x0[i],y0[i],z0[i],ux0[i],uy0[i],uz0[i]])
aveci=np.array([ax0[i],ay0[i],az0[i]])
statevec.append(stateveci)
avec.append(aveci)
statevecnp=np.array(statevec)
avecnp=np.array(avec)
return masses, statevecnp,avecnp