From: Tim Moore Date: Sat, 4 Dec 2010 15:25:39 +0000 (+0100) Subject: Fix bug 182 - make sun direction agree with position of sun from the ephemeris X-Git-Url: https://git.mxchange.org/?a=commitdiff_plain;h=79a2173ef9d6d04af8815d658c0a4daedf553ac4;p=flightgear.git 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 24d1c7caf..348bd524d 100644 --- 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 vec) / (normal 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 a83422759..b903f5e56 100644 --- 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] << ","