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  - curt@flightgear.org
   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
  27 #include <simgear/math/point3d.hxx>
  28 #include <simgear/math/sg_geodesy.hxx>
  29 #include <simgear/timing/sg_time.hxx>
  30
  31 #include <Main/globals.hxx>
  32
  33 #include "sunpos.hxx"
  34
  35 #include "sunsolver.hxx"
  36
  37
  38 static const time_t day_secs = 86400;
  39 static const time_t half_day_secs = day_secs / 2;
  40 static const time_t step_secs = 60;
  41
  42 static double sun_angle( const SGTime &t, sgVec3 world_up,
  43                          double lon_rad, double lat_rad ) {
  44     sgVec3 nup, nsun;
  45     Point3D p, rel_sunpos;
  46
  47     SG_LOG( SG_EVENT, SG_DEBUG, "  Updating Sun position" );
  48     SG_LOG( SG_EVENT, SG_DEBUG, "  Gst = " << t.getGst() );
  49
  50     double sun_lon, sun_gd_lat, sun_gc_lat, sl_radius;
  51     fgSunPositionGST( t.getGst(), &sun_lon, &sun_gd_lat );
  52
  53     sgGeodToGeoc(sun_gd_lat, 0.0, &sl_radius, &sun_gc_lat);
  54
  55     p = Point3D( sun_lon, sun_gc_lat, sl_radius );
  56     Point3D sunpos = sgPolarToCart3d(p);
  57
  58     SG_LOG( SG_EVENT, SG_DEBUG, "    t.cur_time = " << t.get_cur_time() );
  59     SG_LOG( SG_EVENT, SG_DEBUG,
  60             "    Sun Geodetic lat = " << sun_gd_lat
  61             << " Geocentric lat = " << sun_gc_lat );
  62
  63     // calculate the sun's relative angle to local up
  64     sgCopyVec3( nup, world_up );
  65     sgSetVec3( nsun, sunpos.x(), sunpos.y(), sunpos.z() );
  66     sgNormalizeVec3(nup);
  67     sgNormalizeVec3(nsun);
  68     // cout << "nup = " << nup[0] << "," << nup[1] << ","
  69     //      << nup[2] << endl;
  70     // cout << "nsun = " << nsun[0] << "," << nsun[1] << ","
  71     //      << nsun[2] << endl;
  72
  73     double sun_angle = acos( sgScalarProductVec3 ( nup, nsun ) );
  74     double sun_angle_deg = sun_angle * SG_RADIANS_TO_DEGREES;
  75     while ( sun_angle_deg < -180 ) { sun_angle += 360; }
  76     SG_LOG( SG_EVENT, SG_DEBUG, "sun angle relative to current location = "
  77             << sun_angle_deg );
  78
  79     return sun_angle_deg;
  80 }
  81
  82
  83 /**
  84  * Given the current unix time in seconds, calculate seconds to the
  85  * specified sun angle (relative to straight up.)  Also specify if we
  86  * want the angle while the sun is ascending or descending.  For
  87  * instance noon is when the sun angle is 0 (or the closest it can
  88  * get.)  Dusk is when the sun angle is 90 and descending.  Dawn is
  89  * when the sun angle is 90 and ascending.
  90  */
  91 time_t fgTimeSecondsUntilSunAngle( time_t cur_time,
  92                                    double lon_rad,
  93                                    double lat_rad,
  94                                    double target_angle_deg,
  95                                    bool ascending )
  96 {
  97     // cout << "location = " << lon_rad * SG_RADIANS_TO_DEGREES << ", "
  98     //      << lat_rad * SG_RADIANS_TO_DEGREES << endl;
  99     Point3D geod( lon_rad, lat_rad, 0 );
 100     Point3D tmp = sgGeodToCart( geod );
 101     sgVec3 world_up;
 102     sgSetVec3( world_up, tmp.x(), tmp.y(), tmp.z() );
 103     SGTime t = SGTime( lon_rad, lat_rad, "", 0 );
 104
 105     double best_diff = 180.0;
 106     double last_angle = -99999.0;
 107     time_t best_time = cur_time;
 108
 109     for ( time_t secs = cur_time - half_day_secs;
 110           secs < cur_time + half_day_secs;
 111           secs += step_secs )
 112     {
 113         t.update( lon_rad, lat_rad, secs, 0 );
 114         double angle_deg = sun_angle( t, world_up, lon_rad, lat_rad );
 115         double diff = fabs( angle_deg - target_angle_deg );
 116         if ( diff < best_diff ) {
 117             if ( last_angle <= 180.0 && ascending
 118                  && ( last_angle > angle_deg ) ) {
 119                 // cout << "best angle = " << angle << " offset = "
 120                 //      << secs - cur_time << endl;
 121                 best_diff = diff;
 122                 best_time = secs;
 123             } else if ( last_angle <= 180.0 && !ascending
 124                         && ( last_angle < angle_deg ) ) {
 125                 // cout << "best angle = " << angle << " offset = "
 126                 //      << secs - cur_time << endl;
 127                 best_diff = diff;
 128                 best_time = secs;
 129             }
 130         }
 131
 132         last_angle = angle_deg;
 133     }
 134
 135     return best_time - cur_time;
 136 }