2 * sunsolver.cxx - given a location on earth and a time of day/date,
3 * find the number of seconds to various sun positions.
5 * Written by Curtis Olson, started September 2003.
7 * Copyright (C) 2003 Curtis L. Olson - http://www.flightgear.org/~curt
9 * This program is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU General Public License as
11 * published by the Free Software Foundation; either version 2 of the
12 * License, or (at your option) any later version.
14 * This program is distributed in the hope that it will be useful, but
15 * WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * General Public License for more details.
19 * You should have received a copy of the GNU General Public License
20 * along with this program; if not, write to the Free Software
21 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
34 #include <simgear/math/SGMath.hxx>
35 #include <simgear/timing/sg_time.hxx>
37 #include <Main/globals.hxx>
38 #include <Main/fg_props.hxx>
40 #include "sunsolver.hxx"
43 static const time_t day_secs = 86400;
44 static const time_t half_day_secs = day_secs / 2;
45 static const time_t step_secs = 60;
47 /* given a particular time expressed in side real time at prime
48 * meridian (GST), compute position on the earth (lat, lon) such that
49 * sun is directly overhead. (lat, lon are reported in radians */
51 void fgSunPositionGST(double gst, double *lon, double *lat) {
52 /* time_t ssue; seconds since unix epoch */
53 /* double *lat; (return) latitude */
54 /* double *lon; (return) longitude */
58 SGPropertyNode* sun = fgGetNode("/ephemeris/sun");
60 double xs = sun->getDoubleValue("xs");
61 //double ys = sun->getDoubleValue("ys");
62 double ye = sun->getDoubleValue("ye");
63 double ze = sun->getDoubleValue("ze");
64 double ra = atan2(ye, xs);
65 double dec = atan2(ze, sqrt(xs * xs + ye * ye));
67 tmp = ra - (SGD_2PI/24)*gst;
69 double signnedPI = (tmp < 0.0) ? -SGD_PI : SGD_PI;
70 tmp = fmod(tmp+signnedPI, SGD_2PI) - signnedPI;
76 static double sun_angle( const SGTime &t, const SGVec3d& world_up,
77 double lon_rad, double lat_rad ) {
78 SG_LOG( SG_EVENT, SG_DEBUG, " Updating Sun position" );
79 SG_LOG( SG_EVENT, SG_DEBUG, " Gst = " << t.getGst() );
81 double sun_lon, sun_gc_lat;
82 fgSunPositionGST( t.getGst(), &sun_lon, &sun_gc_lat );
83 SGVec3d sunpos = SGVec3d::fromGeoc(SGGeoc::fromRadM(sun_lon, sun_gc_lat,
86 SG_LOG( SG_EVENT, SG_DEBUG, " t.cur_time = " << t.get_cur_time() );
87 SG_LOG( SG_EVENT, SG_DEBUG,
88 " Sun Geocentric lat = " << sun_gc_lat );
90 // calculate the sun's relative angle to local up
91 SGVec3d nup = normalize(world_up);
92 SGVec3d nsun = normalize(sunpos);
93 // cout << "nup = " << nup[0] << "," << nup[1] << ","
95 // cout << "nsun = " << nsun[0] << "," << nsun[1] << ","
96 // << nsun[2] << endl;
98 double sun_angle = acos( dot( nup, nsun ) );
100 double signnedPI = (sun_angle < 0.0) ? -SGD_PI : SGD_PI;
101 sun_angle = fmod(sun_angle+signnedPI, SGD_2PI) - signnedPI;
103 double sun_angle_deg = sun_angle * SG_RADIANS_TO_DEGREES;
104 SG_LOG( SG_EVENT, SG_DEBUG, "sun angle relative to current location = "
107 return sun_angle_deg;
112 * Given the current unix time in seconds, calculate seconds to the
113 * specified sun angle (relative to straight up.) Also specify if we
114 * want the angle while the sun is ascending or descending. For
115 * instance noon is when the sun angle is 0 (or the closest it can
116 * get.) Dusk is when the sun angle is 90 and descending. Dawn is
117 * when the sun angle is 90 and ascending.
119 time_t fgTimeSecondsUntilSunAngle( time_t cur_time,
122 double target_angle_deg,
125 // cout << "location = " << lon_rad * SG_RADIANS_TO_DEGREES << ", "
126 // << lat_rad * SG_RADIANS_TO_DEGREES << endl;
127 SGVec3d world_up = SGVec3d::fromGeod(SGGeod::fromRad(lon_rad, lat_rad));
128 SGTime t = SGTime( lon_rad, lat_rad, "", 0 );
130 double best_diff = 180.0;
131 double last_angle = -99999.0;
132 time_t best_time = cur_time;
134 for ( time_t secs = cur_time - half_day_secs;
135 secs < cur_time + half_day_secs;
138 t.update( lon_rad, lat_rad, secs, 0 );
139 double angle_deg = sun_angle( t, world_up, lon_rad, lat_rad );
140 double diff = fabs( angle_deg - target_angle_deg );
141 if ( diff < best_diff ) {
142 if ( last_angle <= 180.0 && ascending
143 && ( last_angle > angle_deg ) ) {
144 // cout << "best angle = " << angle << " offset = "
145 // << secs - cur_time << endl;
148 } else if ( last_angle <= 180.0 && !ascending
149 && ( last_angle < angle_deg ) ) {
150 // cout << "best angle = " << angle << " offset = "
151 // << secs - cur_time << endl;
157 last_angle = angle_deg;
160 return best_time - cur_time;