00001 /* Copyright (C) 2003-2008 Nikolay Koulev <nkoulev@uni-bremen.de> 00002 00003 This program is free software; you can redistribute it and/or 00004 modify it under the terms of the GNU General Public License as 00005 published by the Free Software Foundation; either version 2 of the 00006 License, or (at your option) any later version. 00007 00008 This program is distributed in the hope that it will be useful, 00009 but WITHOUT ANY WARRANTY; without even the implied warranty of 00010 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00011 GNU General Public License for more details. 00012 00013 You should have received a copy of the GNU General Public License 00014 along with this program; if not, write to the Free Software 00015 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 00016 USA. */ 00017 00027 #include <iostream> 00028 #include <cmath> 00029 #include "arts.h" 00030 #include "matpackI.h" 00031 #include "xml_io.h" 00032 #include "legendre.h" 00033 00034 00035 extern const Numeric PI; 00036 extern const Numeric DEG2RAD; 00037 extern const Numeric EARTH_RADIUS; 00038 00039 void magfield_nk( // Output 00040 Numeric& B_r, // radial component of the geomagnetic field in [nT]. 00041 Numeric& B_th, // latitudinal component of the geomagnetic field in [nT]. 00042 Numeric& B_ph, // longitudinal component of the geomagnetic field in [nT]. 00043 00044 // Input 00045 const Numeric r, // radial distance to the point in [km] 00046 const Numeric theta, // geocentric colatitude of the point in [deg] 00047 const Numeric phi, // longitude of the point in [deg]. 00048 // All coordinates - geocentric! 00049 00050 const Index Ny // number of elapsed years after an epoch year, J - [0,4] 00051 ) 00052 00053 { 00054 00055 Numeric Phi = phi * DEG2RAD; // Longitude angle in radian. 00056 Numeric Theta = PI/2 - theta * DEG2RAD; // Colatitude angle in radian. 00057 00058 // Initializing values of the magnetic field components. 00059 B_r = 0; 00060 B_th = 0; 00061 B_ph = 0; 00062 00063 Matrix M; 00064 00065 xml_read_from_file ("geomag_coefficients.xml", M); 00066 00067 // M(i,0) and M(i,1) - the vectors with the values of the first and second coefficients 00068 // of the IGRF model. 00069 // M(i,2) and M(i,3) - the vectors with the values of the anual rate of change of the 00070 // first and second coefficient of of the IGRF model. 00071 00072 00073 // Loop over the degree number l of the Legendre polynomes. 00074 for (Index l = 1; l <= 10; l++) 00075 { 00076 // Loop over the order number m of the Legendre polynomes. 00077 for (Index m = 0; m <= l; m++) 00078 { 00079 00080 // Relating the row index in M to the coresponding 00081 // degree number l and order number m. 00082 Index j = l * (l + 1) / 2 + m - 1; 00083 00084 // Calculating the associated Schmidt quasi-normalized Legendre 00085 // polynomial for a degree number l and order number m. 00086 Numeric P_lm = 00087 g_legendre_poly_norm_schmidt (l, m, cos(Theta)); 00088 00089 // Calculating the derivative of the associated Schmidt quasi-normalized 00090 // Legendre polynomial for a degree number l and order number m. 00091 Numeric dP_lm = 00092 g_legendre_poly_norm_schmidt_deriv3 (l, m, cos(Theta)); 00093 00094 00095 // Calculating the radial (upward) component of the magnetic field. 00096 B_r += pow((Numeric)(l + 2), EARTH_RADIUS / r) * (Numeric)(l + 1) * 00097 ((M(j,0) + (Numeric)Ny * M(j,2)) * cos((Numeric)m * Phi) 00098 + (M(j,1) + (Numeric)Ny * M(j,3)) * sin((Numeric)m * Phi)) 00099 * P_lm; 00100 00101 // Calculating the latitudinal (southward) component of the magnetic field. 00102 B_th += pow(l + 2, EARTH_RADIUS / r) * 00103 ((M(j,0) + (Numeric)Ny * M(j,2)) * cos((Numeric)m * Phi) 00104 + (M(j,1) + (Numeric)Ny * M(j,3)) * sin((Numeric)m * Phi)) * 00105 dP_lm * sin(Theta); 00106 00107 00108 00109 // Calculating the longitudinal (eastward) component of the magnetic field. 00110 B_ph += pow(l + 2, EARTH_RADIUS / r) * (Numeric)m * 00111 ((M(j,0) + (Numeric)Ny * M(j,2)) * sin((Numeric)m * Phi) 00112 - (M(j,1) + (Numeric)Ny * M(j,3)) * cos((Numeric)m * Phi)) * 00113 P_lm / sin(Theta); 00114 00115 00116 } 00117 } 00118 00119 }
1.5.6