Eccentricity fixes (#1)

[wizzeh]

Hi!

I wanted to do some simulations on synchronously orbiting worlds at higher eccentricities. I increased ORB_ECCEN_MAX and added functionality so the subsolar point would move.

[wizzeh]

Uh I wouldn't accept without checking the signs though lol

[iagocg]

Yes! I have the same problem, I want to set eccentricity = 0.2 and it crashes. How can I use this fix? Will it be implemented in future versions of ExoPlaSim?

[alphaparrot]

Hi folks, let me look into how physics-breaking these changes would be. Keep in mind that eccentricity > 0.2 starts to get very extreme for a habitable planet; you're talking about variations in incident flux that span much of the width of the habitable zone. Depending on length of model year, rotation rate, etc, that can be model-breaking. Most immediately, that limit is there for orbital calculations related to Earth and paleo-Earth that will probably break at higher eccentricities, so setting eccentricity above 0.2 will likely have to be paired with forcing a fixed orbit. If it all checks out, I'll add this to the next update.

[iagocg]

Thank you very much. I've made my calculations by hand and with an eccentricity of 0.2 (the most I dared) and with periapsis at the southern solstice, the northern "summer" is 62.6% of the year and the northern "winter" is 37.4% of the year (understanding summer and winter from equinox to equinox), which is quite strong, but still not too hostile to life. Anyway, these calculations are for written fiction, so I'd would love to see higher eccentricities implemented even if you set a warning explaining that the results may be compromised and that we should interpret them with a pinch of salt. Thank you so much for such an extremely useful and fun program.

[robmod]

I too am interested in doing high-E simulations.I've been struggling to get exoplasim running on my macbook, but if it can't do high-E sims there wouldn't be much point in my continuing my effort.

Is high eccentricity possible currently? Can it be done be editing the code in some way?

I'm trying to planets with rotations rates of 6-48 hrs, around G-F type stars, in orbits of at least 1 AU. Eccentricity should be between 0.1 - 0.4. @wizzeh suggest this can be done through manually editing ORB_ECCEN_MAX?

@alphaparrot Can I do this? If so, any tips? You mention forcing a "fixed orbit." Not sure what you mean there.

I'm okay with "hackey" solutions. I'm not looking for research-grade results, just something that is accurate-enough.

[robmod]

If it matters, the high-E models I want to run would cap out at about 2-4,000 w/m2 at perihelion (mostly 2500 w/m2 or lower). So 2-4x Earth's solar isolation maximum.

[wizzeh]

Hi @robmod,

You should be able to do that just by manually editing ORB_ECCEN_MAX as I do in this PR. I looked into it a little and everything should be robust at 0.4 eccentricity, especially if you are just looking for a hacky result.

There is one exception which is that (please correct me if I'm wrong but) ExoPlaSim can optionally simulate Milankovich cycles and this is the calculation which would stop making sense if the eccentricity becomes non-earthlike. As long as you set fixedorbit=True in your model configuration these calculations won't happen and you shouldn't have to worry about it.

[robmod]

Hi @wizzeh.

That's great! While I'm ok with "hackey" I would like as realistic a result as feasible. Do you foresee any particular things that might be a little weird with the model at eccentricities up to 0.4? It'd be good to know what to take with a grain-of-salt.

I'm totally fine with not using the Milankovich cycles. I'm actually using the planet/orbit data from n-body simulations I've run of full exoplanet systems, so the cyclees from ExoPlaSim just wouldn't correspond to the far more complicated cycles I see in the n-body simulations.

[wizzeh]

I can't speak to the climate model but nothing should go wrong with the orbital calculations.

[hersfeldtn]

If it matters, the high-E models I want to run would cap out at about 2-4,000 w/m2 at perihelion (mostly 2500 w/m2 or lower). So 2-4x Earth's solar isolation maximum.

That high insolation might push the temperature high enough to crash the model even if the eccentricity calculation works fine, so you may have to watch out for that.

[hersfeldtn]

Also while I'm here @wizzeh do you know how the true anomaly calculation works? Wondering if it might be possible to put in some newton-raphson iteration to ensure accuracy up to high eccentricity.

[robmod]

@wizzeh Darn. I tried changing ORB_ECCEN_MAX to 0.4. I then changed the eccentricity of my test model earth to 0.11. The whole thing crashed almost immediately. Didn't even run the first year. Got the below. Any idea what I did wrong? I didn't make those other changes to code, but I presumed those were there for tidal-locking scenarios?

You've run high-E (i.e. up to 0.4) models already with your changes and can confirm they work?

Command '['mpiexec -np 4 most_plasim_t21_l10_p4.x']' returned non-zero exit status 231.
Traceback (most recent call last):
  File "/Users/rob/Library/Python/3.9/lib/python/site-packages/exoplasim/__init__.py", line 521, in runtobalance
    subprocess.run([self._exec+self.executable],shell=True,check=True,
  File "/Applications/Xcode.app/Contents/Developer/Library/Frameworks/Python3.framework/Versions/3.9/lib/python3.9/subprocess.py", line 528, in run
    raise CalledProcessError(retcode, process.args,
subprocess.CalledProcessError: Command '['mpiexec -np 4 most_plasim_t21_l10_p4.x']' returned non-zero exit status 231.

Update: I don't think it has anything to do with temperature, since I just tried a version with a perihelion distance that is a bit further out than 1 AU—so a cold world. Same error.

[robmod]

@wizzeh Huh. Out of curiosity i'm looking at how the eccentricity variable is used in radmod.f90 and...well...its really unusual to my eyes. I'm familiar with n-body orbital calculations, as well as general keplerian calculations, and I'm not seeing anything I recognize. I'm kind of baffled how eccentricity is utilized here (doesn't help I'm unfamiliar with this kind of code).

Would you be cool with explaining your thinking on the changes you've made to eccentricity in this code? I had presumed eccentricity is used to calculate flux at a given time (which is usually a straightforward though long series of calculations), but I get the sense other things are going on...

[robmod]

@hersfeldtn Have you seen in your runs what a max insolation is before the model breaks? Is the breaking insolation dependent or temperature dependent? I have presumed that the model can handle high insolation for a very short time since the temperature would lag the high flux and not reach a crazy value, while it would break at a constant high insolation. But that's a guess.

Also, could you elaborate on your true-anomaly comment? What's the concern?

[alphaparrot]

Hi folks; apologies for the delay on this--having a full time job sadly eats up a lot of time that could go to model development. I'll try to set aside some time for an eccentricity update over the next week. To answer some of the questions that have been raised the last few days:

The fundamental way that ExoPlaSim does its atmospheric dynamics is by transforming the physical variables (wind speed, humidity, temperature, vorticity, etc) into frequency space using a spherical version of a Fourier transform--essentially, decomposing the activity in the atmosphere into spherical modes of different wavelengths. Solving the partial differential equations that govern fluid transport is very easy and efficient in frequency space, which is why the model does this (and is also a big part of why it's so fast at low/moderate resolutions compared to other models--though this method scales as N^3, so it's bad at very high resolution). The results then get transformed back into physical space for physics calculations like radiation, moist physics, etc. If there's extreme stuff going on, like high winds, intense heating, intense cooling, very high temperatures, very low temperatures, etc, that can produce unphysical negative values in the spectral math being done, which can cause crashes later (e.g. a square root or logarithm of a negative number is not good). Small deviations, such as are common in the water vapor field, can be filled by the dynamical solver in order to satisfy rules like mass and energy conservation, but larger excursions from physical reasonableness can't be fixed on the fly. Note that non-spectral models also have these same limitations, but they usually manifest as violations of either energy/mass conservation or violations of the Courant-Friedrich-Levy Condition (material can't travel more than a grid cell in a timestep and still have a unique solution for the fluid transport equations).

Often a shorter timestep will break those extreme changes into less-extreme changes, and allow the solver to simulate that event without too many problems, but not always--this is why at the limits of the model's comfort zone, crashes are more likely and you have to start fiddling with e.g. changing the pressure or insolation by a tiny amount (to get out of whatever local minimum was causing a local problem), and going to progressively shorter and shorter timesteps. There is a point however at which you are just at or beyond the limit of what the model can handle mathematically. With high eccentricity, due to the extreme temperature swings and high solar forcing you're likely to encounter, this sort of thing is going to be more likely--there's just no way around that. A lower-dimensional/lower-complexity model, such as climlab, might be able to run without breaking on problems like that, and give you things like temperature, ice fraction, humidity, and so on, at the cost of granular detail and physical consistency/complexity. Once you're in 3D with lots of physics, there are sadly some problems which require behemoth-scale codes and very large computers in order to simulate without crashing. ExoPlaSim is more robust in extreme parameter spaces than some of its slower cousins, but of course that comes at the cost of physical reliability and accuracy, and it has its own limits.

@hersfeldtn true anomaly isn't really used at the moment in ExoPlaSim; solar declination and current distance from the planet is computed with an equation/parameterization that I don't fully trust at high eccentricity--thus why for something like this I'd probably replace it with a better orbital parameters code. If you check out the 'westeros' branch of the code, you'll actually find some useful code (that I'll likely bring forward into the main code), which does use a Newton-Raphson iterator to compute true anomaly and stellar position in the sky and such. That branch was for an experimental version of the model used to simulate the configuration in the April Fools' paper I wrote a few years ago, but the orbital math is fairly general--the idea was that eventually it would be useful for circumbinary systems, systems hooked up to N-body integrators, etc.

@robmod, @wizzeh is correct -- the fixed orbit thing refers to whether or not Milkankovich cycles and orbital precession are computed. ExoPlaSim has code for doing that for Earth; for exoplanets and fictional planets I generally recommend turning it off with fixedorbit=True. If you've modified the code already to allow high eccentricity, it may just be a matter of playing with the timestep to get it to run without crashing. Check that a standard Earth sim still runs though to make sure you haven't just broken the model.

The way eccentricity is currently used is with the orb_decl routine in radmod.f90, which computes a distance factor (1/r^2) at each timestep in the form of the variable eccf, which is then passed to the radiation calculation as gdist2 to modulate the insolation at the semimajor axis (which is what you prescribe when you configure the model).

Something to consider: there will be a difference between high-eccentricity around higher-mass (F) stars, and high-eccentricity around lower-mass (M and K) stars. A longer orbital period means you spend longer at periapsis and apoapsis, resulting in larger temperature extremes--but smoother changes in insolation. I don't know which would be more numerically robust; you guys are out in relatively uncharted waters, in terms of GCM modelling (sort of thing grad students do to write papers). If you ask "what would a planet with a particular land configuration look like around an F star with a semimajor axis insolation of 0.9 Earth insolation and an eccentricity of 0.35 and a rotation rate of 1.5 days," that's not a question anyone has asked before and it's likely you'll discover new climate regimes by running that experiment.

One of the major caveats to the physical reliability of ExoPlaSim when it comes to modelling out in these uncharted waters: ExoPlaSim does not have dynamic oceans nor dynamic sea and land ice. Glaciers don't flow. With large global changes in forcings, we would expect those things to become even more important. If you're using this for fun and fiction then maybe that doesn't matter, but it also means you have considerable creative license to imagine what the addition of ocean currents and glacial flow would mean for your model results.

[alphaparrot]

I should note--for those of you actively modifying the code and trying new things (I applaud you; this is fantastic); the debug=True flag can be useful (or you can adjust the compilation flags directly in most_compiler_mpi to have -Og -ggdb instead of -O3). You'll lose some speed, but maybe get more verbose output from the Fortran when it crashes (and run the .x directly with mpiexec to avoid all of the Python error handling), and that may help you to understand where the model is crashing--if it's a temperature problem, a humidity problem, etc.

[robmod]

@alphaparrot Great to hear from you! Super looking forward to that eccentricity update. :) The orbital calculation work you did for the Westeros case sounds a lot more sensible to my ears, at least in the case of exoplanets.

FYI – Unrelated to eccentricity, I've sent you an email concerning my install not outputting .nc files. It crashed when any output file type is specified that isn't the default. Don't know why. Hoping you can help at your convenience.

[wizzeh]

Also while I'm here @wizzeh do you know how the true anomaly calculation works? Wondering if it might be possible to put in some newton-raphson iteration to ensure accuracy up to high eccentricity.

@hersfeldtn I believe alphaparrot pretty much answered all the other questions but if you look at line 2911 of radmod.f90 you can see the true anomaly is calculated using a fourier expansion https://en.wikipedia.org/wiki/True_anomaly#From_the_mean_anomaly. You could add more terms if you needed to. In this branch I do sort of use true anomaly to get the subsolar point to move around the planet but ideally a better orbital model would be used rather than this hack.

[robmod]

@wizzeh Off-topic, but do you happen to know how to read the .npy file outputs? I can't figure out how to do it with matplotlib (or anything else). I only ask b/c you seem knowledgeable about ExoPlaSim.

[wizzeh]

@robmod .npy are binary numpy array files. You should be able to np.load them.

[robmod]

@wizzeh Oh, I can load them just fine. But I can't interpret them in a way that makes sense. I don't know what code to use. I tried some variations based of the example code, but those are for .nc files. Everything I try either output gobbleygook graphs or (more often) just produces errors.

[wizzeh]

Oh, I see. No, sorry, I've only really used the hdf5 and nc outputs personally.

[lynnbwilsoniii]

I am genuinely finding this all very fascinating but could someone perhaps remind me how I became subscribed to this discussion?

[wizzeh]

You should be able to unsubscribe but yeah if you have questions that aren't related to the PR feel free to email me contact@elliehigh.com

[alphaparrot]

Hey folks, brief update: I've pushed a sequence of commits (3e9dad5, be66051, and 3b6899e) that might end up satisfying the asks in this PR. Specifically, I've added in a routine for calculating the true anomaly for a generic Keplerian orbit, given the standard parameters and additionally the initial mean anomaly (important because e.g. Jan 1 is not a mean anomaly of 0; it's actually like -3 degrees). It uses a Newton-Raphson iterator to compute the eccentric anomaly and computes the true anomaly from that, so it should be reasonably accurate (to within machine precision) for just about any orbit you want to give it. That does not mean the model will run at any orbit. The usual climate considerations still apply; just because you can now do an orbit of 0.999 eccentricity does not mean that will produce a stable climate (though I'd be fascinated if it did). The only actual orbital limitation now is that the orbit must be bound (eccentricity<1). This mode is enabled with the keplerian=True flag when you configure the model, and setting the meananomaly0 parameter (in degrees). You can also do it directly in planet_namelist with NGENKEPLERIAN = 1 and MEANANOM0 = [whatever]. I haven't had time to test this on standard cases, but the code is committed if anyone wants to risk a try (but until I've had a chance to test, the current state of the github repository is not stable. Tests should come within the next week, hopefully, followed by a stable version release installable with pip.

If all is well with these commits, they will obviate the code in this PR (since that would still run up against the accuracy problems of the existing orbital math at high eccentricity).

These commits will also help with the longer development roadmap, which includes multiple light sources in the sky, and eventually the ability to easily couple to N-body integrators and/or arbitrary sequences of celestial positions.

[hersfeldtn]

Just to be clear, does meananomaly0 have any impact on the relative timing of apoapsis/periapsis and the solstices, or is that still totally controlled by lonvernaleq? (in which case does it basically just control what point of the year the output files start at?)

[alphaparrot]

Yes, meananomaly0 will affect timing of periapse--the mean anomaly is measured from periapse, so if you start at 0, you start at periapse. I believe I have the Earth default set correctly such that periapse is approximately January 4, and vernal equinox is still approximately March 20/21. Lonvernaleq is still relevant; just think of it as setting the spacing between the two. The relative math as far as defaults go (to ensure backwards compatibility) is one of the things I'll be testing this week, to make sure I got it right.

On Mon, Apr 17, 2023, 3:31 AM hersfeldtn @.***> wrote:

Just to be clear, does meananomaly0 have any impact on the relative timing of apoapsis/periapsis and the solstices, or is that still totally controlled by lonvernaleq? (in which case does it basically just control what point of the year the output files start at?)

— Reply to this email directly, view it on GitHub https://github.com/alphaparrot/ExoPlaSim/pull/2#issuecomment-1510846628, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADRQKQBBITH3SSB6ERIZ4D3XBTWWZANCNFSM6AAAAAAQVO43HE . You are receiving this because you were mentioned.Message ID: @.***>

[robmod]

Is anyone else experiencing a bug in which setting rotation rate to almost anything other than "1" changes the timing of the seasons (along with a corresponding change to flux strength/angle etc.)? Most notable is a total inversion of seasons when rotation rate is set to 1.5. (Note this started before the current eccentricity changes.) ...Trying to discern if its the code or something about the setup on my machine.

[robmod]

I'm getting a model crash when I use the updates. The way I implemented the updates was to replace the six modified files in my EPS install with the updated github versions. I then deleted "firstrun" and ran a model, forcing a recompile. I've used that method before without error, so if its a bad method I'd definitely like to know.

The errors I get are:

radmod.f90:857:39:

  857 |      &               ,mypid, nroot,nud)
      |                                       1
Error: Type mismatch in argument 'lambm0' at (1); passed LOGICAL(4) to REAL(8)
radmod.f90:857:39:

  857 |      &               ,mypid, nroot,nud)
      |                                       1
Error: Type mismatch in argument 'mvelpp' at (1); passed INTEGER(4) to REAL(8)
radmod.f90:857:39:

  857 |      &               ,mypid, nroot,nud)
      |                                       1
Error: Type mismatch in argument 'log_print' at (1); passed INTEGER(4) to LOGICAL(4)
make: *** [makefile:23: radmod.o] Error 1

It soon crashes with the following:

Command '['mpiexec -np 4 most_plasim_t21_l10_p4.x']' returned non-zero exit status 134.
Traceback (most recent call last):
  File "/Users/username/Library/Python/3.9/lib/python/site-packages/exoplasim/__init__.py", line 521, in runtobalance
    subprocess.run([self._exec+self.executable],shell=True,check=True,
  File "/Applications/Xcode.app/Contents/Developer/Library/Frameworks/Python3.framework/Versions/3.9/lib/python3.9/subprocess.py", line 528, in run
    raise CalledProcessError(retcode, process.args,
subprocess.CalledProcessError: Command '['mpiexec -np 4 most_plasim_t21_l10_p4.x']' returned non-zero exit status 134.

[robmod]

Could anyone here explain how Lonvernaleq works in relationship to eccentricity? I'm having a real hard time understanding it and making the proper settings. I want to be able to test 4 setups:

1. Winter Solstice coincides with perihelion
2. Spring Equinox coincides with perihelion
3. Summer Solstice coincides with perihelion
4. Autumn Equinox coincides with perihelion

For example I'll set Lonvernaleq so that setup 1 (Winter Solstice at perihelion) is the case (through trial-and-error), but if I add 180 to that Lonvernaleq value it doesn't do I what I hoped for—I don't get the opposite (Summer Solstice at perihelion).

[alphaparrot]

Hey all, brief update--fixed the bugged code with those function arguments. Unfortunately it looks like the orbital math might be bugged; my results for a 0-obliquity e=0.3 aquaplanet and Earth-obliquity e=0.3 aquaplanet don't look right. Will look more into it later, as well as check Rob's query for possible bugs in the code vs just needing to clarify terms. Hold on to the previous stable version before trying to get results with the current code.

On Sat, Apr 22, 2023 at 7:49 PM robmod @.***> wrote:

Could anyone here explain how Lonvernaleq works in relationship to eccentricity? I'm having a real hard time understanding it and making the proper settings. I want to be able to test 4 setups:

1. Winter Solstice coincides with perihelion
2. Spring Equinox coincides with perihelion
3. Summer Solstice coincides with perihelion
4. Autumn Equinox coincides with perihelion

For example I'll set Lonvernaleq so that #1 https://github.com/alphaparrot/ExoPlaSim/pull/1 (Winter Solstice at perihelion) is the case (through trial-and-error), but if I add 180 to that Lonvernaleq value it doesn't do I what I hoped for—I don't get the opposite (Summer Solstice at perihelion).

— Reply to this email directly, view it on GitHub https://github.com/alphaparrot/ExoPlaSim/pull/2#issuecomment-1518898121, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADRQKQGTWIZGQ2HXCPS2ZXDXCRU7RANCNFSM6AAAAAAQVO43HE . You are receiving this because you were mentioned.Message ID: @.***>

[robmod]

Darn. Well I for one appreciate the work you're putting in to get this code running! In the meantime, for those of us who aren't waiting on the new code to try high-E models, do you have an idea in what ways the native Plasim code used in EPS will prove inaccurate when doing high E models? Or how inaccurate?

Regarding my earlier issue (with the code before this eccentricity update)...I think I figured it out. Lonvernaleq is accurately setting angle between perihelion and autumn solstice (documentation says spring solstice...but that does not look correct). However changing rotation rate often changes at what time the simulation begins, and therefore changes the order of the output data. I had assumed that the beginning month was held constant, and so when I examined the data it looked like the seasons were getting messed up.

This quirk doesn't seem to be effecting the accuracy of the modeling, but it does make an exact perfect apple-to-apple comparison of data between some models with different rotation values more difficult (or not possible) since the snapshots/averaged months aren't representing the exact same time periods. For example, one model may (for sake of argument) begin the sim on January 1st, while another may begin it on August 13th...so the 12 averaged months in each model represent different time periods offset by 13 days.

Setting eccentricity with the old code does not seem to mess with the above. Zero Lonvernaleq is still autumn solstice=perihelion at any eccentricity.

[alphaparrot]

It's been a while since I was eyeballs-deep in the calendaring part of the code; that may be a consequence of calendar dates being ill-defined when the length of the day changes--the current code might use 24-hour days to determine dates (basically marking time by Earth sols for consistency). I'll look into potentially adding some of the orbital parameters as time-vector output codes to make apples-to-apples comparisons easier. It'll make debugging easier, that's for sure.

There are two major downsides to the older orbital code for higher eccentricities. One is that the estimate of solar declination and distance from the planet will become increasingly inaccurate at higher eccentricities, as the older code uses an approximation to solve for orbital parameters that only holds for lower eccentricities. The other is that (I think) it'll get the hour angle wrong for more eccentric orbits. The length of the solar day ought to vary with true anomaly given an eccentric orbit and a fixed sidereal day. The newer orbital code should handle that natively; I think the older code has the solar day more or less constant.

On Tue, May 2, 2023, 11:39 AM robmod @.***> wrote:

Darn. Well I for one appreciate the work you're putting in to get this code running! In the meantime, for those of us who aren't waiting on the new code to try high-E models, do you have an idea in what ways the native Plasim code used in EPS will prove inaccurate when doing high E models? Or how inaccurate?

Regarding my earlier issue (with the code before this eccentricity update)...I think I figured it out. Lonvernaleq is accurately setting angle between perihelion and autumn solstice (documentation says spring solstice...but that does not look correct). However changing rotation rate often changes at what time the simulation begins, and therefore changes the order of the output data. I had assumed that the beginning month was held constant, and so when I examined the data it looked like the seasons were getting messed up.

This quirk doesn't seem to be effecting the accuracy of the modeling, but it does make an exact perfect apple-to-apple comparison of data between some models with different rotation values more difficult (or not possible) since the snapshots/averaged months aren't representing the exact same time periods. For example, one model may (for sake of argument) begin the sim on January 1st, while another may begin it on August 13th...so the 12 averaged months in each model represent different time periods offset by 13 days.

Setting eccentricity with the old code does not seem to mess with the above. Zero Lonvernaleq is still autumn solstice=perihelion at any eccentricity.

— Reply to this email directly, view it on GitHub https://github.com/alphaparrot/ExoPlaSim/pull/2#issuecomment-1531695088, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADRQKQBB37SIR2NKTFA6IQTXEETEDANCNFSM6AAAAAAQVO43HE . You are receiving this because you were mentioned.Message ID: @.***>

[hersfeldtn]

I think the spring/autumn thing with lonvernaleq is just a quirk of how longitude of vernal equinox is defined: it is the point of a planet's orbit which is on the far side of the sun when Earth (or other reference planet) is at vernal equinix, which for the Earth itself is the autumn equinox.

I've also run some tests before and it looks like the current model treats the solar and sidereal days as equivalent, or very similar; if you set the rotation to 360 days with a 360-day year, you'll get 1 solar day per year, judging by the temperature and czen outputs

[alphaparrot]

Yes, the math of how solar day is computed assumes the sidereal year is much longer than the sidereal day; it's why tidally-locked models have a special logical switch in the code, rather than just setting day length equal to year length. I've never been a fan of it; going to a more fundamental code from which day length is an emergent property of the orbit and rotation rate will make me very happy.

On Tue, May 2, 2023, 1:20 PM hersfeldtn @.***> wrote:

I think the spring/autumn thing with lonvernaleq is just a quirk of how longitude of vernal equinox is defined: it is the point of a planet's orbit which is on the far side of the sun when Earth (or other reference planet) is at vernal equinix, which for the Earth itself is the autumn equinox.

I've also run some tests before and it looks like the current model treats the solar and sidereal days as equivalent, or very similar; if you set the rotation to 360 days with a 360-day year, you'll get 1 solar day per year, judging by the temperature and czen outputs

— Reply to this email directly, view it on GitHub https://github.com/alphaparrot/ExoPlaSim/pull/2#issuecomment-1531856824, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADRQKQD5EFT5BNIRTQO6R5LXEE65DANCNFSM6AAAAAAQVO43HE . You are receiving this because you were mentioned.Message ID: @.***>

[alphaparrot]

Hey folks, brief update: I've released version 3.2.0, which includes the (now-fixed) new orbital code. By my tests, it's working correctly and will accurately handle sun position in the sky for just about any orbital configuration. Note that you may still encounter crashes if the system you give it results in a sufficiently extreme climate (for example, e=0.7 implies an orbit whose insolation spans from 34.6% the insolation at the semimajor axis to over eleven times the insolation at semimajor axis). There's one final set of guard rails still active, and those limit the obliquity to prograde obliquities (+/- 90 degrees). I'm prepared to take those guard rails off too, but need to do more testing to make sure the code handles retrograde rotation correctly (I suspect it doesn't, but changing that shouldn't be too difficult). Once I've taken those guard rails off I'll close this PR. High-obliquity planets (near 90 degrees or even 90 degrees) should be handled just fine even before those changes (I think).

A note on lonvernaleq (the documentation could probably use some cleaning here): A lot of confusion stems from the fact that we are using heliocentric and geocentric coordinate systems (that are left-handed and right-handed) simultaneously. As Nikolai noted, lonvernaleq is most properly understood to be the longitude of perihelion, which is related to the argument of perihelion, which is the angle of the periapsis measured from the vernal equinox, in heliocentric coordinates. The shift from heliocentric to geocentric introduces an offset of 180 degrees, which is why for Earth this value is 102.7 degrees, even though periapsis is much more than 102.7 degrees from the vernal equinox. It's all a horrendous mess of terminology. Look for some of the related terms to be cleaned up in the next version of ExoPlaSim, to make it easier to predict where your vernal equinox will be. IN GENERAL, lonvernaleq plus 180 degrees will be equal to the angle between the vernal equinox and perihelion.

Let me know if you come across any bugs related to these changes.

Other notable changes in 3.2.0 related to this PR:

• Output files now have time-wise output variables nu (true anomaly in degrees), zdec (solar declination angle in degrees), lambda (solar ecliptic longitude in degrees), and rdist (the planet-star distance modulus). For angles the time-averaging might be a bit wonky when e.g. 355 gets averaged with 5. I'll see if I can come up with a clever way to handle that.
• the inspect method should now allow to specify a list of years and concatenate their outputs together for a given variable, along the time axis (this may be buggy). Similarly, the spatialmath method in gcmt can now allow you to preserve the time axis when computing spatial averages by specifying time='all'.

[hersfeldtn]

I can confirm runs with up to 90 degrees obliquity work, you just might need to start them with either substantial ice cover or lower co2 (or flux) to prevent the model promptly crashing from the extreme summer temperatures at the poles.

[wizzeh]

That's so exciting! Thank you for your work!

[alphaparrot]

Hey folks, final update for this PR/feature request: I've fixed a minor error in the orbital math that would have caused problems for retrograde obliquities, can confirm that obliquity >90 or <-90 should now work just fine (at least with the keplerian mode), and have removed the guard rails, or rather, adjusted them; obliquity is now bounded by [-180,180] degrees, and eccentricity is bounded by [0,1). ExoPlaSim version 3.2.1post2 should be the relevant stable version to use for this. As a side effect, you should be able to get pretty much all expected emergent behavior just from orbital properties, such as tidal-locking, 3:2 resonances, etc. Although for 1:1 synchronicity, it is still safest to use the existing tidally-locked flags and options--there is a possibility for the accumulation of machine error to cause the substellar longitude to drift when rotation period and sidereal orbital period are the same. Since if I'm not mistaken this satisfies the spirit of this PR, I will be closing this PR.