Fix bug 182 - make sun direction agree with position of sun from the ephemeris

http://code.google.com/p/flightgear-bugs/issues/detail?id=182

There were two issues. The biggest is that the sunsolver was accessing
parameters from the ephemeris and doing a different calculation with
them to derive the right ascension and declination of the sun. I'm not
sure who is right, but I changed sunsolver to agree with the
ephemeris. Also, there was an inappropriate use of geodetic
coordinates in calculating  the sun latitude and longitude.

Also, I did some cleanup in updateSunPos().

diff --git a/src/Time/light.cxx b/src/Time/light.cxx
index 24d1c7caf1dd09e23a9ffa49e5fd19af151d7b82..348bd524daa3c8ee399b254cc2517e0fad0bbd40 100644 (file)
--- a/src/Time/light.cxx
+++ b/src/Time/light.cxx
@@ -48,38 +48,6 @@
 #include "light.hxx"
 #include "sunsolver.hxx"
 
-/**
- * Map i.e. project a vector onto a plane.
- * @param normal (in) normal vector for the plane
- * @param v0 (in) a point on the plane
- * @param vec (in) the vector to map onto the plane
- */
-static SGVec3f map_vec_onto_cur_surface_plane(const SGVec3f& normal, 
-                                             const SGVec3f& v0, 
-                                             const SGVec3f& vec)
-{
-    // calculate a vector "u1" representing the shortest distance from
-    // the plane specified by normal and v0 to a point specified by
-    // "vec".  "u1" represents both the direction and magnitude of
-    // this desired distance.
-
-    // u1 = ( (normal <dot> vec) / (normal <dot> normal) ) * normal
-    SGVec3f u1 = (dot(normal, vec) / dot(normal, normal)) * normal;
-
-    // calculate the vector "v" which is the vector "vec" mapped onto
-    // the plane specified by "normal" and "v0".
-
-    // v = v0 + vec - u1
-    SGVec3f v = v0 + vec - u1;
-    
-    // Calculate the vector "result" which is "v" - "v0" which is a
-    // directional vector pointing from v0 towards v
-
-    // result = v - v0
-    return v - v0;
-}
-
-
 // Constructor
 FGLight::FGLight ()
     : _ambient_tbl( NULL ),
@@ -432,16 +400,22 @@ void FGLight::updateSunPos()
     SG_LOG( SG_EVENT, SG_DEBUG, "  Gst = " << t->getGst() );
 
     double sun_l;
-    double sun_gd_lat;
-    fgSunPositionGST(t->getGst(), &sun_l, &sun_gd_lat);
+    double sun_gc_lat;
+    fgSunPositionGST(t->getGst(), &sun_l, &sun_gc_lat);
     set_sun_lon(sun_l);
-    set_sun_lat(sun_gd_lat);
-    SGVec3d sunpos(SGVec3d::fromGeod(SGGeod::fromRad(sun_l, sun_gd_lat)));
+    // It might seem that sun_gc_lat needs to be converted to geodetic
+    // latitude here, but it doesn't. The sun latitude is the latitude
+    // of the point on the earth where the up vector has the same
+    // angle from geocentric Z as the sun direction. But geodetic
+    // latitude is defined as 90 - angle of up vector from Z!
+    set_sun_lat(sun_gc_lat);
+    SGVec3d sunpos(SGVec3d::fromGeoc(SGGeoc::fromRadM(sun_l, sun_gc_lat,
+                                                      SGGeodesy::EQURAD)));
 
     SG_LOG( SG_EVENT, SG_DEBUG, "    t->cur_time = " << t->get_cur_time() );
     SG_LOG( SG_EVENT, SG_DEBUG,
-            "    Sun Geodetic lat = " << sun_gd_lat
-            << " Geodetic lat = " << sun_gd_lat );
+            "    Sun Geocentric lat = " << sun_gc_lat
+            << " Geodcentric lat = " << sun_gc_lat );
 
     // update the sun light vector
     sun_vec() = SGVec4f(toVec3f(normalize(sunpos)), 0);
@@ -450,8 +424,8 @@ void FGLight::updateSunPos()
     // calculate the sun's relative angle to local up
     SGVec3d viewPos = v->get_view_pos();
     SGQuatd hlOr = SGQuatd::fromLonLat(SGGeod::fromCart(viewPos));
-    SGVec3f world_up = toVec3f(hlOr.backTransform(-SGVec3d::e3()));
-    SGVec3f nsun = toVec3f(normalize(sunpos));
+    SGVec3d world_up = hlOr.backTransform(-SGVec3d::e3());
+    SGVec3d nsun = normalize(sunpos);
     // cout << "nup = " << nup[0] << "," << nup[1] << ","
     //      << nup[2] << endl;
     // cout << "nsun = " << nsun[0] << "," << nsun[1] << ","
@@ -461,62 +435,11 @@ void FGLight::updateSunPos()
     SG_LOG( SG_EVENT, SG_DEBUG, "sun angle relative to current location = "
             << get_sun_angle() );
 
-    // calculate vector to sun's position on the earth's surface
-    SGVec3d rel_sunpos = sunpos - v->get_view_pos();
-    // vector in cartesian coordinates from current position to the
-    // postion on the earth's surface the sun is directly over
-    SGVec3f to_sun = toVec3f(rel_sunpos);
-    // printf( "Vector to sun = %.2f %.2f %.2f\n",
-    //         v->to_sun[0], v->to_sun[1], v->to_sun[2]);
-
-    // Given a vector from the view position to the point on the
-    // earth's surface the sun is directly over, map into onto the
-    // local plane representing "horizontal".
-
-    // surface direction to go to head towards sun
-    SGVec3f surface_to_sun;
-    SGVec3f view_pos = toVec3f(v->get_view_pos());
-    surface_to_sun = map_vec_onto_cur_surface_plane(world_up, view_pos, to_sun);
-    surface_to_sun = normalize(surface_to_sun);
-    // cout << "(sg) Surface direction to sun is "
-    //   << surface_to_sun[0] << ","
-    //   << surface_to_sun[1] << ","
-    //   << surface_to_sun[2] << endl;
-    // cout << "Should be close to zero = "
-    //   << sgScalarProductVec3(nup, surface_to_sun) << endl;
-
-    // calculate the angle between surface_to_sun and
-    // v->get_surface_east().  We do this so we can sort out the
-    // acos() ambiguity.  I wish I could think of a more efficient
-    // way. :-(
-    SGVec3f surface_east(toVec3f(hlOr.backTransform(SGVec3d::e2())));
-    float east_dot = dot( surface_to_sun, surface_east );
-    // cout << "  East dot product = " << east_dot << endl;
-
-    // calculate the angle between v->surface_to_sun and
-    // v->surface_south.  this is how much we have to rotate the sky
-    // for it to align with the sun
-    SGVec3f surface_south(toVec3f(hlOr.backTransform(-SGVec3d::e1())));
-    float dot_ = dot( surface_to_sun, surface_south );
-    // cout << "  Dot product = " << dot << endl;
-
-    if (dot_ > 1.0) {
-        SG_LOG( SG_ASTRO, SG_INFO,
-                "Dot product  = " << dot_ << " is greater than 1.0" );
-        dot_ = 1.0;
-    }
-    else if (dot_ < -1.0) {
-         SG_LOG( SG_ASTRO, SG_INFO,
-                 "Dot product  = " << dot_ << " is less than -1.0" );
-         dot_ = -1.0;
-     }
-
-    if ( east_dot >= 0 ) {
-        set_sun_rotation( acos(dot_) );
-    } else {
-        set_sun_rotation( -acos(dot_) );
-    }
+    // Get direction to the sun in the local frame.
+    SGVec3d local_sun_vec = hlOr.transform(nsun);
+    // Angle from south. XXX Is this correct in the southern hemisphere?
+    double angle = atan2(local_sun_vec.x(), -local_sun_vec.y());
+    set_sun_rotation(angle);
     // cout << "  Sky needs to rotate = " << angle << " rads = "
     //      << angle * SGD_RADIANS_TO_DEGREES << " degrees." << endl;
-
 }

diff --git a/src/Time/sunsolver.cxx b/src/Time/sunsolver.cxx
index a83422759ef9bfcbe8252d76bc88f2b8ec6eb9a8..b903f5e568f73244fe868cf83c7dd78f9378bf50 100644 (file)
--- a/src/Time/sunsolver.cxx
+++ b/src/Time/sunsolver.cxx
@@ -58,16 +58,14 @@ void fgSunPositionGST(double gst, double *lon, double *lat) {
 
     SGPropertyNode* sun = fgGetNode("/ephemeris/sun");
     assert(sun);
-    double beta = sun->getDoubleValue("lat-deg");
-    // double r = globals->get_ephem()->get_sun()->getDistance();
     double xs = sun->getDoubleValue("xs");
     double ys = sun->getDoubleValue("ys");
     double ye = sun->getDoubleValue("ye");
     double ze = sun->getDoubleValue("ze");
-    alpha = atan2(ys - tan(beta)*ze/ys, xs);
-    delta = asin(sin(beta)*ye/ys + cos(beta)*ze);
+    double ra = atan2(ye, xs);
+    double dec = atan2(ze, sqrt(xs * xs + ye * ye));
 
-    tmp = alpha - (SGD_2PI/24)*gst;
+    tmp = ra - (SGD_2PI/24)*gst;
     if (tmp < -SGD_PI) {
         do tmp += SGD_2PI;
         while (tmp < -SGD_PI);
@@ -77,7 +75,7 @@ void fgSunPositionGST(double gst, double *lon, double *lat) {
     }
 
     *lon = tmp;
-    *lat = delta;
+    *lat = dec;
 }
 
 static double sun_angle( const SGTime &t, const SGVec3d& world_up,
@@ -85,17 +83,18 @@ static double sun_angle( const SGTime &t, const SGVec3d& world_up,
     SG_LOG( SG_EVENT, SG_DEBUG, "  Updating Sun position" );
     SG_LOG( SG_EVENT, SG_DEBUG, "  Gst = " << t.getGst() );
 
-    double sun_lon, sun_gd_lat;
-    fgSunPositionGST( t.getGst(), &sun_lon, &sun_gd_lat );
-    SGVec3d sunpos = SGVec3d::fromGeod(SGGeod::fromRad(sun_lon, sun_gd_lat));
+    double sun_lon, sun_gc_lat;
+    fgSunPositionGST( t.getGst(), &sun_lon, &sun_gc_lat );
+    SGVec3d sunpos = SGVec3d::fromGeoc(SGGeoc::fromRadM(sun_lon, sun_gc_lat,
+                                                        SGGeodesy::EQURAD));
 
     SG_LOG( SG_EVENT, SG_DEBUG, "    t.cur_time = " << t.get_cur_time() );
     SG_LOG( SG_EVENT, SG_DEBUG, 
-           "    Sun Geodetic lat = " << sun_gd_lat );
+           "    Sun Geocentric lat = " << sun_gc_lat );
 
     // calculate the sun's relative angle to local up
-    SGVec3f nup = normalize(toVec3f(world_up));
-    SGVec3f nsun = normalize(toVec3f(sunpos));
+    SGVec3d nup = normalize(world_up);
+    SGVec3d nsun = normalize(sunpos);
     // cout << "nup = " << nup[0] << "," << nup[1] << "," 
     //      << nup[2] << endl;
     // cout << "nsun = " << nsun[0] << "," << nsun[1] << ","