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., 675 Mass Ave, Cambridge, MA 02139, USA.
  22  *
  23  * $Id$
  24  */
  25
  26 #ifdef SG_HAVE_STD_INCLUDES
  27 #  include <cmath>
  28 // #  include <cstdio>
  29 #  include <ctime>
  30 #  ifdef macintosh
  31      SG_USING_STD(time_t);
  32 #  endif
  33 #else
  34 #  include <math.h>
  35 // #  include <stdio.h>
  36 #  include <time.h>
  37 #endif
  38
  39 #include <simgear/math/point3d.hxx>
  40 #include <simgear/math/sg_geodesy.hxx>
  41 #include <simgear/ephemeris/ephemeris.hxx>
  42 #include <simgear/timing/sg_time.hxx>
  43
  44 #include <Main/globals.hxx>
  45
  46 #include "tmp.hxx"
  47 #include "sunsolver.hxx"
  48
  49
  50 static const time_t day_secs = 86400;
  51 static const time_t half_day_secs = day_secs / 2;
  52 static const time_t step_secs = 60;
  53
  54 /* given a particular time expressed in side real time at prime
  55  * meridian (GST), compute position on the earth (lat, lon) such that
  56  * sun is directly overhead.  (lat, lon are reported in radians */
  57
  58 void fgSunPositionGST(double gst, double *lon, double *lat) {
  59     /* time_t  ssue;           seconds since unix epoch */
  60     /* double *lat;            (return) latitude        */
  61     /* double *lon;            (return) longitude       */
  62
  63     /* double lambda; */
  64     double alpha, delta;
  65     double tmp;
  66
  67     double beta = globals->get_ephem()->get_sun()->getLat();
  68     double r = globals->get_ephem()->get_sun()->getDistance();
  69     double xs = globals->get_ephem()->get_sun()->getxs();
  70     double ys = globals->get_ephem()->get_sun()->getys();
  71     double ye = globals->get_ephem()->get_sun()->getye();
  72     double ze = globals->get_ephem()->get_sun()->getze();
  73     alpha = atan2(ys - tan(beta)*ze/ys, xs);
  74     delta = asin(sin(beta)*ye/ys + cos(beta)*ze);
  75
  76     // tmp = alpha - (SGD_2PI/24)*GST(ssue);
  77     tmp = alpha - (SGD_2PI/24)*gst;
  78     if (tmp < -SGD_PI) {
  79         do tmp += SGD_2PI;
  80         while (tmp < -SGD_PI);
  81     } else if (tmp > SGD_PI) {
  82         do tmp -= SGD_2PI;
  83         while (tmp < -SGD_PI);
  84     }
  85
  86     *lon = tmp;
  87     *lat = delta;
  88 }
  89
  90 static double sun_angle( const SGTime &t, sgVec3 world_up,
  91                          double lon_rad, double lat_rad ) {
  92     sgVec3 nup, nsun;
  93     Point3D p, rel_sunpos;
  94
  95     SG_LOG( SG_EVENT, SG_DEBUG, "  Updating Sun position" );
  96     SG_LOG( SG_EVENT, SG_DEBUG, "  Gst = " << t.getGst() );
  97
  98     double sun_lon, sun_gd_lat;
  99     fgSunPositionGST( t.getGst(), &sun_lon, &sun_gd_lat );
 100     Point3D sunpos = sgGeodToCart(Point3D(sun_lon, sun_gd_lat, 0));
 101
 102     SG_LOG( SG_EVENT, SG_DEBUG, "    t.cur_time = " << t.get_cur_time() );
 103     SG_LOG( SG_EVENT, SG_DEBUG,
 104             "    Sun Geodetic lat = " << sun_gd_lat );
 105
 106     // calculate the sun's relative angle to local up
 107     sgCopyVec3( nup, world_up );
 108     sgSetVec3( nsun, sunpos.x(), sunpos.y(), sunpos.z() );
 109     sgNormalizeVec3(nup);
 110     sgNormalizeVec3(nsun);
 111     // cout << "nup = " << nup[0] << "," << nup[1] << ","
 112     //      << nup[2] << endl;
 113     // cout << "nsun = " << nsun[0] << "," << nsun[1] << ","
 114     //      << nsun[2] << endl;
 115
 116     double sun_angle = acos( sgScalarProductVec3 ( nup, nsun ) );
 117     double sun_angle_deg = sun_angle * SG_RADIANS_TO_DEGREES;
 118     while ( sun_angle_deg < -180 ) { sun_angle += 360; }
 119     SG_LOG( SG_EVENT, SG_DEBUG, "sun angle relative to current location = "
 120             << sun_angle_deg );
 121
 122     return sun_angle_deg;
 123 }
 124
 125
 126 /**
 127  * Given the current unix time in seconds, calculate seconds to the
 128  * specified sun angle (relative to straight up.)  Also specify if we
 129  * want the angle while the sun is ascending or descending.  For
 130  * instance noon is when the sun angle is 0 (or the closest it can
 131  * get.)  Dusk is when the sun angle is 90 and descending.  Dawn is
 132  * when the sun angle is 90 and ascending.
 133  */
 134 time_t fgTimeSecondsUntilSunAngle( time_t cur_time,
 135                                    double lon_rad,
 136                                    double lat_rad,
 137                                    double target_angle_deg,
 138                                    bool ascending )
 139 {
 140     // cout << "location = " << lon_rad * SG_RADIANS_TO_DEGREES << ", "
 141     //      << lat_rad * SG_RADIANS_TO_DEGREES << endl;
 142     Point3D geod( lon_rad, lat_rad, 0 );
 143     Point3D tmp = sgGeodToCart( geod );
 144     sgVec3 world_up;
 145     sgSetVec3( world_up, tmp.x(), tmp.y(), tmp.z() );
 146     SGTime t = SGTime( lon_rad, lat_rad, "", 0 );
 147
 148     double best_diff = 180.0;
 149     double last_angle = -99999.0;
 150     time_t best_time = cur_time;
 151
 152     for ( time_t secs = cur_time - half_day_secs;
 153           secs < cur_time + half_day_secs;
 154           secs += step_secs )
 155     {
 156         t.update( lon_rad, lat_rad, secs, 0 );
 157         double angle_deg = sun_angle( t, world_up, lon_rad, lat_rad );
 158         double diff = fabs( angle_deg - target_angle_deg );
 159         if ( diff < best_diff ) {
 160             if ( last_angle <= 180.0 && ascending
 161                  && ( last_angle > angle_deg ) ) {
 162                 // cout << "best angle = " << angle << " offset = "
 163                 //      << secs - cur_time << endl;
 164                 best_diff = diff;
 165                 best_time = secs;
 166             } else if ( last_angle <= 180.0 && !ascending
 167                         && ( last_angle < angle_deg ) ) {
 168                 // cout << "best angle = " << angle << " offset = "
 169                 //      << secs - cur_time << endl;
 170                 best_diff = diff;
 171                 best_time = secs;
 172             }
 173         }
 174
 175         last_angle = angle_deg;
 176     }
 177
 178     return best_time - cur_time;
 179 }