X-Git-Url: https://git.mxchange.org/?a=blobdiff_plain;f=src%2FTime%2Fsunsolver.cxx;h=b37112816282f1e061219139f600a770087b47ad;hb=da73dd97d3d0e36b6078728ac39f0b98cae46ff7;hp=c04e366dd2a70286c741df704033301f5f898ad8;hpb=9bb782ce720071100e995b017f7a0051a2530568;p=flightgear.git diff --git a/src/Time/sunsolver.cxx b/src/Time/sunsolver.cxx index c04e366dd..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,59 +18,89 @@ * * 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" -const time_t day_secs = 86400; +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 */ - SG_LOG( SG_EVENT, SG_DEBUG, " Updating Sun position" ); - SG_LOG( SG_EVENT, SG_DEBUG, " Gst = " << t.getGst() ); +void fgSunPositionGST(double gst, double *lon, double *lat) { + /* time_t ssue; seconds since unix epoch */ + /* double *lat; (return) latitude */ + /* double *lon; (return) longitude */ - double sun_lon, sun_gd_lat, sun_gc_lat, sl_radius; - fgSunPositionGST( t.getGst(), &sun_lon, &sun_gd_lat ); + double tmp; - sgGeodToGeoc(sun_gd_lat, 0.0, &sl_radius, &sun_gc_lat); + 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() ); - 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 ); @@ -79,147 +109,44 @@ static double sun_angle( const SGTime &t, sgVec3 world_up, /** - * Given the current unix time in seconds, calculate seconds to noon + * 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 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; secs < cur_time + day_secs; secs += 300 ) { - 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; - } - } - - if ( best_time > day_secs / 2 ) { - best_time -= day_secs; - } - - return best_time - cur_time; -} - - -/** - * Given the current unix time in seconds, calculate seconds to midnight - */ -time_t fgTimeSecondsUntilMidnight( time_t cur_time, +time_t fgTimeSecondsUntilSunAngle( 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_angle = 0.0; - time_t best_time = cur_time; - - for ( time_t secs = cur_time; secs < cur_time + day_secs; secs += 300 ) { - 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; - } - } - - if ( best_time > day_secs / 2 ) { - best_time -= day_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 ) -{ - // 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 best_diff = 180.0; double last_angle = -99999.0; time_t best_time = cur_time; - for ( time_t secs = cur_time; secs < cur_time + day_secs; secs += 300 ) { + 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 ); + 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; - } - - if ( best_time > day_secs / 2 ) { - best_time -= day_secs; - } - - 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; secs < cur_time + day_secs; secs += 300 ) { - 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; @@ -227,11 +154,7 @@ time_t fgTimeSecondsUntilDawn( time_t cur_time, } } - last_angle = angle; - } - - if ( best_time > day_secs / 2 ) { - best_time -= day_secs; + last_angle = angle_deg; } return best_time - cur_time;