fly/flightgear
1
0
Fork 0
Code Activity

Merge branch 'next' of gitorious.org:fg/flightgear into next

Browse source
This commit is contained in:
Dave Luff 2010-12-05 10:48:02 +00:00
commit ca7f9dec76
2 changed files with 30 additions and 108 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

115
src/Time/light.cxx
View file

@ -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;
}

23
src/Time/sunsolver.cxx
View file

@ -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] << ","