Skip to content

Poise output fix #345

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
RoryBarnes merged 3 commits into from
Mar 4, 2026
+351 −267
Merged

Poise output fix #345

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2524154
Fixed typo in 4 POISE output options.
RoryBarnes Mar 4, 2026
94dabde
Fixed issue in which some examples were failing with a new astropy de…
RoryBarnes Mar 4, 2026
4107793
Fixed dLuminosity in example/EarthClimate/sun.in
RoryBarnes Mar 4, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 19 additions & 22 deletions examples/EarthClimate/makeplot.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,22 +91,25 @@ def comp2huybers(plname, xrange=False, show=True):
pco2 = float(lines[i].split()[1])

try:
longp = (body.ArgP + body.LongA + body.PrecA + 180)
longp = body.ArgP + body.LongA + body.PrecA + 180
except:
longp = body.PrecA

esinv = ecc * np.sin(longp)

lats = np.unique(body.Latitude)
obl = np.array(obl)
ecc = np.array(ecc)
esinv = np.array(esinv)

lats = np.unique(body.Latitude.value)
nlats = len(lats)
ntimes = len(body.Time)
time = body.Time.value

# plot temperature
temp = np.reshape(body.TempLandLat, (ntimes, nlats))
temp = np.reshape(body.TempLandLat.value, (ntimes, nlats))
ax1 = plt.subplot(7, 1, 5)
c = plt.contourf(
body.Time / 1e6, lats[lats > 58 * u.deg], temp.T[lats > 58 * u.deg], 20
)
c = plt.contourf(time / 1e6, lats[lats > 58], temp.T[lats > 58], 20)
plt.ylabel("Latitude")

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

# plot ice height
ice = np.reshape(body.IceHeight + body.BedrockH, (ntimes, nlats))
ice = np.reshape((body.IceHeight + body.BedrockH).value, (ntimes, nlats))

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

c = plt.contourf(
body.Time / 1e6, lats[lats > 58 * u.deg], ice.T[lats > 58 * u.deg, :], 20
)
c = plt.contourf(time / 1e6, lats[lats > 58], ice.T[lats > 58, :], 20)
plt.ylabel("Latitude")
plt.ylim(60, 83)

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

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

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

plt.subplot(7, 1, 1)
plt.plot(
body.Time / 1e6,
time / 1e6,
obl,
linestyle="solid",
marker="None",
Expand All @@ -213,7 +210,7 @@ def comp2huybers(plname, xrange=False, show=True):

plt.subplot(7, 1, 2)
plt.plot(
body.Time / 1e6,
time / 1e6,
ecc,
linestyle="solid",
marker="None",
Expand All @@ -232,7 +229,7 @@ def comp2huybers(plname, xrange=False, show=True):

plt.subplot(7, 1, 3)
plt.plot(
body.Time / 1e6,
time / 1e6,
esinv,
linestyle="solid",
marker="None",
Expand Down Expand Up @@ -329,7 +326,7 @@ def seasonal_maps(time, show=True):
if p == len(output.bodies) - 1 and ctmp == 0:
raise Exception("Planet %s not found." % plname)

lats = np.unique(body.Latitude)
lats = np.unique(body.Latitude.value)
try:
obl = body.Obliquity[np.where(body.Time == time)[0]]
except:
Expand Down
2 changes: 1 addition & 1 deletion examples/EarthClimate/sun.in
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,6 @@ dMass 1
dSemi 0
dEcc 0
dRadius 0.00135
dLuminosity 3.846e26
dLuminosity -1
sStellarModel none
saModules stellar
26 changes: 13 additions & 13 deletions examples/IceBelts/makeplot.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,15 +92,16 @@ def clim_evol(plname, xrange=False, show=True):

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

lats = np.unique(body.Latitude)
lats = np.unique(body.Latitude.value)
nlats = len(lats)
ntimes = len(body.Time)
time = body.Time.value

# plot temperature
temp = np.reshape(body.TempLat, (ntimes, nlats))
temp = np.reshape(body.TempLat.value, (ntimes, nlats))
ax1 = plt.subplot(5, 1, 1)

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

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

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

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

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

check = 0
for f in glob.glob(str(path / "SeasonalClimateFiles" / "*.DailyInsol.*")):

f1 = f.split(".")

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

lats = np.unique(body.Latitude)
lats = np.unique(body.Latitude.value)
try:
obl = body.Obliquity[np.where(body.Time == time)[0]]
except:
Expand Down
9 changes: 7 additions & 2 deletions examples/StellarEvol/makeplot.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,9 @@
# Colormap
cmap = plt.get_cmap("inferno")

# Rearth per Rsun (output is in Rearth, plot in Rsun)
dRearthPerRsun = 6.957e8 / 6.3781e6

# Star input file template
star = """#
sName s%02d
Expand Down Expand Up @@ -90,7 +93,9 @@ def run(masses):
# Plot stars that survived
for n, m in enumerate(masses):
# Top row: radius, legend
ax[0, 0].plot(time, radius[n], label="%.2f" % m, color=cmap(0.7 * m))
ax[0, 0].plot(
time, radius[n] / dRearthPerRsun, label="%.2f" % m, color=cmap(0.7 * m)
)
# Dummy data for legend
ax[0, 1].plot([101], [100], label="%.2f" % m, color=cmap(0.7 * m))

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

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