Eatdirt: Implement parallax effects for the Moon at the rendering stage and not with the ephemeris.

simgear/ephemeris/ephemeris.cxx

@@ -68,7 +68,8 @@ SGEphemeris::~SGEphemeris( void ) {
 void SGEphemeris::update( double mjd, double lst, double lat ) {
     // update object positions
     our_sun->updatePosition( mjd );
-    moon->updatePosition( mjd, lst, lat, our_sun );
+    //    moon->updatePositionTopo( mjd, lst, lat, our_sun );
+    moon->updatePosition( mjd, our_sun );
     mercury->updatePosition( mjd, our_sun );
     venus->updatePosition( mjd, our_sun );
     mars->updatePosition( mjd, our_sun );

simgear/ephemeris/ephemeris.hxx

@@ -146,6 +146,11 @@ public:
 	return moon->getDeclination();
     }
 
+    /** @return the geocentric distance to the Moon in unit of its semi-mayor axis. */
+    inline double getMoonDistanceInMayorAxis() const {
+      return moon->getDistanceInMayorAxis();
+    }
+
     /** @return the numbers of defined planets. */
     inline int getNumPlanets() const { return nplanets; }
 

simgear/ephemeris/moonpos.cxx

@@ -70,12 +70,14 @@ MoonPos::~MoonPos()
 
 
 /*****************************************************************************
- * void MoonPos::updatePosition(double mjd, Star *ourSun)
+ * void MoonPos::updatePositionTopo(double mjd, double lst, double
- * this member function calculates the actual topocentric position (i.e.) 
+ * lat, Star *ourSun) this member function calculates the actual
- * the position of the moon as seen from the current position on the surface
+ * topocentric position (i.e.)  the position of the moon as seen from
- * of the moon. 
+ * the current position on the surface of the earth. This include
+ * parallax effects, the moon appears on different stars background
+ * for different observers on the surface of the earth.
  ****************************************************************************/
-void MoonPos::updatePosition(double mjd, double lst, double lat, Star *ourSun)
+void MoonPos::updatePositionTopo(double mjd, double lst, double lat, Star *ourSun)
 {
   double
     eccAnom, ecl, actTime,
@@ -154,6 +156,7 @@ void MoonPos::updatePosition(double mjd, double lst, double lat, Star *ourSun)
 	-0.46 * cos(twoD)
 	);
   distance = r;
+  distance_in_a = r/a;
   // SG_LOG(SG_GENERAL, SG_INFO, "Running moon update");
   rcoslatEcl = r * cos(latEcl);
   xg =     cos(lonEcl) * rcoslatEcl;
@@ -234,3 +237,135 @@ void MoonPos::updatePosition(double mjd, double lst, double lat, Star *ourSun)
   I_factor = (log_I - max_loglux) / (max_loglux - min_loglux) + 1.0;
   I_factor = SGMiscd::clip(I_factor, 0, 1);
 }
+
+
+
+/*****************************************************************************
+ * void MoonPos::updatePosition(double mjd, Star *ourSun) this member
+ * function calculates the geocentric position (i.e.)  the position of
+ * the moon as seen from the center of earth. As such, it does not
+ * include any parallax effects. These are taken into account during
+ * the rendering.
+ ****************************************************************************/
+void MoonPos::updatePosition(double mjd, Star *ourSun)
+{
+  double
+    eccAnom, ecl, actTime,
+    xv, yv, v, r, xh, yh, zh, zg, xe,
+    Ls, Lm, D, F, geoRa, geoDec,
+    cosN, sinN, cosvw, sinvw, sinvw_cosi, cosecl, sinecl, rcoslatEcl,
+    FlesstwoD, MlesstwoD, twoD, twoM, twolat, alpha;
+
+  double max_loglux = -0.504030345621;
+  double min_loglux = -4.39964634562;
+  double conv = 1.0319696543787917;    // The log foot-candle to log lux conversion factor.
+  updateOrbElements(mjd);
+  actTime = sgCalcActTime(mjd);
+
+  // calculate the angle between ecliptic and equatorial coordinate system
+  // in Radians
+  ecl = SGD_DEGREES_TO_RADIANS * (23.4393 - 3.563E-7 * actTime);
+  eccAnom = sgCalcEccAnom(M, e);  // Calculate the eccentric anomaly
+  xv = a * (cos(eccAnom) - e);
+  yv = a * (sqrt(1.0 - e*e) * sin(eccAnom));
+  v = atan2(yv, xv);               // the moon's true anomaly
+  r = sqrt (xv*xv + yv*yv);       // and its distance
+
+  // repetitive calculations, minimised for speed
+  cosN = cos(N);
+  sinN = sin(N);
+  cosvw = cos(v+w);
+  sinvw = sin(v+w);
+  sinvw_cosi = sinvw * cos(i);
+  cosecl = cos(ecl);
+  sinecl = sin(ecl);
+
+  // estimate the geocentric rectangular coordinates here
+  xh = r * (cosN * cosvw - sinN * sinvw_cosi);
+  yh = r * (sinN * cosvw + cosN * sinvw_cosi);
+  zh = r * (sinvw * sin(i));
+
+  // calculate the ecliptic latitude and longitude here
+  lonEcl = atan2 (yh, xh);
+  latEcl = atan2(zh, sqrt(xh*xh + yh*yh));
+
+  /* Calculate a number of perturbation, i.e. disturbances caused by the
+   * gravitational influence of the sun and the other major planets.
+   * The largest of these even have a name */
+  Ls = ourSun->getM() + ourSun->getw();
+  Lm = M + w + N;
+  D = Lm - Ls;
+  F = Lm - N;
+
+  twoD = 2 * D;
+  twoM = 2 * M;
+  FlesstwoD = F - twoD;
+  MlesstwoD = M - twoD;
+
+  lonEcl += SGD_DEGREES_TO_RADIANS * (-1.274 * sin(MlesstwoD)
+			  +0.658 * sin(twoD)
+			  -0.186 * sin(ourSun->getM())
+			  -0.059 * sin(twoM - twoD)
+			  -0.057 * sin(MlesstwoD + ourSun->getM())
+			  +0.053 * sin(M + twoD)
+			  +0.046 * sin(twoD - ourSun->getM())
+			  +0.041 * sin(M - ourSun->getM())
+			  -0.035 * sin(D)
+			  -0.031 * sin(M + ourSun->getM())
+			  -0.015 * sin(2*F - twoD)
+			  +0.011 * sin(M - 4*D)
+			  );
+  latEcl += SGD_DEGREES_TO_RADIANS * (-0.173 * sin(FlesstwoD)
+			  -0.055 * sin(M - FlesstwoD)
+			  -0.046 * sin(M + FlesstwoD)
+			  +0.033 * sin(F + twoD)
+			  +0.017 * sin(twoM + F)
+			  );
+  r += (-0.58 * cos(MlesstwoD)
+	-0.46 * cos(twoD)
+	);
+  // from the orbital elements, the unit of the distance is in Earth radius, around 60.
+  distance = r;
+  distance_in_a = r/a;
+  // SG_LOG(SG_GENERAL, SG_INFO, "Running moon update");
+  rcoslatEcl = r * cos(latEcl);
+  xg =     cos(lonEcl) * rcoslatEcl;
+  yg =     sin(lonEcl) * rcoslatEcl;
+  zg = r *               sin(latEcl);
+
+  xe = xg;
+  ye = yg * cosecl -zg * sinecl;
+  ze = yg * sinecl +zg * cosecl;
+
+  geoRa  = atan2(ye, xe);
+  geoDec = atan2(ze, sqrt(xe*xe + ye*ye));
+
+  if (geoRa < 0)
+    geoRa += SGD_2PI;
+
+  rightAscension = geoRa;
+  declination = geoDec;
+
+
+  /* SG_LOG( SG_GENERAL, SG_INFO,
+	  "Ra = (" << (SGD_RADIANS_TO_DEGREES *rightAscension)
+	  << "), Dec= (" << (SGD_RADIANS_TO_DEGREES *declination) << ")" ); */
+
+  // Moon age and phase calculation
+  age = lonEcl - ourSun->getlonEcl();
+  phase = (1 - cos(age)) / 2;
+
+  // The log of the illuminance of the moon outside the atmosphere.
+  // This is the base 10 log of equation 20 from Krisciunas K. and Schaefer B.E.
+  // (1991). A model of the brightness of moonlight, Publ. Astron. Soc. Pacif.
+  // 103(667), 1033-1039 (DOI: http://dx.doi.org/10.1086/132921).
+  alpha = SGD_RADIANS_TO_DEGREES * SGMiscd::normalizeAngle(age + SGMiscd::pi());
+  log_I = -0.4 * (3.84 + 0.026*fabs(alpha) + 4e-9*pow(alpha, 4.0));
+
+  // Convert from foot-candles to lux.
+  log_I += conv;
+
+  // The moon's illuminance factor, bracketed between 0 and 1.
+  I_factor = (log_I - max_loglux) / (max_loglux - min_loglux) + 1.0;
+  I_factor = SGMiscd::clip(I_factor, 0, 1);
+}

simgear/ephemeris/moonpos.hxx

@@ -39,6 +39,7 @@ private:
     double xg, yg;     // the moon's rectangular geocentric coordinates
     double ye, ze;     // the moon's rectangular equatorial coordinates
     double distance;   // the moon's distance to the earth
+    double distance_in_a;  // the moon's distance to the earth in unit of its semi-mayor axis a
     double age;        // the moon's age from 0 to 2pi
     double phase;      // the moon's phase
     double log_I;      // the moon's illuminance outside the atmosphere (logged)
@@ -59,7 +60,8 @@ public:
     MoonPos(double mjd);
     MoonPos();
     ~MoonPos();
-    void updatePosition(double mjd, double lst, double lat, Star *ourSun);
+    void updatePositionTopo(double mjd, double lst, double lat, Star *ourSun);
+    void updatePosition(double mjd, Star *ourSun);
   // void newImage();
     double getM() const;
     double getw() const;
@@ -68,6 +70,7 @@ public:
     double getye() const;
     double getze() const;
     double getDistance() const;
+    double getDistanceInMayorAxis() const;
     double getAge() const;
     double getPhase() const;
     double getLogIlluminance() const;
@@ -109,6 +112,11 @@ inline double MoonPos::getDistance() const
   return distance;
 }
 
+inline double MoonPos::getDistanceInMayorAxis() const
+{
+  return distance_in_a;
+}
+
 inline double MoonPos::getAge() const
 {
   return age;

simgear/math/simd.hxx

@@ -82,7 +82,7 @@ inline T magnitude2(const simd4_t<T,N>& vi) {
 
 template<typename T, int N>
 inline simd4_t<T,N> interpolate(T tau, const simd4_t<T,N>& v1, const simd4_t<T,N>& v2) {
-    return (T(1)-tau)*v1 + tau*v2;
+    return v1 + tau*(v2-v1);
 }
 
 template<typename T, int N>

simgear/scene/sky/moon.cxx

@@ -69,7 +69,7 @@ SGMoon::~SGMoon( void ) {
 osg::Node*
 SGMoon::build( SGPath path, double moon_size ) {
 
-    osg::Node* orb = SGMakeSphere(moon_size, 15, 15);
+    osg::Node* orb = SGMakeSphere(moon_size, 40, 20);
     osg::StateSet* stateSet = orb->getOrCreateStateSet();
     stateSet->setRenderBinDetails(-5, "RenderBin");
 
@@ -162,22 +162,74 @@ bool SGMoon::repaint( double moon_angle ) {
 
 
 // reposition the moon at the specified right ascension and
-// declination, offset by our current position (p) so that it appears
+// declination from the center of Earth. Because the view is actually
-// fixed at a great distance from the viewer.  Also add in an optional
+// offset by our current position (p), we first evaluate our current
-// rotation (i.e. for the current time of day.)
+// position w.r.t the Moon and then shift to the articial center of
-bool SGMoon::reposition( double rightAscension, double declination,
+// earth before shifting to the rendered moon distance. This allows to
-			 double moon_dist )
+// implement any parallax effects.  moon_dist_bare is expected to not
-{
+// change during the rendering, it gives us the normalisation factors
-    osg::Matrix T2, RA, DEC;
+// between real distances and units used in the
+// rendering. moon_dist_fact is any extra factors to put the moon
+// further or closer.
 
+
+bool SGMoon::reposition( double rightAscension, double declination,
+			 double moon_dist_bare, double moon_dist_factor,
+			 double lst, double lat, double alt )
+{
+    osg::Matrix TE, T2, RA, DEC;
+
+    // semi mayor axis
+    const double moon_a_in_rearth = 60.266600;
+    // average (we could also account for equatorial streching like in
+    // oursun.cxx if required)
+    const double earth_radius_in_meters = 6371000.0;
+    //local hour angle in radians
+    double lha;
+    // in unit of the rendering
+    double moon_dist;
+    double earth_radius, viewer_radius;
+    double xp,yp,zp;
+
+    //rendered earth radius according to what has been specified by
+    //moon_dist_bare
+    earth_radius = moon_dist_bare/moon_a_in_rearth;
+
+    //how far are we from the center of Earth
+    viewer_radius = (1.0 + alt/earth_radius_in_meters)*earth_radius;
+
+    // the local hour angle of the moon, .i.e. its angle with respect
+    // to the meridian of the viewer
+    lha = lst - rightAscension;
+
+    // the shift vector of the observer w.r.t. earth center (funny
+    // convention on x?)
+    zp = viewer_radius * sin(lat);
+    yp = viewer_radius * cos(lat)*cos(lha);
+    xp = viewer_radius * cos(lat)*sin(-lha);
+
+    //rotate along the z axis
     RA.makeRotate(rightAscension - 90.0 * SGD_DEGREES_TO_RADIANS,
                   osg::Vec3(0, 0, 1));
-
+    //rotate along the rotated x axis
     DEC.makeRotate(declination, osg::Vec3(1, 0, 0));
 
+    //move to the center of Earth
+    TE.makeTranslate(osg::Vec3(-xp,-yp,-zp));
+
+    //move the moon from the center of Earth to moon_dist
+    moon_dist = moon_dist_bare * moon_dist_factor;
     T2.makeTranslate(osg::Vec3(0, moon_dist, 0));
 
-    moon_transform->setMatrix(T2*DEC*RA);
+    // cout << " viewer radius= " << viewer_radius << endl;
+    // cout << " xp yp zp= " << xp <<" " << yp << " " << zp << endl;
+    // cout << " lha= " << lha << endl;
+    // cout << " moon_dist_bare= " << moon_dist_bare << endl;
+    // cout << " moon_dist_factor= " << moon_dist_factor << endl;
+    // cout << " moon_dist= " << moon_dist << endl;
+
+    moon_transform->setMatrix(T2*TE*DEC*RA);
+
 
     return true;
 }

simgear/scene/sky/moon.hxx

@@ -65,11 +65,12 @@ public:
     bool repaint( double moon_angle );
 
     // reposition the moon at the specified right ascension and
-    // declination, offset by our current position (p) so that it
+    // declination, at moon_dist_bare*moon_dist_factor from the center
-    // appears fixed at a great distance from the viewer.  Also add in
+    // of Earth. lst, lat and alt are need to estimate where is the
-    // an optional rotation (i.e. for the current time of day.)
+    // center of Earth from the current view).
     bool reposition( double rightAscension, double declination,
-		     double moon_dist  );
+		     double moon_dist_bare, double moon_dist_factor,
+		     double lst, double lat, double alt  );
 };
 
 

simgear/scene/sky/sky.cxx

@@ -152,7 +152,7 @@ bool SGSky::reposition( const SGSkyState &st, const SGEphemeris& eph, double dt
     double angleRad = SGMiscd::deg2rad(angle);
 
     SGVec3f zero_elev, view_up;
-    double lon, lat, alt;
+    double lon, lat, alt, lst;
 
     SGGeod geodZeroViewPos = SGGeod::fromGeodM(st.pos_geod, 0);
     zero_elev = toVec3f( SGVec3d::fromGeod(geodZeroViewPos) );
@@ -165,6 +165,8 @@ bool SGSky::reposition( const SGSkyState &st, const SGEphemeris& eph, double dt
     lon = st.pos_geod.getLongitudeRad();
     lat = st.pos_geod.getLatitudeRad();
     alt = st.pos_geod.getElevationM();
+    // Local sidereal time
+    lst = angleRad + lon;
     
     dome->reposition( zero_elev, alt, lon, lat, st.spin );
 
@@ -178,7 +180,12 @@ bool SGSky::reposition( const SGSkyState &st, const SGEphemeris& eph, double dt
 
     double moon_ra = eph.getMoonRightAscension();
     double moon_dec = eph.getMoonDeclination();
-    moon->reposition( moon_ra, moon_dec, st.moon_dist );
+    double moon_r = eph.getMoonDistanceInMayorAxis();
+
+    //this allows to render the moon closer to the viewer when the
+    //moon is closer to the center of Earth, times any articial extra factors
+    double moon_dist_factor = moon_r * st.moon_dist_factor;
+    moon->reposition( moon_ra, moon_dec, st.moon_dist_bare, moon_dist_factor, lst, lat, alt );
 
     for ( unsigned i = 0; i < cloud_layers.size(); ++i ) {
         if ( cloud_layers[i]->getCoverage() != SGCloudLayer::SG_CLOUD_CLEAR ||

simgear/scene/sky/sky.hxx

@@ -63,7 +63,8 @@ struct SGSkyState
   double gst;      //!< GMT side real time.
   double sun_dist; //!< the sun's distance from the current view point
                    //   (to keep it inside your view volume).
-  double moon_dist;//!< The moon's distance from the current view point.
+  double moon_dist_bare ;//!< The moon's semi-mayor axis in the rendering (constant)
+  double moon_dist_factor ;//!< Any factor that are needed to artificially change the moon distance
   double sun_angle;
 };