ARTS  2.3.1285(git:92a29ea9-dirty)
geomag_calc.cc
Go to the documentation of this file.
1 /* Copyright (C) 2003-2012 Nikolay Koulev <nkoulev@uni-bremen.de>
2 
3  This program is free software; you can redistribute it and/or
4  modify it under the terms of the GNU General Public License as
5  published by the Free Software Foundation; either version 2 of the
6  License, or (at your option) any later version.
7 
8  This program is distributed in the hope that it will be useful,
9  but WITHOUT ANY WARRANTY; without even the implied warranty of
10  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11  GNU General Public License for more details.
12 
13  You should have received a copy of the GNU General Public License
14  along with this program; if not, write to the Free Software
15  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
16  USA. */
17 
27 #include <iostream>
28 #include <cmath>
29 #include "arts.h"
30 #include "matpackI.h"
31 #include "xml_io.h"
32 #include "legendre.h"
33 
34 
35 extern const Numeric PI;
36 extern const Numeric DEG2RAD;
37 extern const Numeric EARTH_RADIUS;
38 
39 void magfield_nk(// Output
40  Numeric& B_r, // radial component of the geomagnetic field in [nT].
41  Numeric& B_th, // latitudinal component of the geomagnetic field in [nT].
42  Numeric& B_ph, // longitudinal component of the geomagnetic field in [nT].
43 
44  // Input
45  const Numeric r, // radial distance to the point in [km]
46  const Numeric theta, // geocentric colatitude of the point in [deg]
47  const Numeric phi, // longitude of the point in [deg].
48  // All coordinates - geocentric!
49 
50  const Index Ny, // number of elapsed years after an epoch year, J - [0,4]
51  const Verbosity& verbosity
52  )
53 
54 {
55 
56  Numeric Phi = phi * DEG2RAD; // Longitude angle in radian.
57  Numeric Theta = PI/2 - theta * DEG2RAD; // Colatitude angle in radian.
58 
59  // Initializing values of the magnetic field components.
60  B_r = 0;
61  B_th = 0;
62  B_ph = 0;
63 
64  Matrix M;
65 
66  xml_read_from_file ("geomag_coefficients.xml", M, verbosity);
67 
68  // M(i,0) and M(i,1) - the vectors with the values of the first and second coefficients
69  // of the IGRF model.
70  // M(i,2) and M(i,3) - the vectors with the values of the anual rate of change of the
71  // first and second coefficient of of the IGRF model.
72 
73 
74  // Loop over the degree number l of the Legendre polynomes.
75  for (Index l = 1; l <= 10; l++)
76  {
77  // Loop over the order number m of the Legendre polynomes.
78  for (Index m = 0; m <= l; m++)
79  {
80 
81  // Relating the row index in M to the coresponding
82  // degree number l and order number m.
83  Index j = l * (l + 1) / 2 + m - 1;
84 
85  // Calculating the associated Schmidt quasi-normalized Legendre
86  // polynomial for a degree number l and order number m.
87  Numeric P_lm =
88  g_legendre_poly_norm_schmidt (l, m, cos(Theta));
89 
90  // Calculating the derivative of the associated Schmidt quasi-normalized
91  // Legendre polynomial for a degree number l and order number m.
92  Numeric dP_lm =
93  g_legendre_poly_norm_schmidt_deriv3 (l, m, cos(Theta));
94 
95 
96  // Calculating the radial (upward) component of the magnetic field.
97  B_r += pow((Numeric)(l + 2), EARTH_RADIUS / r) * (Numeric)(l + 1) *
98  ((M(j,0) + (Numeric)Ny * M(j,2)) * cos((Numeric)m * Phi)
99  + (M(j,1) + (Numeric)Ny * M(j,3)) * sin((Numeric)m * Phi))
100  * P_lm;
101 
102  // Calculating the latitudinal (southward) component of the magnetic field.
103  B_th += pow(l + 2, EARTH_RADIUS / r) *
104  ((M(j,0) + (Numeric)Ny * M(j,2)) * cos((Numeric)m * Phi)
105  + (M(j,1) + (Numeric)Ny * M(j,3)) * sin((Numeric)m * Phi)) *
106  dP_lm * sin(Theta);
107 
108 
109 
110  // Calculating the longitudinal (eastward) component of the magnetic field.
111  B_ph += pow(l + 2, EARTH_RADIUS / r) * (Numeric)m *
112  ((M(j,0) + (Numeric)Ny * M(j,2)) * sin((Numeric)m * Phi)
113  - (M(j,1) + (Numeric)Ny * M(j,3)) * cos((Numeric)m * Phi)) *
114  P_lm / sin(Theta);
115 
116 
117  }
118  }
119 
120 }
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
#define M
Definition: rng.cc:165
const Numeric EARTH_RADIUS
Contains the code to calculate Legendre polynomials.
This file contains basic functions to handle XML data files.
Numeric g_legendre_poly_norm_schmidt(Index l, Index m, Numeric x)
g_legendre_poly_norm_schmidt
Definition: legendre.cc:355
The global header file for ARTS.
Numeric g_legendre_poly_norm_schmidt_deriv3(Index l, Index m, Numeric x)
g_legendre_poly_norm_schmidt_deriv3
Definition: legendre.cc:663
void magfield_nk(Numeric &B_r, Numeric &B_th, Numeric &B_ph, const Numeric r, const Numeric theta, const Numeric phi, const Index Ny, const Verbosity &verbosity)
Definition: geomag_calc.cc:39
void xml_read_from_file(const String &filename, T &type, const Verbosity &verbosity)
Reads data from XML file.
Definition: xml_io.cc:901
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
The Matrix class.
Definition: matpackI.h:1193
Implementation of Matrix, Vector, and such stuff.
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:628
const Numeric DEG2RAD
const Numeric PI