I found 3 problems with the GS modeling in flightgear (all my fault originally

I believe.) :-)

- The height of the navaid was not being properly converted to meters
  before being used in our internal calculations.  This caused the GS
  to be placed too high.

- I was using the wrong trig function to calculate the current approach
  angle of the aircraft.  The distance to the GS source is the euclidean
  point to point distance and represents the hypotenuse (not the ground
  distance) so I need to use asin() rather than atan() to calculate the
  angle.

- I was calculating distance directly to the GS source, rather than
  taking into consideration that the GS transmitter projects a plane,
  so I need to take the distance to the line where that plane intersectso
  the ground.  Previously, the way I modeled my distance calculation, the
  GS transmitter effectively formed a 3 degree cone from the source.  The GS
  transmitter is usually placed a 100 meters or so off the runway edge so
  the cone model could never bring you in to the touch down point precisely.

With these changes, the GS will bring you in precisely to the touchdown
point as defined in the default.ils.gz file (it wouldn't before.)  The only
issue that remains is that it will bring you in to the elevation defined
in the ILS database, which doesn't necessarily match the DEM/SRTM terrain
at that point.  Still on average, this will be a big improvement until we
can do a better job of getting the runway end elevations nailed correctly.

diff --git a/src/Cockpit/navcom.cxx b/src/Cockpit/navcom.cxx
index 3f8f21aee6db7666e3ec13ef1e18a98162a51e9b..1fe7b3882c23cf9c0776acd2a67f9db60ba0bd4f 100644 (file)
--- a/src/Cockpit/navcom.cxx
+++ b/src/Cockpit/navcom.cxx
@@ -29,6 +29,7 @@
 
 #include <simgear/compiler.h>
 #include <simgear/math/sg_random.h>
+#include <simgear/math/vector.hxx>
 
 #include <Aircraft/aircraft.hxx>
 #include <Navaids/ilslist.hxx>
@@ -325,9 +326,20 @@ FGNavCom::update(double dt)
        nav_loc_dist = aircraft.distance3D( station );
 
        if ( nav_has_gs ) {
-           station = Point3D( nav_gs_x, nav_gs_y, nav_gs_z );
-           nav_gs_dist = aircraft.distance3D( station );
-            // wgs84 heading to glide slope
+            // find closest distance to the gs base line
+            sgdVec3 p;
+            sgdSetVec3( p, aircraft.x(), aircraft.y(), aircraft.z() );
+            sgdVec3 p0;
+            sgdSetVec3( p0, nav_gs_x, nav_gs_y, nav_gs_z );
+            double dist = sgdClosestPointToLineDistSquared( p, p0,
+                                                            gs_base_vec );
+            nav_gs_dist = sqrt( dist );
+            // cout << nav_gs_dist;
+
+            Point3D tmp( nav_gs_x, nav_gs_y, nav_gs_z );
+            // cout << " (" << aircraft.distance3D( tmp ) << ")" << endl;
+
+            // wgs84 heading to glide slope (to determine sign of distance)
             geo_inverse_wgs_84( elev,
                                 lat * SGD_RADIANS_TO_DEGREES,
                                 lon * SGD_RADIANS_TO_DEGREES, 
@@ -491,6 +503,22 @@ void FGNavCom::search()
            nav_gs_y = ils->get_gs_y();
            nav_gs_z = ils->get_gs_z();
 
+            // derive GS baseline
+            double tlon, tlat, taz;
+            geo_direct_wgs_84 ( 0.0, nav_gslat, nav_gslon, nav_radial + 90,  
+                                100.0, &tlat, &tlon, &taz );
+            cout << nav_gslon << "," << nav_gslat << "  "
+                 << tlon << "," << tlat << "  (" << nav_elev << ")" << endl;
+            Point3D p1 = sgGeodToCart( Point3D(tlon*SGD_DEGREES_TO_RADIANS,
+                                               tlat*SGD_DEGREES_TO_RADIANS,
+                                               nav_elev*SG_FEET_TO_METER) );
+            cout << nav_gs_x << "," << nav_gs_y << "," << nav_gs_z << endl;
+            cout << p1 << endl;
+            sgdSetVec3( gs_base_vec,
+                        p1.x()-nav_gs_x, p1.y()-nav_gs_y, p1.z()-nav_gs_z );
+            cout << gs_base_vec[0] << "," << gs_base_vec[1] << ","
+                 << gs_base_vec[2] << endl;
+
            if ( globals->get_soundmgr()->exists( nav_fx_name ) ) {
                globals->get_soundmgr()->remove( nav_fx_name );
            }
@@ -625,7 +653,8 @@ double FGNavCom::get_nav_gs_needle_deflection() const {
        double x = nav_gs_dist;
        double y = (fgGetDouble("/position/altitude-ft") - nav_elev)
             * SG_FEET_TO_METER;
-       double angle = atan2( y, x ) * SGD_RADIANS_TO_DEGREES;
+        // cout << "dist = " << x << " height = " << y << endl;
+       double angle = asin( y / x ) * SGD_RADIANS_TO_DEGREES;
        return (nav_target_gs - angle) * 5.0;
     } else {
        return 0.0;

diff --git a/src/Cockpit/navcom.hxx b/src/Cockpit/navcom.hxx
index 33a01dcf510606df305a3c3d27ad69233c98fdb1..960d8cdbc9c09da1432a96cbd6dbc06fdddefbff 100644 (file)
--- a/src/Cockpit/navcom.hxx
+++ b/src/Cockpit/navcom.hxx
@@ -102,6 +102,7 @@ class FGNavCom : public FGSubsystem
     double nav_gs_x;
     double nav_gs_y;
     double nav_gs_z;
+    sgdVec3 gs_base_vec;
     double nav_gs_dist;
     double nav_gs_dist_signed;
     SGTimeStamp prev_time;

diff --git a/src/Navaids/ils.hxx b/src/Navaids/ils.hxx
index 953132d9feac244edd7d97c20cd5904418a2bcfe..a9b7da167d1294f3b718468bf2038e070f154046 100644 (file)
--- a/src/Navaids/ils.hxx
+++ b/src/Navaids/ils.hxx
@@ -179,7 +179,9 @@ operator >> ( istream& in, FGILS& i )
     // generate cartesian coordinates
     Point3D geod, cart;
 
-    geod = Point3D( i.loclon * SGD_DEGREES_TO_RADIANS, i.loclat * SGD_DEGREES_TO_RADIANS, i.gselev );
+    geod = Point3D( i.loclon * SGD_DEGREES_TO_RADIANS,
+                    i.loclat * SGD_DEGREES_TO_RADIANS,
+                    i.gselev * SG_FEET_TO_METER );
     cart = sgGeodToCart( geod );
     i.x = cart.x();
     i.y = cart.y();
@@ -190,7 +192,9 @@ operator >> ( istream& in, FGILS& i )
     } else {
        i.has_gs = true;
 
-       geod = Point3D( i.gslon * SGD_DEGREES_TO_RADIANS, i.gslat * SGD_DEGREES_TO_RADIANS, i.gselev );
+       geod = Point3D( i.gslon * SGD_DEGREES_TO_RADIANS,
+                        i.gslat * SGD_DEGREES_TO_RADIANS,
+                        i.gselev * SG_FEET_TO_METER );
        cart = sgGeodToCart( geod );
        i.gs_x = cart.x();
        i.gs_y = cart.y();
@@ -203,7 +207,9 @@ operator >> ( istream& in, FGILS& i )
     } else {
        i.has_dme = true;
 
-       geod = Point3D( i.dmelon * SGD_DEGREES_TO_RADIANS, i.dmelat * SGD_DEGREES_TO_RADIANS, i.gselev);
+       geod = Point3D( i.dmelon * SGD_DEGREES_TO_RADIANS,
+                        i.dmelat * SGD_DEGREES_TO_RADIANS,
+                        i.gselev * SG_FEET_TO_METER );
        cart = sgGeodToCart( geod );
        i.dme_x = cart.x();
        i.dme_y = cart.y();