Physically correct Moon phases

Eugene

22-05-2012 18:20:19

Hello. Some time ago I was Caelum contributor (SpriteSun, root node, proper Sun direction and so on), but my needs were satisfied by Caelum 0.3. Today I decided to upgrade to Caelum 0.6.1, throw away custom SpriteSunWithOppositeMoon class that was not accepted into the main Caelum trunk and use embedded Moon support. To still have Moon opposite to the Sun I selected date of full moon, and it was located properly but not full at all. Choosing dates slightly before or after new moon date shows that phases calculation was totally wrong - lighted side of the moon was not aligned to the Sun, had wrong size and rotates with the camera if Moon and Sun were high in the sky.

I fixed phase calculation (synodic month duration was wrong), applied rotation axis to the moon`s billboard aligned with the real moon rotation axis to make bright side of the moon looking to the sun, replaced saw-shaped movement of terminator over strange phase range [0..2] with proper cosine law over phase range [0..1) where 0 is full moon, 1/2 is new moon and 1 is full moon again.

I can not attach patch here as "board attachment quota has been reached.", so lets go the hard way


# HG changeset patch
# User Eugene Golushkov <eugene@belightsoft.com>
# Date 1337687206 -10800
# Node ID efe480ea5fb0965fff902cd4636771f0517fc195
# Parent 04be98f28298822f2574ad9ff41d448edc46989e
Physically accurate Moon phase

diff -r 04be98f28298 -r efe480ea5fb0 main/include/Astronomy.h
--- a/main/include/Astronomy.h Mon Jan 30 05:40:56 2012 -0800
+++ b/main/include/Astronomy.h Tue May 22 14:46:46 2012 +0300
@@ -132,11 +132,16 @@
LongReal jday,
LongReal longitude, LongReal latitude,
LongReal &azimuth, LongReal &altitude);
- static void getHorizontalMoonPosition (
+ static void getHorizontalMoonPosition (
LongReal jday,
Ogre::Degree longitude, Ogre::Degree latitude,
Ogre::Degree &azimuth, Ogre::Degree &altitude);

+ static void getHorizontalNorthEclipticPolePosition (
+ LongReal jday,
+ Ogre::Degree longitude, Ogre::Degree latitude,
+ Ogre::Degree &azimuth, Ogre::Degree &altitude);
+
/** Get astronomical julian day from normal gregorian calendar.
* From wikipedia: the integer number of days that have elapsed
* since the initial epoch defined as
diff -r 04be98f28298 -r efe480ea5fb0 main/include/CaelumSystem.h
--- a/main/include/CaelumSystem.h Mon Jan 30 05:40:56 2012 -0800
+++ b/main/include/CaelumSystem.h Tue May 22 14:46:46 2012 +0300
@@ -636,14 +636,18 @@
*/
const Ogre::Vector3 getMoonDirection (LongReal jday);

- /** Fake function to get the phase of the moon
+ /** Function to get the phase of the moon
* @param jday Julian day
- * @return the phase of the moon; ranging from 0(full moon) to 2(new moon).
- * The calculations performed by this function are completely fake.
- * It's a triangle wave with a period of 28.5 days.
+ * @return the phase of the moon; ranging from 0(full moon) to 1(again full moon).
*/
const Ogre::Real getMoonPhase (LongReal jday);

+ /** Get the ecliptic's north pole direction at a certain time.
+ * Useful as Moon's north polar axis points within 1.5 degrees of the north ecliptic pole.
+ * @param jday astronomical julian day.
+ */
+ const Ogre::Vector3 getEclipticNorthPoleDirection (LongReal jday);
+
private:
/** Handle FrameListener::frameStarted to call updateSubcomponents every frame.
* If you don't register CaelumSystem as a an ogre frame listener you have to
diff -r 04be98f28298 -r efe480ea5fb0 main/include/Moon.h
--- a/main/include/Moon.h Mon Jan 30 05:40:56 2012 -0800
+++ b/main/include/Moon.h Tue May 22 14:46:46 2012 +0300
@@ -93,6 +93,9 @@

/// Set the moon's phase
void setPhase (Ogre::Real phase);
+
+ /// Set Moon`s north pole direction
+ void setMoonNorthPoleDirection(const Ogre::Vector3& moonNorthPoleDir);

public:
/// Handle camera change.
diff -r 04be98f28298 -r efe480ea5fb0 main/resources/CaelumPhaseMoon.cg
--- a/main/resources/CaelumPhaseMoon.cg Mon Jan 30 05:40:56 2012 -0800
+++ b/main/resources/CaelumPhaseMoon.cg Tue May 22 14:46:46 2012 +0300
@@ -20,26 +20,23 @@

// Get how much of a certain point on the moon is seen (or not) because of the phase.
// uv is the rect position on moon; as seen from the earth.
-// phase ranges from 0 to 2
+// phase ranges from 0 (full moon) to 1 (again fool moon)
float MoonPhaseFactor(float2 uv, float phase)
{
- float alpha = 1.0;
+ // 1. In phase interval [0..1/2) day-to-night terminator appeared on the right side of the moon and moves to the left by cosine law
+ // 2. In phase interval [1/2..1) night-to-day terminator appeared on the right side of the moon and moves to the left by cosine law
+ // but if we swapped for simplicity left and right sides of the moon for the second interval than
+ // 2' In phase interval [1/2..1) night-to-day terminator appeared on the left side of the moon and moves to the right by cosine law
+ // therefore in such coord system terminator would move from right to left and to the right again, and picture would be like this:
+ // Moon => (day | night)

- float srefx = uv.x - 0.5;
- float refx = abs(uv.x - 0.5);
- float refy = abs(uv.y - 0.5);
- float refxfory = sqrt(0.25 - refy * refy);
- float xmin = -refxfory;
- float xmax = refxfory;
- float xmin1 = (xmax - xmin) * (phase / 2) + xmin;
- float xmin2 = (xmax - xmin) * phase + xmin;
- if (srefx < xmin1) {
- alpha = 0;
- } else if (srefx < xmin2 && xmin1 != xmin2) {
- alpha = (srefx - xmin1) / (xmin2 - xmin1);
- }
-
- return alpha;
+ float cY = uv.y - 0.5; // signed
+ float cX = sqrt(0.25 - cY * cY); // positive, half of disk chord
+
+ float termX = cX * cos((2.0 * 3.1416) * phase); // determine terminator position in our tweaked coord system
+ float refX = (uv.x - 0.5) * sign(0.5 - phase); // reverse X axis for phase interval [1/2..1)
+
+ return 0.5 * sign(termX - refX) + 0.5; // day if refX < termX
}

void PhaseMoonFP
diff -r 04be98f28298 -r efe480ea5fb0 main/src/Astronomy.cpp
--- a/main/src/Astronomy.cpp Mon Jan 30 05:40:56 2012 -0800
+++ b/main/src/Astronomy.cpp Tue May 22 14:46:46 2012 +0300
@@ -238,6 +238,23 @@
altitude = Ogre::Degree(al);
}

+ void Astronomy::getHorizontalNorthEclipticPolePosition (
+ LongReal jday,
+ Ogre::Degree longitude, Ogre::Degree latitude,
+ Ogre::Degree &azimuth, Ogre::Degree &altitude)
+ {
+ // Get precise equatorial spherical coordinates of North Ecliptic Pole
+ LongReal rasc = 270.0;
+ LongReal decl = 90.0 - 23.439281;
+
+ // Equatorial to horizontal
+ LongReal az, al;
+ Astronomy::convertEquatorialToHorizontal (
+ jday, longitude.valueDegrees (), latitude.valueDegrees (), rasc, decl, az, al);
+ azimuth = Ogre::Degree(az);
+ altitude = Ogre::Degree(al);
+ }
+
int Astronomy::getJulianDayFromGregorianDate(
int year, int month, int day)
{
@@ -299,22 +316,22 @@
LongReal julianDay, int &year, int &month, int &day,
int &hour, int &minute, LongReal &second)
{
- // Integer julian days are at noon.
- // static_cast<int)(floor( is more precise than Ogre::Math::IFloor.
- // Yes, it does matter.
- int fpmode = enterHighPrecissionFloatingPointMode();
- julianDay += (LongReal)0.5;
- int ijd = static_cast<int>(floor(julianDay));
- getGregorianDateFromJulianDay(ijd, year, month, day);
-
-
- LongReal s = (julianDay - (LongReal)ijd);
- s *= 86400.0;
- hour = static_cast<int>(floor(s / 3600));
- s -= hour * 3600;
- minute = static_cast<int>(floor(s / 60));
- s -= minute * 60;
- second = s;
+ // Integer julian days are at noon.
+ // static_cast<int)(floor( is more precise than Ogre::Math::IFloor.
+ // Yes, it does matter.
+ int fpmode = enterHighPrecissionFloatingPointMode();
+ julianDay += (LongReal)0.5;
+ int ijd = static_cast<int>(floor(julianDay));
+ getGregorianDateFromJulianDay(ijd, year, month, day);
+
+
+ LongReal s = (julianDay - (LongReal)ijd);
+ s *= 86400.0;
+ hour = static_cast<int>(floor(s / 3600));
+ s -= hour * 3600;
+ minute = static_cast<int>(floor(s / 60));
+ s -= minute * 60;
+ second = s;
restoreFloatingPointMode(fpmode);
}

diff -r 04be98f28298 -r efe480ea5fb0 main/src/CaelumSystem.cpp
--- a/main/src/CaelumSystem.cpp Mon Jan 30 05:40:56 2012 -0800
+++ b/main/src/CaelumSystem.cpp Tue May 22 14:46:46 2012 +0300
@@ -514,6 +514,7 @@
moonDir,
moonLightColour,
moonBodyColour);
+ mMoon->setMoonNorthPoleDirection(getEclipticNorthPoleDirection(julDay)); // its not precise, but error is within 1.5 degrees
mMoon->setPhase (moonPhase);
}

@@ -737,6 +738,21 @@
return res;
}

+ const Ogre::Vector3 CaelumSystem::getEclipticNorthPoleDirection (LongReal jday)
+ {
+ Ogre::Degree azimuth, altitude;
+ {
+ ScopedHighPrecissionFloatSwitch precissionSwitch;
+
+ Astronomy::getHorizontalNorthEclipticPolePosition(jday,
+ getObserverLongitude (), getObserverLatitude (),
+ azimuth, altitude);
+ }
+ Ogre::Vector3 res = -makeDirection(azimuth, altitude);
+
+ return res;
+ }
+
const Ogre::Vector3 CaelumSystem::getMoonDirection (LongReal jday)
{
Ogre::Degree azimuth, altitude;
@@ -754,12 +770,13 @@

const Ogre::Real CaelumSystem::getMoonPhase (LongReal jday)
{
- // Calculates julian days since January 22, 2008 13:36 (full moon)
+ // Calculates julian days since June 04, 1993 20:31 (full moon)
// and divides by the time between lunations (synodic month)
- LongReal T = (jday - 2454488.0665L) / 29.531026L;
+ // Formula derived from http://scienceworld.wolfram.com/astronomy/Lunation.html
+ LongReal T = (jday - 2449143.3552943347L) / 29.53058867L;

- T = fabs(fmod(T, 1));
- return -fabs(-4 * T + 2) + 2;
+ T -= floor(T); // [0..1)
+ return T;
}

void CaelumSystem::forceSubcomponentQueryFlags (uint flags)
diff -r 04be98f28298 -r efe480ea5fb0 main/src/Moon.cpp
--- a/main/src/Moon.cpp Mon Jan 30 05:40:56 2012 -0800
+++ b/main/src/Moon.cpp Tue May 22 14:46:46 2012 +0300
@@ -118,6 +118,12 @@
mParams.phase.set(mParams.fpParams, phase);
}

+ void Moon::setMoonNorthPoleDirection(const Ogre::Vector3& moonNorthPoleDir)
+ {
+ mMoonBB->setBillboardType(Ogre::BBT_ORIENTED_COMMON);
+ mMoonBB->setCommonDirection(moonNorthPoleDir);
+ }
+
void Moon::setQueryFlags (uint flags) {
mMoonBB->setQueryFlags (flags);
mBackBB->setQueryFlags (flags);

duststorm

17-07-2012 20:26:21

Very nice!
Is there a reason it was never applied to the official repository? Did you ask forum member tdev whether this patch could be applied?

Eugene

14-11-2012 13:53:08

Very nice!
Is there a reason it was never applied to the official repository? Did you ask forum member tdev whether this patch could be applied?


OK, patch is applied.