]> git.mxchange.org Git - flightgear.git/blobdiff - src/Time/sunsolver.cxx
toggle fullscreen: also adapt GUI plane when resizing
[flightgear.git] / src / Time / sunsolver.cxx
index a83422759ef9bfcbe8252d76bc88f2b8ec6eb9a8..cc8a64e8f266fd2c58598dc60e0e0209fa556682 100644 (file)
@@ -53,57 +53,53 @@ void fgSunPositionGST(double gst, double *lon, double *lat) {
     /* double *lat;            (return) latitude        */
     /* double *lon;            (return) longitude       */
 
-    double alpha, delta;
     double tmp;
 
     SGPropertyNode* sun = fgGetNode("/ephemeris/sun");
     assert(sun);
-    double beta = sun->getDoubleValue("lat-deg");
-    // double r = globals->get_ephem()->get_sun()->getDistance();
     double xs = sun->getDoubleValue("xs");
-    double ys = sun->getDoubleValue("ys");
+    //double ys = sun->getDoubleValue("ys");
     double ye = sun->getDoubleValue("ye");
     double ze = sun->getDoubleValue("ze");
-    alpha = atan2(ys - tan(beta)*ze/ys, xs);
-    delta = asin(sin(beta)*ye/ys + cos(beta)*ze);
-
-    tmp = alpha - (SGD_2PI/24)*gst;
-    if (tmp < -SGD_PI) {
-        do tmp += SGD_2PI;
-        while (tmp < -SGD_PI);
-    } else if (tmp > SGD_PI) {
-        do tmp -= SGD_2PI;
-        while (tmp < -SGD_PI);
-    }
+    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 = delta;
+    *lat = dec;
 }
 
-static double sun_angle( const SGTime &t, const SGVec3d& world_up,
-                         double lon_rad, double lat_rad ) {
+static double sun_angle( const SGTime &t, const SGVec3d& world_up) {
     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 );
-    SGVec3d sunpos = SGVec3d::fromGeod(SGGeod::fromRad(sun_lon, sun_gd_lat));
+    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
-    SGVec3f nup = normalize(toVec3f(world_up));
-    SGVec3f nsun = normalize(toVec3f(sunpos));
+    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( 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 );
 
@@ -120,15 +116,12 @@ static double sun_angle( const SGTime &t, const SGVec3d& world_up,
  * when the sun angle is 90 and ascending.
  */
 time_t fgTimeSecondsUntilSunAngle( time_t cur_time,
-                                   double lon_rad,
-                                   double lat_rad,
+                                   const SGGeod& loc,
                                    double target_angle_deg,
                                    bool ascending )
 {
-    // cout << "location = " << lon_rad * SG_RADIANS_TO_DEGREES << ", "
-    //      << lat_rad * SG_RADIANS_TO_DEGREES << endl;
-    SGVec3d world_up = SGVec3d::fromGeod(SGGeod::fromRad(lon_rad, lat_rad));
-    SGTime t = SGTime( lon_rad, lat_rad, "", 0 );
+    SGVec3d world_up = SGVec3d::fromGeod(loc);
+    SGTime t = SGTime( loc, SGPath(), 0 );
 
     double best_diff = 180.0;
     double last_angle = -99999.0;
@@ -138,8 +131,8 @@ time_t fgTimeSecondsUntilSunAngle( time_t cur_time,
           secs < cur_time + half_day_secs;
           secs += step_secs )
     {
-        t.update( lon_rad, lat_rad, secs, 0 );
-        double angle_deg = sun_angle( t, world_up, lon_rad, lat_rad );
+        t.update( loc, secs, 0 );
+        double angle_deg = sun_angle( t, world_up );
         double diff = fabs( angle_deg - target_angle_deg );
         if ( diff < best_diff ) {
             if ( last_angle <= 180.0 && ascending