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 const time_t day_secs = 86400;
  39
  40 static double sun_angle( const SGTime &t, sgVec3 world_up,
  41                          double lon_rad, double lat_rad ) {
  42     sgVec3 nup, nsun;
  43     Point3D p, rel_sunpos;
  44
  45     SG_LOG( SG_EVENT, SG_DEBUG, "  Updating Sun position" );
  46     SG_LOG( SG_EVENT, SG_DEBUG, "  Gst = " << t.getGst() );
  47
  48     double sun_lon, sun_gd_lat, sun_gc_lat, sl_radius;
  49     fgSunPositionGST( t.getGst(), &sun_lon, &sun_gd_lat );
  50
  51     sgGeodToGeoc(sun_gd_lat, 0.0, &sl_radius, &sun_gc_lat);
  52
  53     p = Point3D( sun_lon, sun_gc_lat, sl_radius );
  54     Point3D sunpos = sgPolarToCart3d(p);
  55
  56     SG_LOG( SG_EVENT, SG_DEBUG, "    t.cur_time = " << t.get_cur_time() );
  57     SG_LOG( SG_EVENT, SG_DEBUG,
  58             "    Sun Geodetic lat = " << sun_gd_lat
  59             << " Geocentric lat = " << sun_gc_lat );
  60
  61     // calculate the sun's relative angle to local up
  62     sgCopyVec3( nup, world_up );
  63     sgSetVec3( nsun, sunpos.x(), sunpos.y(), sunpos.z() );
  64     sgNormalizeVec3(nup);
  65     sgNormalizeVec3(nsun);
  66     // cout << "nup = " << nup[0] << "," << nup[1] << ","
  67     //      << nup[2] << endl;
  68     // cout << "nsun = " << nsun[0] << "," << nsun[1] << ","
  69     //      << nsun[2] << endl;
  70
  71     double sun_angle = acos( sgScalarProductVec3 ( nup, nsun ) );
  72     double sun_angle_deg = sun_angle * SG_RADIANS_TO_DEGREES;
  73     while ( sun_angle_deg < -180 ) { sun_angle += 360; }
  74     SG_LOG( SG_EVENT, SG_DEBUG, "sun angle relative to current location = "
  75             << sun_angle_deg );
  76
  77     return sun_angle_deg;
  78 }
  79
  80
  81 /**
  82  * Given the current unix time in seconds, calculate seconds to noon
  83  */
  84 time_t fgTimeSecondsUntilNoon( time_t cur_time,
  85                                double lon_rad,
  86                                double lat_rad )
  87 {
  88     // cout << "location = " << lon_rad * SG_RADIANS_TO_DEGREES << ", "
  89     //      << lat_rad * SG_RADIANS_TO_DEGREES << endl;
  90     Point3D geod( lon_rad, lat_rad, 0 );
  91     Point3D tmp = sgGeodToCart( geod );
  92     sgVec3 world_up;
  93     sgSetVec3( world_up, tmp.x(), tmp.y(), tmp.z() );
  94     SGTime t = SGTime( lon_rad, lat_rad, "", 0 );
  95
  96     double best_angle = 180.0;
  97     time_t best_time = cur_time;
  98
  99     for ( time_t secs = cur_time; secs < cur_time + day_secs; secs += 300 ) {
 100         t.update( lon_rad, lat_rad, secs, 0 );
 101         double angle = sun_angle( t, world_up, lon_rad, lat_rad );
 102         if ( angle < best_angle ) {
 103             // cout << "best angle = " << angle << " offset = "
 104             // << secs - cur_time << endl;
 105             best_angle = angle;
 106             best_time = secs;
 107         }
 108     }
 109
 110     if ( best_time > day_secs / 2 ) {
 111         best_time -= day_secs;
 112     }
 113
 114     return best_time - cur_time;
 115 }
 116
 117
 118 /**
 119  * Given the current unix time in seconds, calculate seconds to midnight
 120  */
 121 time_t fgTimeSecondsUntilMidnight( time_t cur_time,
 122                                    double lon_rad,
 123                                    double lat_rad )
 124 {
 125     // cout << "location = " << lon_rad * SG_RADIANS_TO_DEGREES << ", "
 126     //      << lat_rad * SG_RADIANS_TO_DEGREES << endl;
 127     Point3D geod( lon_rad, lat_rad, 0 );
 128     Point3D tmp = sgGeodToCart( geod );
 129     sgVec3 world_up;
 130     sgSetVec3( world_up, tmp.x(), tmp.y(), tmp.z() );
 131     SGTime t = SGTime( lon_rad, lat_rad, "", 0 );
 132
 133     double best_angle = 0.0;
 134     time_t best_time = cur_time;
 135
 136     for ( time_t secs = cur_time; secs < cur_time + day_secs; secs += 300 ) {
 137         t.update( lon_rad, lat_rad, secs, 0 );
 138         double angle = sun_angle( t, world_up, lon_rad, lat_rad );
 139         if ( angle > best_angle ) {
 140             // cout << "best angle = " << angle << " offset = "
 141             //      << secs - cur_time << endl;
 142             best_angle = angle;
 143             best_time = secs;
 144         }
 145     }
 146
 147     if ( best_time > day_secs / 2 ) {
 148         best_time -= day_secs;
 149     }
 150
 151     return best_time - cur_time;
 152 }
 153
 154
 155 /**
 156  * Given the current unix time in seconds, calculate seconds to dusk
 157  */
 158 time_t fgTimeSecondsUntilDusk( time_t cur_time,
 159                                double lon_rad,
 160                                double lat_rad )
 161 {
 162     // cout << "location = " << lon_rad * SG_RADIANS_TO_DEGREES << ", "
 163     //      << lat_rad * SG_RADIANS_TO_DEGREES << endl;
 164     Point3D geod( lon_rad, lat_rad, 0 );
 165     Point3D tmp = sgGeodToCart( geod );
 166     sgVec3 world_up;
 167     sgSetVec3( world_up, tmp.x(), tmp.y(), tmp.z() );
 168     SGTime t = SGTime( lon_rad, lat_rad, "", 0 );
 169
 170     double best_diff = 90.0;
 171     double last_angle = -99999.0;
 172     time_t best_time = cur_time;
 173
 174     for ( time_t secs = cur_time; secs < cur_time + day_secs; secs += 300 ) {
 175         t.update( lon_rad, lat_rad, secs, 0 );
 176         double angle = sun_angle( t, world_up, lon_rad, lat_rad );
 177         double diff = fabs( angle - 90.0 );
 178         if ( diff < best_diff ) {
 179             if ( last_angle <= 180.0 && ( last_angle < angle ) ) {
 180                 // cout << "best angle = " << angle << " offset = "
 181                 //      << secs - cur_time << endl;
 182                 best_diff = diff;
 183                 best_time = secs;
 184             }
 185         }
 186
 187         last_angle = angle;
 188     }
 189
 190     if ( best_time > day_secs / 2 ) {
 191         best_time -= day_secs;
 192     }
 193
 194     return best_time - cur_time;
 195 }
 196
 197
 198 /**
 199  * Given the current unix time in seconds, calculate seconds to dawn
 200  */
 201 time_t fgTimeSecondsUntilDawn( time_t cur_time,
 202                                double lon_rad,
 203                                double lat_rad )
 204 {
 205     // cout << "location = " << lon_rad * SG_RADIANS_TO_DEGREES << ", "
 206     //      << lat_rad * SG_RADIANS_TO_DEGREES << endl;
 207     Point3D geod( lon_rad, lat_rad, 0 );
 208     Point3D tmp = sgGeodToCart( geod );
 209     sgVec3 world_up;
 210     sgSetVec3( world_up, tmp.x(), tmp.y(), tmp.z() );
 211     SGTime t = SGTime( lon_rad, lat_rad, "", 0 );
 212
 213     double best_diff = 90.0;
 214     double last_angle = -99999.0;
 215     time_t best_time = cur_time;
 216
 217     for ( time_t secs = cur_time; secs < cur_time + day_secs; secs += 300 ) {
 218         t.update( lon_rad, lat_rad, secs, 0 );
 219         double angle = sun_angle( t, world_up, lon_rad, lat_rad );
 220         double diff = fabs( angle - 90.0 );
 221         if ( diff < best_diff ) {
 222             if ( last_angle <= 180.0 && ( last_angle > angle ) ) {
 223                 // cout << "best angle = " << angle << " offset = "
 224                 //      << secs - cur_time << endl;
 225                 best_diff = diff;
 226                 best_time = secs;
 227             }
 228         }
 229
 230         last_angle = angle;
 231     }
 232
 233     if ( best_time > day_secs / 2 ) {
 234         best_time -= day_secs;
 235     }
 236
 237     return best_time - cur_time;
 238 }