Skip to content

Commit 9086d64

Browse files
RoryBarnesRoryBarnes
authored
Merge pull request #345 from VirtualPlanetaryLaboratory/PoiseOutputFix
Poise output fix
2 parents 294a6a2 + 4107793 commit 9086d64

5 files changed

Lines changed: 351 additions & 267 deletions

File tree

‎examples/EarthClimate/makeplot.py‎

Lines changed: 19 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -91,22 +91,25 @@ def comp2huybers(plname, xrange=False, show=True):
9191
pco2 = float(lines[i].split()[1])
9292

9393
try:
94-
longp = (body.ArgP + body.LongA + body.PrecA + 180)
94+
longp = body.ArgP + body.LongA + body.PrecA + 180
9595
except:
9696
longp = body.PrecA
9797

9898
esinv = ecc * np.sin(longp)
9999

100-
lats = np.unique(body.Latitude)
100+
obl = np.array(obl)
101+
ecc = np.array(ecc)
102+
esinv = np.array(esinv)
103+
104+
lats = np.unique(body.Latitude.value)
101105
nlats = len(lats)
102106
ntimes = len(body.Time)
107+
time = body.Time.value
103108

104109
# plot temperature
105-
temp = np.reshape(body.TempLandLat, (ntimes, nlats))
110+
temp = np.reshape(body.TempLandLat.value, (ntimes, nlats))
106111
ax1 = plt.subplot(7, 1, 5)
107-
c = plt.contourf(
108-
body.Time / 1e6, lats[lats > 58 * u.deg], temp.T[lats > 58 * u.deg], 20
109-
)
112+
c = plt.contourf(time / 1e6, lats[lats > 58], temp.T[lats > 58], 20)
110113
plt.ylabel("Latitude")
111114

112115
plt.ylim(60, 83)
@@ -128,13 +131,11 @@ def comp2huybers(plname, xrange=False, show=True):
128131
clb.set_label("Surface Temp.\n($^{\circ}$C)", fontsize=12)
129132

130133
# plot ice height
131-
ice = np.reshape(body.IceHeight + body.BedrockH, (ntimes, nlats))
134+
ice = np.reshape((body.IceHeight + body.BedrockH).value, (ntimes, nlats))
132135

133136
ax3 = plt.subplot(7, 1, 4)
134137

135-
c = plt.contourf(
136-
body.Time / 1e6, lats[lats > 58 * u.deg], ice.T[lats > 58 * u.deg, :], 20
137-
)
138+
c = plt.contourf(time / 1e6, lats[lats > 58], ice.T[lats > 58, :], 20)
138139
plt.ylabel("Latitude")
139140
plt.ylim(60, 83)
140141

@@ -146,10 +147,8 @@ def comp2huybers(plname, xrange=False, show=True):
146147
clb.set_label("Ice sheet\nheight (m)", fontsize=12)
147148

148149
ax4 = plt.subplot(7, 1, 6)
149-
acc = np.reshape(body.IceAccum, (ntimes, nlats))
150-
c = plt.contourf(
151-
body.Time / 1e6, lats[lats > 58 * u.deg], acc.T[lats > 58 * u.deg], 20
152-
)
150+
acc = np.reshape(body.IceAccum.value, (ntimes, nlats))
151+
c = plt.contourf(time / 1e6, lats[lats > 58], acc.T[lats > 58], 20)
153152
plt.ylabel("Latitude")
154153
plt.ylim(60, 83)
155154
plt.yticks([60, 70, 80])
@@ -171,10 +170,8 @@ def comp2huybers(plname, xrange=False, show=True):
171170
clb.set_label("Ice Accum.\n(m year$^{-1}$)", fontsize=12)
172171

173172
ax5 = plt.subplot(7, 1, 7)
174-
abl = np.reshape(body.IceAblate, (ntimes, nlats))
175-
c = plt.contourf(
176-
body.Time / 1e6, lats[lats > 58 * u.deg], -abl.T[lats > 58 * u.deg], 20
177-
)
173+
abl = np.reshape(body.IceAblate.value, (ntimes, nlats))
174+
c = plt.contourf(time / 1e6, lats[lats > 58], -abl.T[lats > 58], 20)
178175
plt.ylabel("Latitude")
179176
plt.ylim(60, 83)
180177
plt.yticks([60, 70, 80])
@@ -194,7 +191,7 @@ def comp2huybers(plname, xrange=False, show=True):
194191

195192
plt.subplot(7, 1, 1)
196193
plt.plot(
197-
body.Time / 1e6,
194+
time / 1e6,
198195
obl,
199196
linestyle="solid",
200197
marker="None",
@@ -213,7 +210,7 @@ def comp2huybers(plname, xrange=False, show=True):
213210

214211
plt.subplot(7, 1, 2)
215212
plt.plot(
216-
body.Time / 1e6,
213+
time / 1e6,
217214
ecc,
218215
linestyle="solid",
219216
marker="None",
@@ -232,7 +229,7 @@ def comp2huybers(plname, xrange=False, show=True):
232229

233230
plt.subplot(7, 1, 3)
234231
plt.plot(
235-
body.Time / 1e6,
232+
time / 1e6,
236233
esinv,
237234
linestyle="solid",
238235
marker="None",
@@ -329,7 +326,7 @@ def seasonal_maps(time, show=True):
329326
if p == len(output.bodies) - 1 and ctmp == 0:
330327
raise Exception("Planet %s not found." % plname)
331328

332-
lats = np.unique(body.Latitude)
329+
lats = np.unique(body.Latitude.value)
333330
try:
334331
obl = body.Obliquity[np.where(body.Time == time)[0]]
335332
except:

‎examples/EarthClimate/sun.in‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,6 @@ dMass 1
44
dSemi 0
55
dEcc 0
66
dRadius 0.00135
7-
dLuminosity 3.846e26
7+
dLuminosity -1
88
sStellarModel none
99
saModules stellar

‎examples/IceBelts/makeplot.py‎

Lines changed: 13 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -92,15 +92,16 @@ def clim_evol(plname, xrange=False, show=True):
9292

9393
esinv = ecc * np.sin(longp) * np.sin(obl * np.pi / 180.0)
9494

95-
lats = np.unique(body.Latitude)
95+
lats = np.unique(body.Latitude.value)
9696
nlats = len(lats)
9797
ntimes = len(body.Time)
98+
time = body.Time.value
9899

99100
# plot temperature
100-
temp = np.reshape(body.TempLat, (ntimes, nlats))
101+
temp = np.reshape(body.TempLat.value, (ntimes, nlats))
101102
ax1 = plt.subplot(5, 1, 1)
102103

103-
c = plt.contourf(body.Time, lats, temp.T, cmap="plasma")
104+
c = plt.contourf(time, lats, temp.T, cmap="plasma")
104105
plt.ylabel(r"Latitude [$^\circ$]", fontsize=10)
105106
plt.title(r"Surface Temp [$^{\circ}$C]", fontsize=12)
106107
plt.ylim(-90, 90)
@@ -116,9 +117,9 @@ def clim_evol(plname, xrange=False, show=True):
116117
plt.setp(cbar.ax.yaxis.get_ticklabels(), fontsize=9)
117118

118119
# plot albedo
119-
alb = np.reshape(body.AlbedoLat, (ntimes, nlats))
120+
alb = np.reshape(body.AlbedoLat.value, (ntimes, nlats))
120121
ax2 = plt.subplot(5, 1, 3)
121-
c = plt.contourf(body.Time, lats, alb.T, cmap="Blues_r")
122+
c = plt.contourf(time, lats, alb.T, cmap="Blues_r")
122123
plt.ylabel(r"Latitude [$^\circ$]", fontsize=10)
123124
plt.title("Albedo [TOA]", fontsize=12)
124125
plt.ylim(-90, 90)
@@ -130,9 +131,9 @@ def clim_evol(plname, xrange=False, show=True):
130131
plt.setp(cbar.ax.yaxis.get_ticklabels(), fontsize=9)
131132

132133
# plot ice height
133-
ice = np.reshape(body.IceHeight, (ntimes, nlats))
134+
ice = np.reshape(body.IceHeight.value, (ntimes, nlats))
134135
ax3 = plt.subplot(5, 1, 4)
135-
c = plt.contourf(body.Time, lats, ice.T, cmap="Blues_r")
136+
c = plt.contourf(time, lats, ice.T, cmap="Blues_r")
136137
plt.ylabel(r"Latitude [$^\circ$]", fontsize=10)
137138
plt.title("Ice sheet height [m]", fontsize=12)
138139
plt.ylim(-90, 90)
@@ -144,9 +145,9 @@ def clim_evol(plname, xrange=False, show=True):
144145
plt.setp(cbar.ax.yaxis.get_ticklabels(), fontsize=9)
145146

146147
# plot bedrock
147-
brock = np.reshape(body.BedrockH, (ntimes, nlats))
148+
brock = np.reshape(body.BedrockH.value, (ntimes, nlats))
148149
ax4 = plt.subplot(5, 1, 5)
149-
c = plt.contourf(body.Time, lats, brock.T, cmap="Reds_r")
150+
c = plt.contourf(time, lats, brock.T, cmap="Reds_r")
150151
plt.ylabel(r"Latitude [$^\circ$]", fontsize=10)
151152
plt.title("Bedrock height [m]", fontsize=12)
152153
plt.ylim(-90, 90)
@@ -159,9 +160,9 @@ def clim_evol(plname, xrange=False, show=True):
159160
plt.setp(cbar.ax.yaxis.get_ticklabels(), fontsize=9)
160161

161162
# plot insolation
162-
insol = np.reshape(body.AnnInsol, (ntimes, nlats))
163+
insol = np.reshape(body.AnnInsol.value, (ntimes, nlats))
163164
ax5 = plt.subplot(5, 1, 2)
164-
c = plt.contourf(body.Time, lats, insol.T, cmap="plasma")
165+
c = plt.contourf(time, lats, insol.T, cmap="plasma")
165166
plt.ylabel(r"Latitude [$^\circ$]", fontsize=10)
166167
plt.title(r"Annual average insolation [W/m$^2$]", fontsize=12)
167168
plt.ylim(-90, 90)
@@ -201,7 +202,6 @@ def seasonal_maps(time, show=True):
201202

202203
check = 0
203204
for f in glob.glob(str(path / "SeasonalClimateFiles" / "*.DailyInsol.*")):
204-
205205
f1 = f.split(".")
206206

207207
if len(f1) == 4:
@@ -249,7 +249,7 @@ def seasonal_maps(time, show=True):
249249
if p == len(output.bodies) - 1 and ctmp == 0:
250250
raise Exception("Planet %s not found" % plname)
251251

252-
lats = np.unique(body.Latitude)
252+
lats = np.unique(body.Latitude.value)
253253
try:
254254
obl = body.Obliquity[np.where(body.Time == time)[0]]
255255
except:

‎examples/StellarEvol/makeplot.py‎

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,9 @@
1818
# Colormap
1919
cmap = plt.get_cmap("inferno")
2020

21+
# Rearth per Rsun (output is in Rearth, plot in Rsun)
22+
dRearthPerRsun = 6.957e8 / 6.3781e6
23+
2124
# Star input file template
2225
star = """#
2326
sName s%02d
@@ -90,7 +93,9 @@ def run(masses):
9093
# Plot stars that survived
9194
for n, m in enumerate(masses):
9295
# Top row: radius, legend
93-
ax[0, 0].plot(time, radius[n], label="%.2f" % m, color=cmap(0.7 * m))
96+
ax[0, 0].plot(
97+
time, radius[n] / dRearthPerRsun, label="%.2f" % m, color=cmap(0.7 * m)
98+
)
9499
# Dummy data for legend
95100
ax[0, 1].plot([101], [100], label="%.2f" % m, color=cmap(0.7 * m))
96101

@@ -109,7 +114,7 @@ def run(masses):
109114
# Top row: radius, legend
110115
m = 1.3
111116
time = dead[:, 0]
112-
ax[0, 0].plot(time, dead[:, 4], label="%.2f" % m, color=cmap(0.7 * m))
117+
ax[0, 0].plot(time, dead[:, 4] / dRearthPerRsun, label="%.2f" % m, color=cmap(0.7 * m))
113118

114119
# Dummy data for legend
115120
ax[0, 1].plot([101], [100], label="%.2f" % m, color=cmap(0.7 * m))

0 commit comments

Comments
 (0)