src/thgeomag.c

   1 /**
   2 * @file thgeomag.cxx
   3 */
   4
   5 /* Copyright (C) 2006 Martin Budaj
   6 *
   7 * based on GPL-licensed code by
   8 * Copyright (C) 2000  Edward A Williams <Ed_Williams@compuserve.com>
   9 *
  10 * --------------------------------------------------------------------
  11 * This program is free software; you can redistribute it and/or modify
  12 * it under the terms of the GNU General Public License as published by
  13 * the Free Software Foundation; either version 2 of the License, or
  14 * any later version.
  15 *
  16 * This program is distributed in the hope that it will be useful,
  17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19 * GNU General Public License for more details.
  20 *
  21 * You should have received a copy of the GNU General Public License
  22 * along with this program; if not, write to the Free Software
  23 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  24 * --------------------------------------------------------------------
  25 */
  26
  27 #include "thgeomag.h"
  28
  29 #include <math.h>
  30 #include "thgeomagdata.h"
  31
  32 #define max(a,b)        (((a) > (b)) ? (a) : (b))
  33
  34 /*struct magfield_ {
  35 double X, Y, Z;
  36 };
  37 magfield_ magfield;*/
  38
  39 #define nmax thgeomag_maxdeg
  40 #define nmaxl thgeomag_maxdeg
  41
  42 #define pi 3.14159265358979
  43 #define a 6378.137
  44 #define b 6356.7523142
  45 #define r_0 6371.2
  46
  47 double thgeomag(double lat, double lon, double h, double dat) {
  48
  49   int n,m;
  50
  51   static double P[nmax+1][nmax+1];
  52   static double DP[nmax+1][nmax+1];
  53   static double gnm[nmax+1][nmax+1];
  54   static double hnm[nmax+1][nmax+1];
  55   static double sm[nmax+1];
  56   static double cm[nmax+1];
  57
  58   static double root[nmax+1];
  59   static double roots[nmax+1][nmax+1][2];
  60
  61
  62   double yearfrac,sr,r,theta,c,s,psi,fn,fn_0,B_r,B_theta,B_phi,X,Y; /* Z */
  63   double sinpsi, cospsi, inv_s;
  64
  65   static int been_here = 0;
  66
  67   double sinlat = sin(lat);
  68   double coslat = cos(lat);
  69
  70   h = h / 1000;
  71
  72   /* convert to geocentric */
  73   sr = sqrt(a*a*coslat*coslat + b*b*sinlat*sinlat);
  74   /* sr is effective radius */
  75   theta = atan2(coslat * (h*sr + a*a), sinlat * (h*sr + b*b));
  76
  77   /* theta is geocentric co-latitude */
  78
  79   r = h*h + 2.0*h * sr +
  80     (a*a*a*a - ( a*a*a*a - b*b*b*b ) * sinlat*sinlat ) /
  81     (a*a - (a*a - b*b) * sinlat*sinlat );
  82
  83   r = sqrt(r);
  84
  85   /* r is geocentric radial distance */
  86   c = cos(theta);
  87   s = sin(theta);
  88   /* protect against zero divide at geographic poles */
  89   inv_s =  1.0 / (s + (s == 0.)*1.0e-8);
  90
  91   /*zero out arrays */
  92   for ( n = 0; n <= nmax; n++ ) {
  93     for ( m = 0; m <= n; m++ ) {
  94       P[n][m] = 0;
  95       DP[n][m] = 0;
  96     }
  97   }
  98
  99   /* diagonal elements */
 100   P[0][0] = 1;
 101   P[1][1] = s;
 102   DP[0][0] = 0;
 103   DP[1][1] = c;
 104   P[1][0] = c ;
 105   DP[1][0] = -s;
 106
 107   /* these values will not change for subsequent function calls */
 108   if( !been_here ) {
 109     for ( n = 2; n <= nmax; n++ ) {
 110       root[n] = sqrt((2.0*n-1) / (2.0*n));
 111     }
 112
 113     for ( m = 0; m <= nmax; m++ ) {
 114       double mm = m*m;
 115       for ( n = max(m + 1, 2); n <= nmax; n++ ) {
 116         roots[m][n][0] = sqrt((n-1)*(n-1) - mm);
 117         roots[m][n][1] = 1.0 / sqrt( n*n - mm);
 118       }
 119     }
 120     been_here = 1;
 121   }
 122
 123   for ( n=2; n <= nmax; n++ ) {
 124     /*  double root = sqrt((2.0*n-1) / (2.0*n)); */
 125     P[n][n] = P[n-1][n-1] * s * root[n];
 126     DP[n][n] = (DP[n-1][n-1] * s + P[n-1][n-1] * c) * root[n];
 127   }
 128
 129   /* lower triangle */
 130   for ( m = 0; m <= nmax; m++ ) {
 131     /*  double mm = m*m;  */
 132     for ( n = max(m + 1, 2); n <= nmax; n++ ) {
 133       /* double root1 = sqrt((n-1)*(n-1) - mm); */
 134       /* double root2 = 1.0 / sqrt( n*n - mm);  */
 135       P[n][m] = (P[n-1][m] * c * (2.0*n-1) -
 136         P[n-2][m] * roots[m][n][0]) * roots[m][n][1];
 137       DP[n][m] = ((DP[n-1][m] * c - P[n-1][m] * s) *
 138         (2.0*n-1) - DP[n-2][m] * roots[m][n][0]) * roots[m][n][1];
 139     }
 140   }
 141
 142   /* compute gnm, hnm at dat */
 143
 144   int mindex = (int)((dat - thgeomag_minyear) / thgeomag_step);
 145   if (mindex < 0) mindex = 0;
 146   if (mindex > thgeomag_maxmindex) mindex = thgeomag_maxmindex;
 147   yearfrac = dat - thgeomag_step*mindex - thgeomag_minyear;
 148
 149   for (n=1;n<=nmaxl;n++) {
 150     for (m = 0;m<=nmaxl;m++) {
 151       if (mindex == thgeomag_maxmindex) {
 152         gnm[n][m] = thgeomag_GNM[mindex][n][m] + yearfrac * thgeomag_GNMD[n][m];
 153         hnm[n][m] = thgeomag_HNM[mindex][n][m] + yearfrac * thgeomag_HNMD[n][m];
 154       } else {
 155         gnm[n][m] = thgeomag_GNM[mindex][n][m] + yearfrac / thgeomag_step * (thgeomag_GNM[mindex+1][n][m] - thgeomag_GNM[mindex][n][m]);
 156         hnm[n][m] = thgeomag_HNM[mindex][n][m] + yearfrac / thgeomag_step * (thgeomag_HNM[mindex+1][n][m] - thgeomag_HNM[mindex][n][m]);
 157       }
 158     }
 159   }
 160
 161   /* compute sm (sin(m lon) and cm (cos(m lon)) */
 162   for (m = 0;m<=nmaxl;m++) {
 163     sm[m] = sin(m * lon);
 164     cm[m] = cos(m * lon);
 165   }
 166
 167   /* compute B fields */
 168   B_r = 0.0;
 169   B_theta = 0.0;
 170   B_phi = 0.0;
 171   fn_0 = r_0/r;
 172   fn = fn_0 * fn_0;
 173
 174   for ( n = 1; n <= nmaxl; n++ ) {
 175     double c1_n=0;
 176     double c2_n=0;
 177     double c3_n=0;
 178     for ( m = 0; m <= n; m++ ) {
 179       double tmp = (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]);
 180       c1_n += tmp * P[n][m];
 181       c2_n += tmp * DP[n][m];
 182       c3_n +=  m * (gnm[n][m] * sm[m] - hnm[n][m] * cm[m]) * P[n][m];
 183     }
 184     /* fn=pow(r_0/r,n+2.0);   */
 185     fn *= fn_0;
 186     B_r += (n + 1) * c1_n * fn;
 187     B_theta -= c2_n * fn;
 188     B_phi += c3_n * fn * inv_s;
 189   }
 190
 191
 192
 193   /* Find geodetic field components: */
 194   psi = theta - (pi / 2.0 - lat);
 195   sinpsi = sin(psi);
 196   cospsi = cos(psi);
 197   X = -B_theta * cospsi - B_r * sinpsi;
 198   Y = B_phi;
 199   /* Z = B_theta * sinpsi - B_r * cospsi; */
 200
 201   /*    field[0]=B_r;
 202   field[1]=B_theta;
 203   field[2]=B_phi;
 204   field[3]=X;
 205   field[4]=Y;
 206   field[5]=Z;*/   /* output fields */
 207   /* find variation in radians */
 208   /* return zero variation at magnetic pole X=Y=0. */
 209   /* E is positive */
 210
 211   /*    magfield.X = X;
 212   magfield.Y = Y;
 213   magfield.Z = Z; */
 214
 215   return (X != 0. || Y != 0.) ? atan2(Y, X) : (double) 0.;
 216 }
 217