-
Notifications
You must be signed in to change notification settings - Fork 55
/
Copy pathmakeplot.py
119 lines (99 loc) · 3.59 KB
/
makeplot.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
import pathlib
import sys
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import vplot
import vplanet
# Path hacks
path = pathlib.Path(__file__).parents[0].absolute()
sys.path.insert(1, str(path.parents[0]))
from get_args import get_args
# Run vplanet
rd4 = vplanet.run(path / "RD4" / "vpl.in")
ll2 = vplanet.run(path / "LL2" / "vpl.in")
time = rd4.TGstar.Time / 1e3
#################### 4th order DistOrb ####################
fig = plt.figure(figsize=(6.5, 8))
plt.subplot(3, 2, 1)
plt.plot(time, rd4.TGb.Obliquity, color=vplot.colors.dark_blue)
plt.plot(time, rd4.TGc.Obliquity, color=vplot.colors.red)
plt.xlabel("Time (kyr)")
plt.ylabel(r"Obliquity ($^\circ$)")
plt.legend('', frameon=False) # Why do we need these!??
plt.subplot(3, 2, 2)
plt.plot(time, rd4.TGb.Eccentricity, color=vplot.colors.dark_blue)
plt.plot(time, rd4.TGc.Eccentricity, color=vplot.colors.red)
plt.xlabel("Time (kyr)")
plt.ylabel("Eccentricity")
plt.legend('', frameon=False)
plt.subplot(3, 2, 3)
plt.plot(time, rd4.TGb.RotPer, color=vplot.colors.dark_blue)
plt.plot(time, rd4.TGc.RotPer, color=vplot.colors.red)
plt.xlabel("Time (kyr)")
plt.ylabel("Rotation Period (days)")
plt.legend('', frameon=False)
plt.subplot(3, 2, 4)
plt.plot(time, rd4.TGb.DynEllip, color=vplot.colors.dark_blue, label="b")
plt.plot(time, rd4.TGc.DynEllip, color=vplot.colors.red, label="c")
plt.xlabel("Time (kyr)")
plt.ylabel("Dynamical Ellipticity")
plt.legend(loc="upper right", fontsize=12, ncol=1)
plt.subplot(3, 2, 5)
plt.plot(time, rd4.TGb.CassiniOne, color=vplot.colors.dark_blue)
plt.plot(time, rd4.TGc.CassiniOne, color=vplot.colors.red)
plt.xlabel("Time (kyr)")
plt.ylabel("$\sin{\Psi}$")
plt.legend('', frameon=False)
plt.subplot(3, 2, 6)
plt.plot(time, rd4.TGb.CassiniTwo, color=vplot.colors.dark_blue)
plt.plot(time, rd4.TGc.CassiniTwo, color=vplot.colors.red)
plt.xlabel("Time (kyr)")
plt.ylabel("$\cos{\Psi}$")
plt.legend('', frameon=False)
# Save the figure
ext = get_args().ext
fig.tight_layout()
fig.savefig(path / f"CassiniMulti4.{ext}")
#################### 2nd order DistOrb ####################
fig = plt.figure(figsize=(6.5, 8))
plt.subplot(3, 2, 1)
plt.plot(time, ll2.TGb.Obliquity, color=vplot.colors.dark_blue)
plt.plot(time, ll2.TGc.Obliquity, color=vplot.colors.red)
plt.xlabel("Time (kyr)")
plt.ylabel(r"Obliquity ($^\circ$)")
plt.legend('', frameon=False)
plt.subplot(3, 2, 2)
plt.plot(time, ll2.TGb.Eccentricity, color=vplot.colors.dark_blue)
plt.plot(time, ll2.TGc.Eccentricity, color=vplot.colors.red)
plt.xlabel("Time (kyr)")
plt.ylabel("Eccentricity")
plt.legend('', frameon=False)
plt.subplot(3, 2, 3)
plt.plot(time, ll2.TGb.RotPer, color=vplot.colors.dark_blue)
plt.plot(time, ll2.TGc.RotPer, color=vplot.colors.red)
plt.xlabel("Time (kyr)")
plt.ylabel("Rotation Period (days)")
plt.legend('', frameon=False)
plt.subplot(3, 2, 4)
plt.plot(time, ll2.TGb.DynEllip, color=vplot.colors.dark_blue, label="b")
plt.plot(time, ll2.TGc.DynEllip, color=vplot.colors.red, label="c")
plt.xlabel("Time (kyr)")
plt.ylabel("Dynamical Ellipticity")
plt.legend(loc="upper right", fontsize=12, ncol=1)
plt.subplot(3, 2, 5)
plt.plot(time, ll2.TGb.CassiniOne, color=vplot.colors.dark_blue)
plt.plot(time, ll2.TGc.CassiniOne, color=vplot.colors.red)
plt.xlabel("Time (kyr)")
plt.ylabel("$\sin{\Psi}$")
plt.legend('', frameon=False)
plt.subplot(3, 2, 6)
plt.plot(time, ll2.TGb.CassiniTwo, color=vplot.colors.dark_blue)
plt.plot(time, ll2.TGc.CassiniTwo, color=vplot.colors.red)
plt.xlabel("Time (kyr)")
plt.ylabel("$\cos{\Psi}$")
plt.legend('', frameon=False)
# Save the figure
ext = get_args().ext
fig.tight_layout()
fig.savefig(path / f"CassiniMulti2.{ext}")