X-Git-Url: https://git.mxchange.org/?a=blobdiff_plain;f=src%2FTime%2Fsunsolver.cxx;h=b37112816282f1e061219139f600a770087b47ad;hb=da73dd97d3d0e36b6078728ac39f0b98cae46ff7;hp=384761c47b21eb60f7cda0f55bb01f686efd365a;hpb=d05121ef4689d2b50b3fe1848cbb0d1f5a1db877;p=flightgear.git

diff --git a/src/Time/sunsolver.cxx b/src/Time/sunsolver.cxx
index 384761c47..b37112816 100644
--- a/src/Time/sunsolver.cxx
+++ b/src/Time/sunsolver.cxx
@@ -18,19 +18,24 @@
  *
  * You should have received a copy of the GNU General Public License
  * along with this program; if not, write to the Free Software
- * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  */
 
+#ifdef HAVE_CONFIG_H
+#      include <config.h>
+#endif
 
-#include <simgear/math/point3d.hxx>
-#include <simgear/math/sg_geodesy.hxx>
+#include <cmath>
+#include <ctime>
+#include <cassert>
+
+#include <simgear/math/SGMath.hxx>
 #include <simgear/timing/sg_time.hxx>
 
 #include <Main/globals.hxx>
-
-#include "sunpos.hxx"
+#include <Main/fg_props.hxx>
 
 #include "sunsolver.hxx"
 
@@ -39,35 +44,63 @@ static const time_t day_secs = 86400;
 static const time_t half_day_secs = day_secs / 2;
 static const time_t step_secs = 60;
 
-static double sun_angle( const SGTime &t, sgVec3 world_up,
-                         double lon_rad, double lat_rad ) {
-    sgVec3 nup, nsun;
-    Point3D p, rel_sunpos;
+/* given a particular time expressed in side real time at prime
+ * meridian (GST), compute position on the earth (lat, lon) such that
+ * sun is directly overhead.  (lat, lon are reported in radians */
+
+void fgSunPositionGST(double gst, double *lon, double *lat) {
+    /* time_t  ssue;           seconds since unix epoch */
+    /* double *lat;            (return) latitude        */
+    /* double *lon;            (return) longitude       */
+
+    double tmp;
+
+    SGPropertyNode* sun = fgGetNode("/ephemeris/sun");
+    assert(sun);
+    double xs = sun->getDoubleValue("xs");
+    //double ys = sun->getDoubleValue("ys");
+    double ye = sun->getDoubleValue("ye");
+    double ze = sun->getDoubleValue("ze");
+    double ra = atan2(ye, xs);
+    double dec = atan2(ze, sqrt(xs * xs + ye * ye));
+
+    tmp = ra - (SGD_2PI/24)*gst;
+    
+    double signedPI = (tmp < 0.0) ? -SGD_PI : SGD_PI;
+    tmp = fmod(tmp+signedPI, SGD_2PI) - signedPI;
+
+    *lon = tmp;
+    *lat = dec;
+}
 
+static double sun_angle( const SGTime &t, const SGVec3d& world_up,
+                         double lon_rad, double lat_rad ) {
     SG_LOG( SG_EVENT, SG_DEBUG, "  Updating Sun position" );
     SG_LOG( SG_EVENT, SG_DEBUG, "  Gst = " << t.getGst() );
 
-    double sun_lon, sun_gd_lat;
-    fgSunPositionGST( t.getGst(), &sun_lon, &sun_gd_lat );
-    Point3D sunpos = sgGeodToCart(Point3D(sun_lon, sun_gd_lat, 0));
+    double sun_lon, sun_gc_lat;
+    fgSunPositionGST( t.getGst(), &sun_lon, &sun_gc_lat );
+    SGVec3d sunpos = SGVec3d::fromGeoc(SGGeoc::fromRadM(sun_lon, sun_gc_lat,
+                                                        SGGeodesy::EQURAD));
 
     SG_LOG( SG_EVENT, SG_DEBUG, "    t.cur_time = " << t.get_cur_time() );
     SG_LOG( SG_EVENT, SG_DEBUG, 
-	    "    Sun Geodetic lat = " << sun_gd_lat );
+	    "    Sun Geocentric lat = " << sun_gc_lat );
 
     // calculate the sun's relative angle to local up
-    sgCopyVec3( nup, world_up );
-    sgSetVec3( nsun, sunpos.x(), sunpos.y(), sunpos.z() );
-    sgNormalizeVec3(nup);
-    sgNormalizeVec3(nsun);
+    SGVec3d nup = normalize(world_up);
+    SGVec3d nsun = normalize(sunpos);
     // cout << "nup = " << nup[0] << "," << nup[1] << "," 
     //      << nup[2] << endl;
     // cout << "nsun = " << nsun[0] << "," << nsun[1] << "," 
     //      << nsun[2] << endl;
 
-    double sun_angle = acos( sgScalarProductVec3 ( nup, nsun ) );
+    double sun_angle = acos( dot( nup, nsun ) );
+
+    double signedPI = (sun_angle < 0.0) ? -SGD_PI : SGD_PI;
+    sun_angle = fmod(sun_angle+signedPI, SGD_2PI) - signedPI;
+
     double sun_angle_deg = sun_angle * SG_RADIANS_TO_DEGREES;
-    while ( sun_angle_deg < -180 ) { sun_angle += 360; }
     SG_LOG( SG_EVENT, SG_DEBUG, "sun angle relative to current location = "
 	    << sun_angle_deg );
 
@@ -91,10 +124,7 @@ time_t fgTimeSecondsUntilSunAngle( time_t cur_time,
 {
     // cout << "location = " << lon_rad * SG_RADIANS_TO_DEGREES << ", "
     //      << lat_rad * SG_RADIANS_TO_DEGREES << endl;
-    Point3D geod( lon_rad, lat_rad, 0 );
-    Point3D tmp = sgGeodToCart( geod );
-    sgVec3 world_up;
-    sgSetVec3( world_up, tmp.x(), tmp.y(), tmp.z() );
+    SGVec3d world_up = SGVec3d::fromGeod(SGGeod::fromRad(lon_rad, lat_rad));
     SGTime t = SGTime( lon_rad, lat_rad, "", 0 );
 
     double best_diff = 180.0;