-
Notifications
You must be signed in to change notification settings - Fork 55
/
Copy pathmakeplot.py
145 lines (108 loc) · 4.17 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
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
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
# Typical plot parameters that make for pretty plot
mpl.rcParams["figure.figsize"] = (10, 8)
mpl.rcParams["font.size"] = 15.0
# Run vplanet
lc17 = vplanet.run(path / "LehmerCatling17" / "vpl.in", units=False)
dyn = vplanet.run(path / "Dynamic" / "vpl.in", units=False)
### First figure: Comparative atmospheric escape ###
# Plot
fig, axes = plt.subplots(nrows=4, ncols=3, sharex=True, figsize=(16, 11))
time = lc17.star.Time
## Upper left: Envelope mass ##
axes[0, 0].plot(time, lc17.planet.EnvelopeMass, color="k")
axes[0, 0].plot(time, dyn.planet.EnvelopeMass, color=vplot.colors.dark_blue)
# Format
axes[0, 0].set_ylim(-0.01, 0.07)
axes[0, 0].set_ylabel(r"Env. Mass ($M_\oplus$)")
## Upper middle: Total Mass ##
axes[0, 1].plot(time, lc17.planet.Mass, color="k")
axes[0, 1].plot(time, dyn.planet.Mass, color=vplot.colors.dark_blue)
# Format
axes[0, 1].set_ylim(1.93, 2)
axes[0, 1].set_ylabel(r"Total Mass ($M_\oplus$)")
## Upper right: Envelope mass time derivative ##
axes[0, 2].plot(time, lc17.planet.DEnvMassDt / 1e9, color="k")
axes[0, 2].plot(time, dyn.planet.DEnvMassDt / 1e9, color=vplot.colors.dark_blue)
axes[0, 2].set_ylim(-6.5, 0.1)
# axes[0,2].set_yscale("symlog", linthreshy=0.1)
axes[0, 2].set_ylabel(r"$\dot{M}_{envelope}$ [$10^9$ kg/s]")
## Upper Middle left: Stellar Luminosity
axes[1, 0].plot(time, lc17.star.Luminosity, color="k")
axes[1, 0].plot(time, dyn.star.Luminosity, color=vplot.colors.dark_blue)
# Format
axes[1, 0].set_ylim(0.25, 2)
axes[1, 0].set_ylabel(r"Luminosity ($L_\odot$)")
## Upper Middle middle: XUV Luminosity
axes[1, 1].plot(time, lc17.star.LXUVTot * 1000, color="k")
axes[1, 1].plot(time, dyn.star.LXUVTot * 1000, color=vplot.colors.dark_blue)
# Format
axes[1, 1].set_ylim(0.1, 2)
axes[1, 1].set_ylabel(r"$L_{XUV}$ ($10^{-3}L_\odot$)")
## Upper Right middle: Incident XUV flux
axes[1, 2].plot(time, lc17.planet.FXUV, color="k", label="L+C (2017)")
axes[1, 2].plot(
time, dyn.planet.FXUV, color=vplot.colors.dark_blue, label="Auto AtmEsc"
)
axes[1, 2].legend(loc="upper right")
# Format
axes[1, 2].set_ylim(50, 250)
axes[1, 2].set_ylabel(r"$F_{XUV}$ ($W/m^2$)")
## Lower Middle left: Planetary radius
axes[2, 0].plot(time, lc17.planet.PlanetRadius, color="k")
axes[2, 0].plot(time, dyn.planet.PlanetRadius, color=vplot.colors.dark_blue)
# Format
axes[2, 0].set_ylim(0, 30)
axes[2, 0].set_ylabel(r"Radius ($R_\oplus$)")
## Lower Middle middle: Planetary XUV radius
axes[2, 1].plot(time, lc17.planet.RadXUV, color="k")
axes[2, 1].plot(time, dyn.planet.RadXUV, color=vplot.colors.dark_blue)
# Format
axes[2, 1].set_ylim(0, 30)
axes[2, 1].set_ylabel(r"XUV Radius ($R_\oplus$)")
## Lower Middle right: Roche radius
axes[2, 2].plot(time, lc17.planet.RocheRadius, color="k")
axes[2, 2].plot(time, dyn.planet.RocheRadius, color=vplot.colors.dark_blue)
# Format
axes[2, 2].set_ylim(29.1, 29.7)
axes[2, 2].set_ylabel(r"Roche Radius ($R_\oplus$)")
## Lower left: Scale height
axes[3, 0].plot(time, lc17.planet.ScaleHeight, color="k")
axes[3, 0].plot(time, dyn.planet.ScaleHeight, color=vplot.colors.dark_blue)
# Format
axes[3, 0].set_ylim(190, 460)
axes[3, 0].set_ylabel("Scale Height (km)")
## Lower middle: Surface Pressure
axes[3, 1].plot(time, lc17.planet.PresSurf, color="k")
axes[3, 1].plot(time, dyn.planet.PresSurf, color=vplot.colors.dark_blue)
# Format
axes[3, 1].set_ylim(-0.2, 7.9)
axes[3, 1].set_ylabel("Surf. Press. (GPa)")
## Lower right: Thermal Temperature
axes[3, 2].plot(time, lc17.planet.ThermTemp, color="k")
axes[3, 2].plot(time, dyn.planet.ThermTemp, color=vplot.colors.dark_blue)
# Format
axes[3, 2].set_ylim(780, 1500)
axes[3, 2].set_ylabel("Therm. Temp. (K)")
# Format
axes[3, 1].set_xlabel("Time [yr]")
# Format all axes
for ax in axes.flatten():
# Format x axis
ax.set_xlim(time[1], time.max())
ax.set_xscale("log")
# Set rasterization
ax.set_rasterization_zorder(0)
# Save the figure
ext = get_args().ext
fig.savefig(path / f"MiniNeptuneEvap.{ext}", bbox_inches="tight", dpi=600)