frogleap ok but not hermite
This commit is contained in:
@@ -3,7 +3,9 @@
|
|||||||
"""
|
"""
|
||||||
Class definition for physical atribute
|
Class definition for physical atribute
|
||||||
"""
|
"""
|
||||||
|
from os import system
|
||||||
import numpy as np
|
import numpy as np
|
||||||
|
from lib.plots import DynamicUpdate
|
||||||
|
|
||||||
globals()['G'] = 6.67e-11 #Gravitational constant in SI units
|
globals()['G'] = 6.67e-11 #Gravitational constant in SI units
|
||||||
globals()['Ms'] = 2e30 #Solar mass in kg
|
globals()['Ms'] = 2e30 #Solar mass in kg
|
||||||
@@ -150,15 +152,38 @@ class System:
|
|||||||
return 1
|
return 1
|
||||||
|
|
||||||
def HPC(self, dt): # update position and velocities of bodies in system with hermite predictor corrector
|
def HPC(self, dt): # update position and velocities of bodies in system with hermite predictor corrector
|
||||||
self.update_a()
|
self.COMShift()
|
||||||
self.update_j()
|
self.Update_a()
|
||||||
self.predict(dt)
|
self.Update_j()
|
||||||
self.update_ap()
|
self.Predict(dt)
|
||||||
self.update_jp()
|
self.Update_ap()
|
||||||
self.update(dt)
|
self.Update_jp()
|
||||||
|
self.Correct(dt)
|
||||||
self.time = self.time + dt
|
self.time = self.time + dt
|
||||||
for body in self.bodylist:
|
for body in self.bodylist:
|
||||||
body.p = body.v*body.m
|
body.p = body.v*body.m
|
||||||
|
|
||||||
|
def hermite(self, duration, dt, display=False):
|
||||||
|
if display:
|
||||||
|
try:
|
||||||
|
system("mkdir tmp")
|
||||||
|
except IOError:
|
||||||
|
system("rm tmp/*")
|
||||||
|
d = DynamicUpdate()
|
||||||
|
d.on_launch()
|
||||||
|
|
||||||
|
N = np.ceil(duration/dt).astype(int)
|
||||||
|
for j in range(N):
|
||||||
|
self.HPC(dt)
|
||||||
|
|
||||||
|
if display:
|
||||||
|
# display progression
|
||||||
|
q_array = self.get_positions()
|
||||||
|
if len(self.bodylist) == 1:
|
||||||
|
d.on_running(q_array[0], q_array[1], q_array[2], step=j, label="step {0:d}/{1:d}".format(j,N))
|
||||||
|
else:
|
||||||
|
d.on_running(q_array[:,0], q_array[:,1], q_array[:,2], step=j, label="step {0:d}/{1:d}".format(j,N))
|
||||||
|
#time.sleep(1e-5)
|
||||||
|
|
||||||
return 1
|
return 1
|
||||||
|
|
||||||
|
|||||||
3
main.py
3
main.py
@@ -33,8 +33,9 @@ def main():
|
|||||||
dyn_syst = System(bodylist)
|
dyn_syst = System(bodylist)
|
||||||
dyn_syst.COMShift()
|
dyn_syst.COMShift()
|
||||||
|
|
||||||
duration, step = 100*3e7, 1e5
|
duration, step = 100*3e7, 1e4
|
||||||
E, L = frogleap(duration, step, dyn_syst, recover_param=True)#, display=True)
|
E, L = frogleap(duration, step, dyn_syst, recover_param=True)#, display=True)
|
||||||
|
#dyn_syst.hermite(duration,step, display=True)
|
||||||
|
|
||||||
fig1 = plt.figure(figsize=(30,15))
|
fig1 = plt.figure(figsize=(30,15))
|
||||||
ax1 = fig1.add_subplot(111)
|
ax1 = fig1.add_subplot(111)
|
||||||
|
|||||||
BIN
plots/Em.png
Normal file
BIN
plots/Em.png
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 295 KiB |
BIN
plots/L2.png
Normal file
BIN
plots/L2.png
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 126 KiB |
Reference in New Issue
Block a user