simgear/magvar/magvar.cxx

   1 // magvar.cxx -- compute local magnetic variation given position,
   2 //               altitude, and date
   3 //
   4 // This is an implimentation of the NIMA WMM 2000
   5 //
   6 //    http://www.nima.mil/GandG/ngdc-wmm2000.html
   7 //
   8 // Copyright (C) 2000  Edward A Williams <Ed_Williams@compuserve.com>
   9 //
  10 // Adapted from Excel 3.0 version 3/27/94 EAW
  11 // Recoded in C++ by Starry Chan
  12 // WMM95 added and rearranged in ANSI-C EAW 7/9/95
  13 // Put shell around program and made Borland & GCC compatible EAW 11/22/95
  14 // IGRF95 added 2/96 EAW
  15 // WMM2000 IGR2000 added 2/00 EAW
  16 // Released under GPL 3/26/00 EAW
  17 // Adaptions and modifications for the SimGear project  3/27/2000 CLO
  18 //
  19 // This program is free software; you can redistribute it and/or
  20 // modify it under the terms of the GNU General Public License as
  21 // published by the Free Software Foundation; either version 2 of the
  22 // License, or (at your option) any later version.
  23 //
  24 // This program is distributed in the hope that it will be useful, but
  25 // WITHOUT ANY WARRANTY; without even the implied warranty of
  26 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  27 // General Public License for more details.
  28 //
  29 // You should have received a copy of the GNU General Public License
  30 // along with this program; if not, write to the Free Software
  31 // Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  32 //
  33 // $Id$
  34
  35
  36 #include <stdio.h>
  37 #include <stdlib.h>
  38 #include <math.h>
  39
  40 #include "magvar.hxx"
  41
  42
  43 #define max(a,b)        (((a) > (b)) ? (a) : (b))
  44
  45 static const double pi = 3.14159265358979;
  46 static const double a = 6378.16;        /* major radius (km) IAU66 ellipsoid */
  47 static const double f = 1.0 / 298.25;   /* inverse flattening IAU66 ellipsoid */
  48 static const double b = 6378.16 * (1.0 -1.0 / 298.25 );
  49 /* minor radius b=a*(1-f) */
  50 static const double r_0 = 6371.2;       /* "mean radius" for spherical harmonic expansion */
  51
  52 static double gnm_wmm2000[13][13] =
  53 {
  54     {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},
  55     {-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},
  56     {-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},
  57     {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},
  58     {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},
  59     {-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},
  60     {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},
  61     {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},
  62     {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},
  63     {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},
  64     {-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},
  65     {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},
  66     {-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},
  67 };
  68
  69 static double hnm_wmm2000[13][13]=
  70 {
  71     {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},
  72     {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},
  73     {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},
  74     {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},
  75     {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},
  76     {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},
  77     {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},
  78     {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},
  79     {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},
  80     {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},
  81     {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},
  82     {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},
  83     {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},
  84 };
  85
  86 static double gtnm_wmm2000[13][13]=
  87 {
  88     {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},
  89     {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},
  90     {-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},
  91     {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},
  92     {-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},
  93     {-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},
  94     {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},
  95     {-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},
  96     {-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},
  97     {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},
  98     {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},
  99     {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},
 100     {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},
 101 };
 102
 103 static double htnm_wmm2000[13][13]=
 104 {
 105     {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},
 106     {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},
 107     {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},
 108     {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},
 109     {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},
 110     {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},
 111     {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},
 112     {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},
 113     {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},
 114     {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},
 115     {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},
 116     {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},
 117     {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},
 118 };
 119
 120 static const int nmax = 12;
 121
 122 static double P[13][13];
 123 static double DP[13][13];
 124 static double gnm[13][13];
 125 static double hnm[13][13];
 126 static double sm[13];
 127 static double cm[13];
 128
 129
 130 /* Convert date to Julian day    1950-2049 */
 131 unsigned long int yymmdd_to_julian_days( int yy, int mm, int dd )
 132 {
 133     unsigned long jd;
 134
 135     yy = (yy < 50) ? (2000 + yy) : (1900 + yy);
 136     jd = dd - 32075L + 1461L * (yy + 4800L + (mm - 14) / 12 ) / 4;
 137     jd = jd + 367L * (mm - 2 - (mm - 14) / 12*12) / 12;
 138     jd = jd - 3 * ((yy + 4900L + (mm - 14) / 12) / 100) / 4;
 139
 140     /* printf("julian date = %d\n", jd ); */
 141     return jd;
 142 }
 143
 144
 145 /* Convert degrees to radians */
 146 double deg_to_rad( double deg )
 147 {
 148     return deg*pi/180.;
 149 }
 150
 151
 152 /* Convert radians to degrees */
 153 double rad_to_deg( double rad )
 154 {
 155     return rad*180./pi;
 156 }
 157
 158
 159 /* return variation (in degrees) given geodetic latitude (radians), longitude
 160 (radians) ,height (km) and (Julian) date
 161 N and E lat and long are positive, S and W negative
 162 */
 163
 164 double SGMagVar( double lat, double lon, double h, long dat, double* field )
 165 {
 166     /* output field B_r,B_th,B_phi,B_x,B_y,B_z */
 167     int n,m;
 168     /* reference dates */
 169     long date0_wmm2000 = yymmdd_to_julian_days(0,1,1);
 170
 171     double yearfrac,sr,r,theta,c,s,psi,fn,B_r,B_theta,B_phi,X,Y,Z;
 172
 173     /* convert to geocentric coords: */
 174     sr = sqrt(pow(a*cos(lat),2.0)+pow(b*sin(lat),2.0));
 175     /* sr is effective radius */
 176     theta = atan2(cos(lat) * (h * sr + a * a),
 177                   sin(lat) * (h * sr + b * b));
 178     /* theta is geocentric co-latitude */
 179
 180     r = h * h + 2.0*h * sr +
 181         (pow(a,4.0) - (pow(a,4.0) - pow(b,4.0)) * pow(sin(lat),2.0)) /
 182         (a * a - (a * a - b * b) * pow(sin(lat),2.0));
 183
 184     r = sqrt(r);
 185
 186     /* r is geocentric radial distance */
 187     c = cos(theta);
 188     s = sin(theta);
 189
 190     /*zero out arrays */
 191     for ( n = 0; n <= nmax; n++ ) {
 192         for ( m = 0; m <= n; m++ ) {
 193             P[n][m] = 0;
 194             DP[n][m] = 0;
 195         }
 196     }
 197
 198     /* diagonal elements */
 199     P[0][0] = 1;
 200     P[1][1] = s;
 201     DP[0][0] = 0;
 202     DP[1][1] = c;
 203     P[1][0] = c ;
 204     DP[1][0] = -s;
 205
 206     for ( n = 2; n <= nmax; n++ ) {
 207         P[n][n] = P[n-1][n-1] * s * sqrt((2.0*n-1) / (2.0*n));
 208         DP[n][n] = (DP[n-1][n-1] * s + P[n-1][n-1] * c) *
 209             sqrt((2.0*n-1) / (2.0*n));
 210     }
 211
 212     /* lower triangle */
 213     for ( m = 0; m <= nmax; m++ ) {
 214         for ( n = max(m + 1, 2); n <= nmax; n++ ) {
 215             P[n][m] = (P[n-1][m] * c * (2.0*n-1) - P[n-2][m] *
 216                        sqrt(1.0*(n-1)*(n-1) - m * m)) /
 217                 sqrt(1.0* n * n - m * m);
 218             DP[n][m] = ((DP[n-1][m] * c - P[n-1][m] * s) *
 219                         (2.0*n-1) - DP[n-2][m] *
 220                         sqrt(1.0*(n-1) * (n-1) - m * m)) /
 221                 sqrt(1.0* n * n - m * m);
 222         }
 223     }
 224
 225     /* compute gnm, hnm at dat */
 226     /* WMM2000 */
 227     yearfrac = (dat - date0_wmm2000) / 365.25;
 228     for ( n = 1; n <= nmax; n++ ) {
 229         for ( m = 0; m <= nmax; m++ ) {
 230             gnm[n][m] = gnm_wmm2000[n][m] + yearfrac * gtnm_wmm2000[n][m];
 231             hnm[n][m] = hnm_wmm2000[n][m] + yearfrac * htnm_wmm2000[n][m];
 232         }
 233     }
 234
 235     /* compute sm (sin(m lon) and cm (cos(m lon)) */
 236     for ( m = 0; m <= nmax; m++ ) {
 237         sm[m] = sin(m * lon);
 238         cm[m] = cos(m * lon);
 239     }
 240
 241     /* compute B fields */
 242     B_r = 0.0;
 243     B_theta = 0.0;
 244     B_phi = 0.0;
 245
 246     for ( n = 1; n <= nmax; n++ ) {
 247         double c1_n=0;
 248         double c2_n=0;
 249         double c3_n=0;
 250         for ( m = 0; m <= n; m++ ) {
 251             c1_n=c1_n + (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]) * P[n][m];
 252             c2_n=c2_n + (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]) * DP[n][m];
 253             c3_n=c3_n + m * (gnm[n][m] * sm[m] - hnm[n][m] * cm[m]) * P[n][m];
 254         }
 255         fn=pow(r_0/r,n+2.0);
 256         B_r = B_r + (n + 1) * c1_n * fn;
 257         B_theta = B_theta - c2_n * fn;
 258         B_phi = B_phi + c3_n * fn / s;
 259     }
 260
 261     /* Find geodetic field components: */
 262     psi = theta - (pi / 2.0 - lat);
 263     X = -B_theta * cos(psi) - B_r * sin(psi);
 264     Y = B_phi;
 265     Z = B_theta * sin(psi) - B_r * cos(psi);
 266
 267     field[0]=B_r;
 268     field[1]=B_theta;
 269     field[2]=B_phi;
 270     field[3]=X;
 271     field[4]=Y;
 272     field[5]=Z;   /* output fields */
 273
 274     /* find variation, convert to degrees! */
 275     return rad_to_deg(atan2(Y, X));  /* E is positive */
 276 }
 277