Lunar Eclipse Scripting Question

[NotSoRandomOne]

This is a question that turns into a feature request if it doesn't exist. The feature request may be related to "Display of Earth's Shadow" request, at #430.

Does Stellarium expose any way to tell if the sun and moon are in opposition, or close to it, via scripting? For solar eclipses, getObjectInfo can be used to find the sun and moon's RA and Dec, to find their separation via formula. If they are less than half a degree or so, they are eclipsing.

But for lunar eclipses, the center of the earth must be taken into account, so RA and Dec for a given location, in conjunction with separation via formula, won't often work, unless the observer's location happens to be close to inline with the sun/earth vector.

Now that I think about it more, it is possible to check for separation > 177.5 or so, and get the moon's apparent magnitude. Once it drops below a certain value, you know an eclipse is occurring. Is that the only way?

Thanks,
David

PS - I don't know if it is just my system, but executing the 'solar_eclipse.ssc' script included with Stellarium in the 19.1.17113 build results in serious screen flicker, like double-buffering isn't occurring properly, during the heart of the eclipse event.

[gzotti]

There are eclipse-finding formulae which are still (only) in the literature...

[NotSoRandomOne]

If you add a function to locate the center of the sun's shadow, it makes adopting this script for lunar eclipse use a piece of cake. (This one is fun to watch, at least for me!)

// ******************************************************************************************
// Stellarium script to find the next next solar eclipse, to within specifiable radian
// difference of 'maxDiff'.
//
// Created by David O'Neil, Aug, 2019. Released to the public domain, for anyone interested.
// Just leave this 'created by' (or 'originally created by') verbage in the derived work.
//
// To use:
//   * Set the date via Stellarium's user interface close to a time you are looking for a
//     solar eclipse to occur.
//   * To go backwards in time, change the 'dayDiff' variable accordingly.
//   * Run the script. It will stop executing when the 'maxDiff' condition is met.
// ******************************************************************************************

//Set the following to "-29 days" to go backwards:
dayDiff = "+29 days";
maxDiff = 0.004;

//Don't change anything below here:
obj1 = "Sun";
obj2 = "Moon";
a1 = 0.0;
a2 = 1.0;
pi = Math.PI;
degToRad = pi/180;

core.selectObjectByName(obj1, true);
StelMovementMgr.setFlagTracking(true);
justFindSeparation = false;
diff = 1.0;
if (justFindSeparation) {
   d = getRadiansBetweenObjects();
   d = d/degToRad;
   degrees = Math.floor(d);
   minutes = d - degrees;
   minutes = minutes * 60;
   seconds = minutes - Math.floor(minutes);
   core.debug(degrees + " " + Math.floor(minutes) + " " + seconds*60);
   }
else {
   while (diff > maxDiff) {
      findClosestConjunction();
      diff = getRadiansBetweenObjects();
      core.debug(core.getDate() + ", Angle Diff: " + diff);
      if (diff >= maxDiff) core.setDate(dayDiff);
      }
   core.debug("Success!");
   core.debug(core.getDate());
   core.debug("Angle Diff: " + a2);
   }



function findClosestConjunction() {
   core.selectObjectByName(obj1, true);
   StelMovementMgr.setFlagTracking(true);
   core.setTimeRate(0); //Lock the sky in place, so we can see it unmoving!

   //Get the initial 2 points:
   a1 = getRadiansBetweenObjects();

   //Do initial iteration in 1-day increment:
   a2 = getTheAngleAfter("+1 days");

   direction = "+";
   if (a2 < a1) {
      direction = "+";
      }
   else {
      direction = "-";
      a2 = getTheAngleAfter("-1 days");
      a2 = getTheAngleAfter("-1 days");   //Don't know why can't use "-2 days" reliably!
      }
      
   timeStepString = direction + "1 days";

   //Iterate:
   while (a2 < a1) {
      a1 = a2;
      a2 = getTheAngleAfter(timeStepString);
      }
      
   //Assume we've gone through several days, and have passed the last good one,
   //so step back (or forward) a day to fix the overstep:
   if (direction == "+") {
      a1 = getTheAngleAfter("-1 days");
      }
   else {
      a1 = getTheAngleAfter("+1 days");
      }
      

   //Now we are doing it by the hours. As before, we must check whether we need to
   //go back or forward.

   //Now determine whether to go forward or back:
   a2 = getTheAngleAfter("+1 hour");

   if (a2 < a1) {
      direction = "+";
      }
   else {
      direction = "-";
      a2 = getTheAngleAfter("-1 hours");
      a2 = getTheAngleAfter("-1 hours");  //As before, don't know why can't just -2 this.
      }

   timeStepString = direction + "1 hour";

   while (a2 < a1) {
      a1 = a2;
      a2 = getTheAngleAfter(timeStepString);
      }

   //Assume we've gone through several hours, and have passed the last good one,
   //so step back (or forward) an hour to fix the overstep:
   if (direction == "+") {
      a1 = getTheAngleAfter("-1 hours");
      }
   else {
      a1 = getTheAngleAfter("+1 hours");
      }



   //And finally, for the minutes:
   //First, set a1 to a2, as a2 was last measurement taken at current time:
   a1 = a2;
   a2 = getTheAngleAfter("+1 minute");

   if (a2 < a1) {
      direction = "+";
      }
   else {
      direction = "-";
      a2 = getTheAngleAfter("-1 minutes");
      a2 = getTheAngleAfter("-1 minutes");
      }

   timeStepString = direction + "1 minute";

   while (a2 < a1) {
      a1 = a2;
      a2 = getTheAngleAfter(timeStepString);
      }

   core.setTimeRate(0); //Lock the sky in place, so we can see it unmoving!
   }


function getTheAngleAfter(timeStep) {
   obj2Azi = core.getObjectInfo(obj2).azimuth;
   core.setDate(timeStep);
   while (core.getObjectInfo(obj2).azimuth == obj2Azi) {
      core.wait(0.01);
      }
   angleDiff = getRadiansBetweenObjects();
   return angleDiff;
   }


function getRadiansBetweenObjects() {
   obj2Azi = core.getObjectInfo(obj2).azimuth;
   obj2Azi = obj2Azi * degToRad;
   obj2Alt = core.getObjectInfo(obj2).altitude;
   obj2Alt = obj2Alt * degToRad;

   obj1Azi = core.getObjectInfo(obj1).azimuth;
   obj1Azi = obj1Azi * degToRad;
   obj1Alt = core.getObjectInfo(obj1).altitude;
   obj1Alt = obj1Alt * degToRad;

   d = Math.acos(Math.sin(obj2Alt)*Math.sin(obj1Alt)+Math.cos(obj2Alt)*Math.cos(obj1Alt)*Math.cos(obj2Azi-obj1Azi));
   
   //The following formula is for when the objects are very close together, or very close to 180 degrees apart:
   //I may have Azi and Alt mixed up.
   
   // delAzi = obj1Azi - obj2Azi;
   // delAlt = obj1Alt - obj2Alt;

   //d = sqrt ( (cos(aziA-aziB) * (altA-altB)) ^2 + (aziA-aziB)^2 )
   //d = Math.sqrt ( Math.cos(delAzi)*(delAlt) * Math.cos(delAzi)*(delAlt) + (delAzi)*(delAzi) );
   
   return d;   //In radians
   }

[alex-w]

@NotSoRandomOne are you using Stellarium on Linux?

[NotSoRandomOne]

@alex-w No. Win10.

[NotSoRandomOne]

@alex-w I have thought about trying to build Stellarium on my system, and add a function to allow grabbing the earth's shadow. Your code is clean enough I would feel comfortable doing so, but I haven't had the time to master Git.

[NotSoRandomOne]

Can a script get the text that 'Control-C' currently copies into the clipboard? After lots of digging, I see no easy way to determine if a lunar eclipse is occurring via script.

The simplest solution in my mind right now is to somehow access the 'Control - C' text when the sun is selected, search it for 'Eclipse magnitude', and then take appropriate action if it is found.

But 'getSelectedObjectInfo' seems to be the closest action to 'Ctrl-C,' and they do different things! So there doesn't seem to be a way by script to get whether 'Eclipse obscuration/magnitude' are being displayed.

[alex-w]

Please check a fresh beta - https://github.com/Stellarium/stellarium-data/releases/tag/beta

[NotSoRandomOne]

Alex,

I got a script put together for lunar eclipses, and was wondering why the new build wasn't showing me the information. Realized you got it to add the info for solar eclipses, not lunar ones! In other words, we need eclipse info for the moon.

Here's the script I've got so far, if you want to play with it. It isn't complete, but will stop when it finds a key named "eclipse-magnitude" in the map.

Thanks,
David

// ******************************************************************************************
// This is currently the best Stellarium can do for finding lunar eclipses.
// Go the next full moon, and let the observer decide if it is an eclipse:
// ******************************************************************************************

"use strict";

var moon = "Moon";
var sun = "Sun";

var secondsInMoonOrbit = 29.53 * 60 * 60 * 24;

var stepForward = false;

if (stepForward) separation = getTheSeparationAfter("+29 days");
else separation = getTheSeparationAfter("-29 days");

core.setTimeRate(0);

var separation = sunMoonEclipticSeparation();

if (separation > 180 && stepForward) {
   var degToMove = 360 - separation + 180;
   var daysToMove = Math.floor(degToMove / 360 * 29.5);
   if (daysToMove != 0) separation = getTheSeparationAfter("+" + daysToMove + " days");
   }
else if (separation < 180 && !stepForward) {
   var degToMove = 180 + separation;
   var daysToMove = Math.floor(degToMove / 360 * 29.5);
   if (daysToMove != 0) separation = getTheSeparationAfter("-" + daysToMove + " days");
   }
else {
   var degToMove = 180 - separation;
   var daysToMove = Math.floor(degToMove / 360 * 29.5);
   if (daysToMove > 0) separation = getTheSeparationAfter("+" + daysToMove + " days");
   else if (daysToMove < 0) separation = getTheSeparationAfter("-" + daysToMove + " days");
   }

var foundEclipse = false;

//while (!foundEclipse) {
   //Now find opposition:
   var secondsToMove = 100;   //Initial dummy value
   do {
      var degToMove = 180 - separation;
      var percentToMove = degToMove / 360;
      secondsToMove = Math.floor(percentToMove * secondsInMoonOrbit);
      if (secondsToMove == 0) break;
      core.debug(secondsToMove);
      //core.pauseScript();
      if (secondsToMove > 0) separation = getTheSeparationAfter("+" + secondsToMove + " seconds");
      else {
         separation = getTheSeparationAfter("-" + Math.abs(secondsToMove) + " seconds");
         }
      }
      while (secondsToMove > 1);
      
   // //Now see if we have a lunar eclipse:
   // core.selectObjectByName("Moon", true);
   // map=core.getSelectedObjectInfo();

   // var t = map["eclipse-magnitude"];
   // if (t !== undefined) {
      // foundEclipse = true;
      // core.debug("Here!");
      // }
   // else {
      // if (stepForward) separation = getTheSeparationAfter("+29 days");
      // else separation = getTheSeparationAfter("-29 days");
      // }
      
      
      
//   }
   
   core.selectObjectByName("Moon", true);
   StelMovementMgr.setFlagTracking(true);

   
   

function getTheSeparationAfter(timeStep) {
   var initVal = core.getObjectInfo("Moon").azimuth;
   timeStep = cleanupString(timeStep);
   core.setDate(timeStep);
   while (core.getObjectInfo("Moon").azimuth == initVal) {
      core.wait(0.01);
      }
   var result = sunMoonEclipticSeparation();
   //core.debug(result);
   return result;
   }


function cleanupString(timeStep) {
   if (timeStep.substring(0, 2) == "--") {
      var len = timeStep.length;
      return timeStep.substring(1, len);
      }
   return timeStep;
   }
   
   
function sunMoonEclipticSeparation() {
   //This function returns a value from 0 to 360, representing the angular ecliptic
   //between the sun and moon. Small values represent new moons, close to 360
   //represent old moons, getting ready to disappear into the morning glow.
   //180° represents full moons.
   var moonDelta = core.getObjectInfo(moon).elong - core.getObjectInfo(sun).elong
   while (moonDelta < 0) moonDelta += 360; //Returned in degrees.
   return moonDelta;
   }

[NotSoRandomOne]

(My previous script already works for solar eclipses without using the map.)

[alex-w]

@NotSoRandomOne We do not calculate eclipse magnitude and obscuration for lunar eclipses :-/ Could you suggest to me the maths/formulae for lunar eclipses stuff?

[NotSoRandomOne]

Looking at your source code, my first thought would be to add an 'eclipse-illumination' item to the map instead of obscuration/magnitude. You could update it in the 'SolarSystem::getEclipseFactor' function. Would that work? 1.00 = no eclipse, 0.0 is fully eclipsed?

[gzotti]

No, the eclipseFactor is only for Solar eclipses, where we can compute that easily with the geometries of the spherical bodies.

For Lunar eclipses, we need to compute the diameter of the shadow in the Moon's distance which should include an enlargement factor by Earth's atmosphere, which however differs among scientists (USNO vs. Danjon?). It has been on my longer-term list to track down the respective formulae, but this all takes time. At least the visual simulation quality is now much better than a year ago.

[NotSoRandomOne]

Where in the code is the lunar eclipse processing taking place? A quick look doesn't reveal the location to me. Whichever variables are controlling that processing is related to the attributes needed to check via scripting.

[gzotti]

Lunar position and Solar/Earth position, and both bodies' geometries. The rest is vector math. And with some original code still underdocumented, this makes life not easier. Read Planet.cpp and search for "eclipse".
I have no time for further hints this week, sorry.

[NotSoRandomOne]

Looking at the code briefly, the quickest thing that comes to mind is modifying the Planet::getCandidatesForShadow() function. Something like:

QVector<const Planet*> Planet::getCandidatesForShadow() const
{
   QVector<const Planet*> res;
   const SolarSystem *ssystem=GETSTELMODULE(SolarSystem);
   const Planet* sun = ssystem->getSun().data();
   if (this==sun || (parent.data()==sun && satellites.empty()))
      return res;
   
   for (const auto& planet : satellites)
   {
   
   
   
      //NEW:
      if (this == earth && planet==moon) {
         if (willCastShadow(this, planet.data())) {
            "lunar-eclipse-occurring" = true;   //Or however you change a map value
            res.append(planet.data());
         }
         else {
            "lunar-eclipse-occurring" = false;
         }
      }
      else {
         if (willCastShadow(this, planet.data()))
            res.append(planet.data());
      }
   }
   
   
   
   if (willCastShadow(this, parent.data()))
      res.append(parent.data());
   // Test satellites mutual occultations.
   if (parent.data() != sun)
   {
      for (const auto& planet : parent->satellites)
      {
         //skip self-shadowing
         if(planet.data() == this )
            continue;
         if (willCastShadow(this, planet.data()))
            res.append(planet.data());
      }
   }
   
   return res;
}

[NotSoRandomOne]

Of course, that is more pseudo-code than real, because I doubt I've done the correct thing to get the earth and moon, and the map value doesn't exist, and I haven't looked into how map values are created (which this one would have to be) or changed.

[NotSoRandomOne]

(It looks like OpenGL is taking care of performing the shadows? So Stellarium code is only checking if the shadow could occur, and putting it into OpenGL's shadowing routines. It may be possible to get rough info on the amount of lunar eclipse from the existing code in the Planet::drawSphere routine, through the 'push' variable:

   if (this==ssm->getMoon())
   {
      GL(normalMap->bind(2));
      GL(moonShaderProgram->setUniformValue(moonShaderVars.normalMap, 2));
      if (!rData.shadowCandidates.isEmpty())
      {
         GL(texEarthShadow->bind(3));
         GL(moonShaderProgram->setUniformValue(moonShaderVars.earthShadow, 3));
         // Ad-hoc visibility improvement during lunar eclipses:
         // During partial umbra phase, make moon brighter so that the bright limb and umbra border has more visibility.
         // When the moon is about half in umbra (geoc.elong 179.4), we start to raise its brightness.
         GLfloat push=1.0f;
         const double elong=getElongation(ssm->getEarth()->getEclipticPos()) * (180.0/M_PI);
         const float x=static_cast<float>(elong) - 179.5f;
         if (x>0.0f)
            push+=20.0f * x;
         if (x>0.1f)
            push=3.0f;

         GL(moonShaderProgram->setUniformValue(moonShaderVars.eclipsePush, push)); // constant for now...
      }

[NotSoRandomOne]

To give a percent value, I believe the 'willCastShadow' routine can be refactored to the following, and the 'penumbraTouchPercent' can be used. I don't know if any astronomic references have created such a usage, but I don't see why it couldn't be expressed that way:

bool willCastShadow(const Planet* shadowedPlanet, const Planet* castingPlanet)
{
   Vec3d shadowedPlanetPosition = shadowedPlanet->getHeliocentricEclipticPos();
   Vec3d castingPlanetPosition = castingPlanet->getHeliocentricEclipticPos();
   
   // If the planet castingPlanet is farther from the sun than this planet, it can't cast shadow on it.
   if (castingPlanetPosition.lengthSquared()>shadowedPlanetPosition.lengthSquared())
      return false;

   //Vec3d normalizedCastPosition = castingPlanetPosition();
   //normalizedCastPosition.normalize();
   //Can the previous two just be:
   Vec3d normalizedCastPosition = castingPlanetPosition.normalize();
   
   double shadowDistance = normalizedCastPosition * shadowedPlanetPosition;
   static const double sunRadius = 696000./AU;
   double castShadowDiameter = castingPlanetPosition.length() / (castingPlanet->getEquatorialRadius()/sunRadius+1);
   double penumbraRadius = (shadowDistance-castShadowDiameter)/castShadowDiameter*sunRadius;
   // TODO: Note that Earth's shadow should be enlarged a bit. (6-7% following Danjon?)
   
   double distanceBetweenPenumbraCenterAndShadowedPlanetCenter = (normalizedCastPosition*shadowDistance-shadowedPlanetPosition).length();
   
   double maxPenumbraOffset = penumbraRadius+shadowedPlanet->getEquatorialRadius(); //Penumbra just touching planet.
   
   if (distanceBetweenPenumbraCenterAndShadowedPlanetCenter < maxPenumbraOffset) {
      penumbraTouchPercent = 0.0;
      return false;
   }
   //If I've done the math correctly, you can express how close the penumbra is to the planetary center as a percentage:
   penumbraTouchPercent = (maxPenumbraOffset-distanceBetweenPenumbraCenterAndShadowedPlanetCenter)/maxPenumbraOffset;
   return true;
}

[gzotti]

I did not look at the code. But when we not only let OpenGL do its thing but want to give numerical values, we definitely must deal with penumbra and umbra sizes, shadow magnification and other subtleties. Probably also deal with visibility criteria (e.g. 70% in penumbra is needed for being just visible). A "deep" eclipse can be more than 100%, depending on various definitions. I don't want to invent our own workarounds here, but something comparable to common almanachs. These things should be available in the ephemeris literature, just did not find the time so far.

[NotSoRandomOne]

The last code listing has a refactoring that might be helpful understanding the code in the future. An older Fortran listing for calculating these charactisterics is in Appendix B of https://eclipse.gsfc.nasa.gov/SEpubs/RP1216.html.

edit - And I suspect the proper place to calculate those values is right where I added the penumbraTouchPercent variable.

[NotSoRandomOne]

I don't want to invent our own workarounds here, ...

That is understandable. Would you be agreeable to replacing my penumbraTouchPercent approach with a simple "lunar-eclipse-possible" variable in the meantime, so Stellarium can at least be queried via scripting for the possibility of lunar eclipses occurring? That change should be simple, and won't take the time the in-depth approach is going to require. It also brings some parity between solar an lunar eclipse scripting abilities to Stellarium.

[gzotti]

I can try to have a look at that topic hopefully next week. Maybe it's easy enough to solve faster than originally planned. But please write "penumbra" ;-)

[gzotti]

Actually I have written an eclipse finder for the HP48 pocket calculator, >20 years ago... This should one day be redone as plugin for Stellarium, of course.

[NotSoRandomOne]

But please write "penumbra" ;-)

Oops! I fixed them, but also found I had written it correctly a couple times! The joys of dyslexia! Don't know where that phoneticism came into being at...

HP48...

I still have my old 41CX in the drawer from my college engineering days (batteries out). Now I use Oliver de Smet's emulator on my phone: https://play.google.com/store/apps/details?id=o2s.emul.hp41cx. Liked it so much I bought it. It can even run the old programs!

Maybe it's easy enough to solve faster than originally planned...

Hope so! But I understand if not. Have a great day!

(I updated the previous lunar eclipse script to do all that can be done in Stellarium right now: advance to the next or previous full moon and let the user see it for themselves.)

[alex-w]

Hmm... the code for compute penumbral and umbral eclipse magnitudes are exist, but the data is not available in scripting engine.

[github-actions]

Hello @NotSoRandomOne! Please check the fresh version (development snapshot) of Stellarium:
https://github.com/Stellarium/stellarium-data/releases/tag/weekly-snapshot

[github-actions]

Hello @NotSoRandomOne! Please check the latest stable version of Stellarium:
https://github.com/Stellarium/stellarium/releases/latest