source: git/src/thgeomag.c @ 76cf7f1

RELEASE/1.2debug-cidebug-ci-sanitiserswalls-data
Last change on this file since 76cf7f1 was ac9afd6, checked in by Olly Betts <olly@…>, 9 years ago

src/thgeomag.c,src/thgeomag.h: Make thgeomag() return declination in
radians, since that's what we want and it's more consistent with it
taking lat and lon in radians.

  • Property mode set to 100644
File size: 5.7 KB
Line 
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_ {
35double X, Y, Z;
36};
37magfield_ 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
47double 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
Note: See TracBrowser for help on using the repository browser.