src/Time/sunsolver.cxx

   1 /*
   2  * sunsolver.cxx - given a location on earth and a time of day/date,
   3  *                 find the number of seconds to various sun positions.
   4  *
   5  * Written by Curtis Olson, started September 2003.
   6  *
   7  * Copyright (C) 2003  Curtis L. Olson  - http://www.flightgear.org/~curt
   8  *
   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.
  13  *
  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.
  18  *
  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.
  22  *
  23  * $Id$
  24  */
  25
  26 #ifdef HAVE_CONFIG_H
  27 #      include <config.h>
  28 #endif
  29
  30 #include <cmath>
  31 #include <ctime>
  32 #include <cassert>
  33
  34 #include <simgear/math/SGMath.hxx>
  35 #include <simgear/timing/sg_time.hxx>
  36
  37 #include <Main/globals.hxx>
  38 #include <Main/fg_props.hxx>
  39
  40 #include "sunsolver.hxx"
  41
  42
  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;
  46
  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 */
  50
  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       */
  55
  56     double tmp;
  57
  58     SGPropertyNode* sun = fgGetNode("/ephemeris/sun");
  59     assert(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));
  66
  67     tmp = ra - (SGD_2PI/24)*gst;
  68
  69     double signedPI = (tmp < 0.0) ? -SGD_PI : SGD_PI;
  70     tmp = fmod(tmp+signedPI, SGD_2PI) - signedPI;
  71
  72     *lon = tmp;
  73     *lat = dec;
  74 }
  75
  76 static double sun_angle( const SGTime &t, const SGVec3d& world_up) {
  77     SG_LOG( SG_EVENT, SG_DEBUG, "  Updating Sun position" );
  78     SG_LOG( SG_EVENT, SG_DEBUG, "  Gst = " << t.getGst() );
  79
  80     double sun_lon, sun_gc_lat;
  81     fgSunPositionGST( t.getGst(), &sun_lon, &sun_gc_lat );
  82     SGVec3d sunpos = SGVec3d::fromGeoc(SGGeoc::fromRadM(sun_lon, sun_gc_lat,
  83                                                         SGGeodesy::EQURAD));
  84
  85     SG_LOG( SG_EVENT, SG_DEBUG, "    t.cur_time = " << t.get_cur_time() );
  86     SG_LOG( SG_EVENT, SG_DEBUG,
  87             "    Sun Geocentric lat = " << sun_gc_lat );
  88
  89     // calculate the sun's relative angle to local up
  90     SGVec3d nup = normalize(world_up);
  91     SGVec3d nsun = normalize(sunpos);
  92     // cout << "nup = " << nup[0] << "," << nup[1] << ","
  93     //      << nup[2] << endl;
  94     // cout << "nsun = " << nsun[0] << "," << nsun[1] << ","
  95     //      << nsun[2] << endl;
  96
  97     double sun_angle = acos( dot( nup, nsun ) );
  98
  99     double signedPI = (sun_angle < 0.0) ? -SGD_PI : SGD_PI;
 100     sun_angle = fmod(sun_angle+signedPI, SGD_2PI) - signedPI;
 101
 102     double sun_angle_deg = sun_angle * SG_RADIANS_TO_DEGREES;
 103     SG_LOG( SG_EVENT, SG_DEBUG, "sun angle relative to current location = "
 104             << sun_angle_deg );
 105
 106     return sun_angle_deg;
 107 }
 108
 109
 110 /**
 111  * Given the current unix time in seconds, calculate seconds to the
 112  * specified sun angle (relative to straight up.)  Also specify if we
 113  * want the angle while the sun is ascending or descending.  For
 114  * instance noon is when the sun angle is 0 (or the closest it can
 115  * get.)  Dusk is when the sun angle is 90 and descending.  Dawn is
 116  * when the sun angle is 90 and ascending.
 117  */
 118 time_t fgTimeSecondsUntilSunAngle( time_t cur_time,
 119                                    const SGGeod& loc,
 120                                    double target_angle_deg,
 121                                    bool ascending )
 122 {
 123     SGVec3d world_up = SGVec3d::fromGeod(loc);
 124     SGTime t = SGTime( loc, SGPath(), 0 );
 125
 126     double best_diff = 180.0;
 127     double last_angle = -99999.0;
 128     time_t best_time = cur_time;
 129
 130     for ( time_t secs = cur_time - half_day_secs;
 131           secs < cur_time + half_day_secs;
 132           secs += step_secs )
 133     {
 134         t.update( loc, secs, 0 );
 135         double angle_deg = sun_angle( t, world_up );
 136         double diff = fabs( angle_deg - target_angle_deg );
 137         if ( diff < best_diff ) {
 138             if ( last_angle <= 180.0 && ascending
 139                  && ( last_angle > angle_deg ) ) {
 140                 // cout << "best angle = " << angle << " offset = "
 141                 //      << secs - cur_time << endl;
 142                 best_diff = diff;
 143                 best_time = secs;
 144             } else if ( last_angle <= 180.0 && !ascending
 145                         && ( last_angle < angle_deg ) ) {
 146                 // cout << "best angle = " << angle << " offset = "
 147                 //      << secs - cur_time << endl;
 148                 best_diff = diff;
 149                 best_time = secs;
 150             }
 151         }
 152
 153         last_angle = angle_deg;
 154     }
 155
 156     return best_time - cur_time;
 157 }