]> git.mxchange.org Git - simgear.git/blobdiff - simgear/magvar/coremag.cxx
projects / simgear.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
Performance optimization: empty() instead of size()>0
[simgear.git] / simgear / magvar / coremag.cxx
diff --git a/simgear/magvar/coremag.cxx b/simgear/magvar/coremag.cxx
index b29634a2882b08c11a658e150a8f0d5bf0fa7072..e6d42c1d1a5bf45341a2b476c19ea31ee9913337 100644 (file)
--- a/simgear/magvar/coremag.cxx
+++ b/simgear/magvar/coremag.cxx
@@ -23,6 +23,19 @@
 //
 // Put in some bullet-proofing to handle magnetic and geographic poles.
 // 3/28/2000 EAW
+//
+// Updated coefficient arrays to use the current WMM2005 model,
+// (valid between 2005.0 and 2010.0)
+// Also removed unused variables and corrected earth radii constants
+// to the values for WGS84 and WMM2005.
+// Reference:
+//     McLean, S., S. Macmillan, S. Maus, V. Lesur, A.
+//     Thomson, and D. Dater, December 2004, The
+//     US/UK World Magnetic Model for 2005-2010,
+//     NOAA Technical Report NESDIS/NGDC-1.
+//
+// 25/10/2006  Wim Van Hoydonck -- wim.van.hoydonck@gmail.com
+//
 
 //  The routine uses a spherical harmonic expansion of the magnetic
 // potential up to twelfth order, together with its time variation, as
@@ -57,10 +70,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// 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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 
@@ -70,81 +82,76 @@
 #include <math.h>
 
 #include <simgear/constants.h>
+#include <simgear/sg_inlines.h>
 
 #include "coremag.hxx"
 
+static const double a = 6378.137;       /* semi-major axis (equatorial radius) of WGS84 ellipsoid */
+static const double b = 6356.7523142;   /* semi-minor axis referenced to the WGS84 ellipsoid */
+static const double r_0 = 6371.2;      /* standard Earth magnetic reference radius  */
 
-#define        max(a,b)        (((a) > (b)) ? (a) : (b))
-
-static const double pi = 3.14159265358979;
-static const double a = 6378.16;       /* major radius (km) IAU66 ellipsoid */
-static const double f = 1.0 / 298.25;  /* inverse flattening IAU66 ellipsoid */
-static const double b = 6378.16 * (1.0 -1.0 / 298.25 );
-/* minor radius b=a*(1-f) */
-static const double r_0 = 6371.2;      /* "mean radius" for spherical harmonic expansion */
-
-static double gnm_wmm2000[13][13] =
+static double gnm_wmm2005[13][13] =
 {
     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {-29616.0, -1722.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {-2266.7, 3070.2, 1677.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {1322.4, -2291.5, 1255.9, 724.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {932.1, 786.3, 250.6, -401.5, 106.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {-211.9, 351.6, 220.8, -134.5, -168.8, -13.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {73.8, 68.2, 74.1, -163.5, -3.8, 17.1, -85.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {77.4, -73.9, 2.2, 35.7, 7.3, 5.2, 8.4, -1.5, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {23.3, 7.3, -8.5, -6.6, -16.9, 8.6, 4.9, -7.8, -7.6, 0.0, 0.0, 0.0, 0.0},
-    {5.7, 8.5, 2.0, -9.8, 7.6, -7.0, -2.0, 9.2, -2.2, -6.6, 0.0, 0.0, 0.0},
-    {-2.2, -5.7, 1.6, -3.7, -0.6, 4.1, 2.2, 2.2, 4.6, 2.3, 0.1, 0.0, 0.0},
-    {3.3, -1.1, -2.4, 2.6, -1.3, -1.7, -0.6, 0.4, 0.7, -0.3, 2.3, 4.2, 0.0},
-    {-1.5, -0.2, -0.3, 0.5, 0.2, 0.9, -1.4, 0.6, -0.6, -1.0, -0.3, 0.3, 0.4},
+    {-29556.8, -1671.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {-2340.6, 3046.9, 1657.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {1335.4, -2305.1, 1246.7, 674.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {919.8, 798.1, 211.3, -379.4, 100.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {-227.4, 354.6, 208.7, -136.5, -168.3, -14.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {73.2, 69.7, 76.7, -151.2, -14.9, 14.6, -86.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {80.1, -74.5, -1.4, 38.5, 12.4, 9.5, 5.7, 1.8, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {24.9, 7.7, -11.6, -6.9, -18.2, 10.0, 9.2, -11.6, -5.2, 0.0, 0.0, 0.0, 0.0},
+    {5.6, 9.9, 3.5, -7.0, 5.1, -10.8, -1.3, 8.8, -6.7, -9.1, 0.0, 0.0, 0.0},
+    {-2.3, -6.3, 1.6, -2.6, 0.0, 3.1, 0.4, 2.1, 3.9, -0.1, -2.3, 0.0, 0.0},
+    {2.8, -1.6, -1.7, 1.7, -0.1, 0.1, -0.7, 0.7, 1.8, 0.0, 1.1, 4.1, 0.0},
+    {-2.4, -0.4, 0.2, 0.8, -0.3, 1.1, -0.5, 0.4, -0.3, -0.3, -0.1, -0.3, -0.1},
 };
 
-static double hnm_wmm2000[13][13]=
+static double hnm_wmm2005[13][13]=
 {
     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {0.0, 5194.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {0.0, -2484.8, -467.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {0.0, -224.7, 293.0, -486.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {0.0, 273.3, -227.9, 120.9, -302.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {0.0, 42.0, 173.8, -135.0, -38.6, 105.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {0.0, -17.4, 61.2, 63.2, -62.9, 0.2, 43.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {0.0, -62.3, -24.5, 8.9, 23.4, 15.0, -27.6, -7.8, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {0.0, 12.4, -20.8, 8.4, -21.2, 15.5, 9.1, -15.5, -5.4, 0.0, 0.0, 0.0, 0.0},
-    {0.0, -20.4, 13.9, 12.0, -6.2, -8.6, 9.4, 5.0, -8.4, 3.2, 0.0, 0.0, 0.0},
-    {0.0, 0.9, -0.7, 3.9, 4.8, -5.3, -1.0, -2.4, 1.3, -2.3, -6.4, 0.0, 0.0},
-    {0.0, -1.5, 0.7, -1.1, -2.3, 1.3, -0.6, -2.8, -1.6, -0.1, -1.9, 1.4, 0.0},
-    {0.0, -1.0, 0.7, 2.2, -2.5, -0.2, 0.0, -0.2, 0.0, 0.2, -0.9, -0.2, 1.0},
+    {0.0, 5079.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {0.0, -2594.7, -516.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {0.0, -199.9, 269.3, -524.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {0.0, 281.5, -226.0, 145.8, -304.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {0.0, 42.4, 179.8, -123.0, -19.5, 103.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {0.0, -20.3, 54.7, 63.6, -63.4, -0.1, 50.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {0.0, -61.5, -22.4, 7.2, 25.4, 11.0, -26.4, -5.1, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {0.0, 11.2, -21.0, 9.6, -19.8, 16.1, 7.7, -12.9, -0.2, 0.0, 0.0, 0.0, 0.0},
+    {0.0, -20.1, 12.9, 12.6, -6.7, -8.1, 8.0, 2.9, -7.9, 6.0, 0.0, 0.0, 0.0},
+    {0.0, 2.4, 0.2, 4.4, 4.8, -6.5, -1.1, -3.4, -0.8, -2.3, -7.9, 0.0, 0.0},
+    {0.0, 0.3, 1.2, -0.8, -2.5, 0.9, -0.6, -2.7, -0.9, -1.3, -2.0, -1.2, 0.0},
+    {0.0, -0.4, 0.3, 2.4, -2.6, 0.6, 0.3, 0.0, 0.0, 0.3, -0.9, -0.4, 0.8},
 };
 
-static double gtnm_wmm2000[13][13]=
+static double gtnm_wmm2005[13][13]=
 {
     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {14.7, 11.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {-13.6, -0.7, -1.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {0.3, -4.3, 0.9, -8.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {-1.6, 0.9, -7.6, 2.2, -3.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {-0.9, -0.2, -2.5, -2.7, -0.9, 1.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {1.2, 0.2, 1.7, 1.6, -0.1, -0.3, 0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {-0.4, -0.8, -0.2, 1.1, 0.4, 0.0, -0.2, -0.2, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {-0.3, 0.6, -0.8, 0.3, -0.2, 0.5, 0.0, -0.6, 0.1, 0.0, 0.0, 0.0, 0.0},
+    {8.0, 10.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {-15.1, -7.8, -0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {0.4, -2.6, -1.2, -6.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {-2.5, 2.8, -7.0, 6.2, -3.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {-2.8, 0.7, -3.2, -1.1, 0.1, -0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {-0.7, 0.4, -0.3, 2.3, -2.1, -0.6, 1.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {0.2, -0.1, -0.3, 1.1, 0.6, 0.5, -0.4, 0.6, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {0.1, 0.3, -0.4, 0.3, -0.3, 0.2, 0.4, -0.7, 0.4, 0.0, 0.0, 0.0, 0.0},
     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
 };
 
-static double htnm_wmm2000[13][13]=
+static double htnm_wmm2005[13][13]=
 {
     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {0.0, -20.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {0.0, -21.5, -9.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {0.0, 6.4, -1.3, -13.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {0.0, 2.3, 0.7, 3.7, -0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {0.0, 0.0, 2.1, 2.3, 3.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {0.0, -0.3, -1.7, -0.9, -1.0, -0.1, 1.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {0.0, 1.4, 0.2, 0.7, 0.4, -0.3, -0.8, -0.1, 0.0, 0.0, 0.0, 0.0, 0.0},
-    {0.0, -0.5, 0.1, -0.2, 0.0, 0.1, -0.1, 0.3, 0.2, 0.0, 0.0, 0.0, 0.0},
+    {0.0, -20.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {0.0, -23.2, -14.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {0.0, 5.0, -7.0, -0.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {0.0, 2.2, 1.6, 5.8, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {0.0, 0.0, 1.7, 2.1, 4.8, -1.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {0.0, -0.6, -1.9, -0.4, -0.5, -0.3, 0.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {0.0, 0.6, 0.4, 0.2, 0.3, -0.8, -0.2, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0},
+    {0.0, -0.2, 0.1, 0.3, 0.4, 0.1, -0.2, 0.4, 0.4, 0.0, 0.0, 0.0, 0.0},
     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
     {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
@@ -188,8 +195,8 @@ double calc_magvar( double lat, double lon, double h, long dat, double* field )
 {
     /* output field B_r,B_th,B_phi,B_x,B_y,B_z */
     int n,m;
-    /* reference dates */
-    long date0_wmm2000 = yymmdd_to_julian_days(0,1,1);
+    /* reference date for current model is 1 januari 2005 */
+    long date0_wmm2005 = yymmdd_to_julian_days(5,1,1);
 
     double yearfrac,sr,r,theta,c,s,psi,fn,fn_0,B_r,B_theta,B_phi,X,Y,Z;
     double sinpsi, cospsi, inv_s;
@@ -243,7 +250,7 @@ double calc_magvar( double lat, double lon, double h, long dat, double* field )
 
        for ( m = 0; m <= nmax; m++ ) {
            double mm = m*m;
-           for ( n = max(m + 1, 2); n <= nmax; n++ ) {
+           for ( n = SG_MAX2(m + 1, 2); n <= nmax; n++ ) {
                roots[m][n][0] = sqrt((n-1)*(n-1) - mm);
                roots[m][n][1] = 1.0 / sqrt( n*n - mm);
            }
@@ -261,7 +268,7 @@ double calc_magvar( double lat, double lon, double h, long dat, double* field )
     /* lower triangle */
     for ( m = 0; m <= nmax; m++ ) {
        // double mm = m*m;
-       for ( n = max(m + 1, 2); n <= nmax; n++ ) {
+       for ( n = SG_MAX2(m + 1, 2); n <= nmax; n++ ) {
            // double root1 = sqrt((n-1)*(n-1) - mm);
            // double root2 = 1.0 / sqrt( n*n - mm);
            P[n][m] = (P[n-1][m] * c * (2.0*n-1) -
@@ -274,13 +281,14 @@ double calc_magvar( double lat, double lon, double h, long dat, double* field )
        }
     }
 
-    /* compute gnm, hnm at dat */
-    /* WMM2000 */
-    yearfrac = (dat - date0_wmm2000) / 365.25;
+    /* compute Gauss coefficients gnm and hnm of degree n and order m for the desired time
+       achieved by adjusting the coefficients at time t0 for linear secular variation */
+    /* WMM2005 */
+    yearfrac = (dat - date0_wmm2005) / 365.25;
     for ( n = 1; n <= nmax; n++ ) {
        for ( m = 0; m <= nmax; m++ ) {
-           gnm[n][m] = gnm_wmm2000[n][m] + yearfrac * gtnm_wmm2000[n][m];
-           hnm[n][m] = hnm_wmm2000[n][m] + yearfrac * htnm_wmm2000[n][m];
+           gnm[n][m] = gnm_wmm2005[n][m] + yearfrac * gtnm_wmm2005[n][m];
+           hnm[n][m] = hnm_wmm2005[n][m] + yearfrac * htnm_wmm2005[n][m];
        }
     }
 
@@ -315,7 +323,7 @@ double calc_magvar( double lat, double lon, double h, long dat, double* field )
     }
 
     /* Find geodetic field components: */
-    psi = theta - ((pi / 2.0) - lat);
+    psi = theta - ((M_PI / 2.0) - lat);
     sinpsi = sin(psi);
     cospsi = cos(psi);
     X = -B_theta * cospsi - B_r * sinpsi;
@@ -342,7 +350,7 @@ double SGMagVarOrig( double lat, double lon, double h, long dat, double* field )
     /* output field B_r,B_th,B_phi,B_x,B_y,B_z */
     int n,m;
     /* reference dates */
-    long date0_wmm2000 = yymmdd_to_julian_days(0,1,1);
+    long date0_wmm2005 = yymmdd_to_julian_days(5,1,1);
 
     double yearfrac,sr,r,theta,c,s,psi,fn,B_r,B_theta,B_phi,X,Y,Z;
 
@@ -387,7 +395,7 @@ double SGMagVarOrig( double lat, double lon, double h, long dat, double* field )
 
     /* lower triangle */
     for ( m = 0; m <= nmax; m++ ) {
-       for ( n = max(m + 1, 2); n <= nmax; n++ ) {
+       for ( n = SG_MAX2(m + 1, 2); n <= nmax; n++ ) {
            P[n][m] = (P[n-1][m] * c * (2.0*n-1) - P[n-2][m] *
                       sqrt(1.0*(n-1)*(n-1) - m * m)) /
                sqrt(1.0* n * n - m * m);
@@ -399,12 +407,12 @@ double SGMagVarOrig( double lat, double lon, double h, long dat, double* field )
     }
 
     /* compute gnm, hnm at dat */
-    /* WMM2000 */
-    yearfrac = (dat - date0_wmm2000) / 365.25;
+    /* WMM2005 */
+    yearfrac = (dat - date0_wmm2005) / 365.25;
     for ( n = 1; n <= nmax; n++ ) {
        for ( m = 0; m <= nmax; m++ ) {
-           gnm[n][m] = gnm_wmm2000[n][m] + yearfrac * gtnm_wmm2000[n][m];
-           hnm[n][m] = hnm_wmm2000[n][m] + yearfrac * htnm_wmm2000[n][m];
+           gnm[n][m] = gnm_wmm2005[n][m] + yearfrac * gtnm_wmm2005[n][m];
+           hnm[n][m] = hnm_wmm2005[n][m] + yearfrac * htnm_wmm2005[n][m];
        }
     }
 
Personal clone of SimGear, the library for FlightGear.