X-Git-Url: https://git.mxchange.org/?a=blobdiff_plain;f=src%2FTime%2Fsunsolver.cxx;h=b37112816282f1e061219139f600a770087b47ad;hb=da73dd97d3d0e36b6078728ac39f0b98cae46ff7;hp=c4e48c43403e21fe376b737f3ff2ed3c5958f56c;hpb=b10715833c0749343549d7169e3a7b7c6422d56e;p=flightgear.git diff --git a/src/Time/sunsolver.cxx b/src/Time/sunsolver.cxx index c4e48c434..b37112816 100644 --- a/src/Time/sunsolver.cxx +++ b/src/Time/sunsolver.cxx @@ -4,7 +4,7 @@ * * Written by Curtis Olson, started September 2003. * - * Copyright (C) 2003 Curtis L. Olson - curt@flightgear.org + * Copyright (C) 2003 Curtis L. Olson - http://www.flightgear.org/~curt * * This program is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License as @@ -18,19 +18,24 @@ * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software - * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. * * $Id$ */ +#ifdef HAVE_CONFIG_H +# include +#endif -#include -#include +#include +#include +#include + +#include #include #include
- -#include "sunpos.hxx" +#include
#include "sunsolver.hxx" @@ -39,40 +44,63 @@ static const time_t day_secs = 86400; static const time_t half_day_secs = day_secs / 2; static const time_t step_secs = 60; -static double sun_angle( const SGTime &t, sgVec3 world_up, - double lon_rad, double lat_rad ) { - sgVec3 nup, nsun; - Point3D p, rel_sunpos; +/* given a particular time expressed in side real time at prime + * meridian (GST), compute position on the earth (lat, lon) such that + * sun is directly overhead. (lat, lon are reported in radians */ + +void fgSunPositionGST(double gst, double *lon, double *lat) { + /* time_t ssue; seconds since unix epoch */ + /* double *lat; (return) latitude */ + /* double *lon; (return) longitude */ + + double tmp; + + SGPropertyNode* sun = fgGetNode("/ephemeris/sun"); + assert(sun); + double xs = sun->getDoubleValue("xs"); + //double ys = sun->getDoubleValue("ys"); + double ye = sun->getDoubleValue("ye"); + double ze = sun->getDoubleValue("ze"); + double ra = atan2(ye, xs); + double dec = atan2(ze, sqrt(xs * xs + ye * ye)); + + tmp = ra - (SGD_2PI/24)*gst; + + double signedPI = (tmp < 0.0) ? -SGD_PI : SGD_PI; + tmp = fmod(tmp+signedPI, SGD_2PI) - signedPI; + + *lon = tmp; + *lat = dec; +} +static double sun_angle( const SGTime &t, const SGVec3d& world_up, + double lon_rad, double lat_rad ) { SG_LOG( SG_EVENT, SG_DEBUG, " Updating Sun position" ); SG_LOG( SG_EVENT, SG_DEBUG, " Gst = " << t.getGst() ); - double sun_lon, sun_gd_lat, sun_gc_lat, sl_radius; - fgSunPositionGST( t.getGst(), &sun_lon, &sun_gd_lat ); - - sgGeodToGeoc(sun_gd_lat, 0.0, &sl_radius, &sun_gc_lat); - - p = Point3D( sun_lon, sun_gc_lat, sl_radius ); - Point3D sunpos = sgPolarToCart3d(p); + 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 - << " Geocentric lat = " << sun_gc_lat ); + " Sun Geocentric lat = " << sun_gc_lat ); // calculate the sun's relative angle to local up - sgCopyVec3( nup, world_up ); - sgSetVec3( nsun, sunpos.x(), sunpos.y(), sunpos.z() ); - sgNormalizeVec3(nup); - sgNormalizeVec3(nsun); + SGVec3d nup = normalize(world_up); + SGVec3d nsun = normalize(sunpos); // cout << "nup = " << nup[0] << "," << nup[1] << "," // << nup[2] << endl; // cout << "nsun = " << nsun[0] << "," << nsun[1] << "," // << nsun[2] << endl; - double sun_angle = acos( sgScalarProductVec3 ( nup, nsun ) ); + double sun_angle = acos( dot( nup, nsun ) ); + + double signedPI = (sun_angle < 0.0) ? -SGD_PI : SGD_PI; + sun_angle = fmod(sun_angle+signedPI, SGD_2PI) - signedPI; + double sun_angle_deg = sun_angle * SG_RADIANS_TO_DEGREES; - while ( sun_angle_deg < -180 ) { sun_angle += 360; } SG_LOG( SG_EVENT, SG_DEBUG, "sun angle relative to current location = " << sun_angle_deg ); @@ -81,93 +109,25 @@ static double sun_angle( const SGTime &t, sgVec3 world_up, /** - * Given the current unix time in seconds, calculate seconds to noon - */ -time_t fgTimeSecondsUntilNoon( time_t cur_time, - double lon_rad, - double lat_rad ) -{ - // cout << "location = " << lon_rad * SG_RADIANS_TO_DEGREES << ", " - // << lat_rad * SG_RADIANS_TO_DEGREES << endl; - Point3D geod( lon_rad, lat_rad, 0 ); - Point3D tmp = sgGeodToCart( geod ); - sgVec3 world_up; - sgSetVec3( world_up, tmp.x(), tmp.y(), tmp.z() ); - SGTime t = SGTime( lon_rad, lat_rad, "", 0 ); - - double best_angle = 180.0; - time_t best_time = cur_time; - - for ( time_t secs = cur_time - half_day_secs; - secs < cur_time + half_day_secs; - secs += step_secs ) - { - t.update( lon_rad, lat_rad, secs, 0 ); - double angle = sun_angle( t, world_up, lon_rad, lat_rad ); - if ( angle < best_angle ) { - // cout << "best angle = " << angle << " offset = " - // << secs - cur_time << endl; - best_angle = angle; - best_time = secs; - } - } - - return best_time - cur_time; -} - - -/** - * Given the current unix time in seconds, calculate seconds to midnight + * Given the current unix time in seconds, calculate seconds to the + * specified sun angle (relative to straight up.) Also specify if we + * want the angle while the sun is ascending or descending. For + * instance noon is when the sun angle is 0 (or the closest it can + * get.) Dusk is when the sun angle is 90 and descending. Dawn is + * when the sun angle is 90 and ascending. */ -time_t fgTimeSecondsUntilMidnight( time_t cur_time, +time_t fgTimeSecondsUntilSunAngle( time_t cur_time, double lon_rad, - double lat_rad ) -{ - // cout << "location = " << lon_rad * SG_RADIANS_TO_DEGREES << ", " - // << lat_rad * SG_RADIANS_TO_DEGREES << endl; - Point3D geod( lon_rad, lat_rad, 0 ); - Point3D tmp = sgGeodToCart( geod ); - sgVec3 world_up; - sgSetVec3( world_up, tmp.x(), tmp.y(), tmp.z() ); - SGTime t = SGTime( lon_rad, lat_rad, "", 0 ); - - double best_angle = 0.0; - time_t best_time = cur_time; - - for ( time_t secs = cur_time - half_day_secs; - secs < cur_time + half_day_secs; - secs += step_secs ) - { - t.update( lon_rad, lat_rad, secs, 0 ); - double angle = sun_angle( t, world_up, lon_rad, lat_rad ); - if ( angle > best_angle ) { - // cout << "best angle = " << angle << " offset = " - // << secs - cur_time << endl; - best_angle = angle; - best_time = secs; - } - } - - return best_time - cur_time; -} - - -/** - * Given the current unix time in seconds, calculate seconds to dusk - */ -time_t fgTimeSecondsUntilDusk( time_t cur_time, - double lon_rad, - double lat_rad ) + double lat_rad, + double target_angle_deg, + bool ascending ) { // cout << "location = " << lon_rad * SG_RADIANS_TO_DEGREES << ", " // << lat_rad * SG_RADIANS_TO_DEGREES << endl; - Point3D geod( lon_rad, lat_rad, 0 ); - Point3D tmp = sgGeodToCart( geod ); - sgVec3 world_up; - sgSetVec3( world_up, tmp.x(), tmp.y(), tmp.z() ); + SGVec3d world_up = SGVec3d::fromGeod(SGGeod::fromRad(lon_rad, lat_rad)); SGTime t = SGTime( lon_rad, lat_rad, "", 0 ); - double best_diff = 90.0; + double best_diff = 180.0; double last_angle = -99999.0; time_t best_time = cur_time; @@ -176,52 +136,17 @@ time_t fgTimeSecondsUntilDusk( time_t cur_time, secs += step_secs ) { t.update( lon_rad, lat_rad, secs, 0 ); - double angle = sun_angle( t, world_up, lon_rad, lat_rad ); - double diff = fabs( angle - 90.0 ); + double angle_deg = sun_angle( t, world_up, lon_rad, lat_rad ); + double diff = fabs( angle_deg - target_angle_deg ); if ( diff < best_diff ) { - if ( last_angle <= 180.0 && ( last_angle < angle ) ) { + if ( last_angle <= 180.0 && ascending + && ( last_angle > angle_deg ) ) { // cout << "best angle = " << angle << " offset = " // << secs - cur_time << endl; best_diff = diff; best_time = secs; - } - } - - last_angle = angle; - } - - return best_time - cur_time; -} - - -/** - * Given the current unix time in seconds, calculate seconds to dawn - */ -time_t fgTimeSecondsUntilDawn( time_t cur_time, - double lon_rad, - double lat_rad ) -{ - // cout << "location = " << lon_rad * SG_RADIANS_TO_DEGREES << ", " - // << lat_rad * SG_RADIANS_TO_DEGREES << endl; - Point3D geod( lon_rad, lat_rad, 0 ); - Point3D tmp = sgGeodToCart( geod ); - sgVec3 world_up; - sgSetVec3( world_up, tmp.x(), tmp.y(), tmp.z() ); - SGTime t = SGTime( lon_rad, lat_rad, "", 0 ); - - double best_diff = 90.0; - double last_angle = -99999.0; - time_t best_time = cur_time; - - for ( time_t secs = cur_time - half_day_secs; - secs < cur_time + half_day_secs; - secs += step_secs ) - { - t.update( lon_rad, lat_rad, secs, 0 ); - double angle = sun_angle( t, world_up, lon_rad, lat_rad ); - double diff = fabs( angle - 90.0 ); - if ( diff < best_diff ) { - if ( last_angle <= 180.0 && ( last_angle > angle ) ) { + } else if ( last_angle <= 180.0 && !ascending + && ( last_angle < angle_deg ) ) { // cout << "best angle = " << angle << " offset = " // << secs - cur_time << endl; best_diff = diff; @@ -229,7 +154,7 @@ time_t fgTimeSecondsUntilDawn( time_t cur_time, } } - last_angle = angle; + last_angle = angle_deg; } return best_time - cur_time;