From 5701379d88ed30e1d21c200dc654dc0b4a761d85 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andr=C3=A9=20Schmidt?= <upwli@student.kit.edu>
Date: Wed, 17 Nov 2021 10:42:19 +0000
Subject: [PATCH] Update WMM.hpp
---
corsika/media/WMM.hpp | 60 +++++++++++++++----------------------------
1 file changed, 21 insertions(+), 39 deletions(-)
diff --git a/corsika/media/WMM.hpp b/corsika/media/WMM.hpp
index 9c4a790b9..6d15c51c1 100644
--- a/corsika/media/WMM.hpp
+++ b/corsika/media/WMM.hpp
@@ -6,27 +6,14 @@
* the license.
*/
+#include <boost/filesystem.hpp>
+#include <boost/math/special_functions/factorials.hpp>
#include <boost/math/special_functions/legendre.hpp>
+
#include <cmath>
#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/utility/CorsikaData.hpp>
#include <fstream>
-#include <iostream>
-using namespace std;
-
-double factorial(int n) {
- double factorial = 1;
- if (n < 0) {
- CORSIKA_LOG_ERROR("Factorial of a negative number does not exist.");
- abort();
- }
- else {
- for(int i = 1; i <= n; ++i) {
- factorial *= i;
- }
- }
- return factorial;
-}
namespace corsika {
MagneticFieldVector get_wmm(const CoordinateSystemPtr Cs, const double year, const LengthType altitude, const double latitude, const double longitude) {
@@ -34,13 +21,13 @@ namespace corsika {
CORSIKA_LOG_ERROR("Year has to be between 2020 and 2025.");
abort();
}
- if (altitude < -10_km || altitude > 3000_km) {
+ if (altitude < -1_km || altitude > 850_km) {
CORSIKA_LOG_WARN("Altitude should be between -10_km and 3000_km.");
}
if (latitude < -90 || latitude > 90) {
CORSIKA_LOG_ERROR("Latitude has to be between -90 and 90 degree.");
abort();
- } else if (latitude < -82 || latitude > 82) {
+ } else if (latitude < -89.992 || latitude > 89.992) {
CORSIKA_LOG_WARN("Latitude is close to the poles.");
}
if (longitude < -180 || longitude > 180) {
@@ -48,14 +35,13 @@ namespace corsika {
abort();
}
- double lat_geo = latitude * constants::pi / 180;
- double lon = longitude * constants::pi / 180;
+ const double lat_geo = latitude * constants::pi / 180;
+ const double lon = longitude * constants::pi / 180;
// Transform into spherical coordinates
- const LengthType A = 6378137_m;
const double f = 1 / 298.257223563;
const double e_squared = f * (2 - f);
- LengthType R_c = A / sqrt(1 - e_squared * pow(sin(lat_geo), 2));
+ LengthType R_c = constants::EarthRadius::Equatorial / sqrt(1 - e_squared * pow(sin(lat_geo), 2));
LengthType p = (R_c + altitude) * cos(lat_geo);
LengthType z = sin(lat_geo) * (altitude + R_c * (1 - e_squared));
LengthType r = sqrt(p * p + z * z);
@@ -63,49 +49,45 @@ namespace corsika {
const int length = 90;
double epoch, g[length], h[length], g_dot[length], h_dot[length];
- char model_name[7];
- char release_date[10];
+ std::string model_name;
+ std::string release_date;
int n[length], m[length];
- ifstream infile;
// Read in coefficients
boost::filesystem::path const path = corsika::corsika_data("GeoMag/WMM.COF");
- char mdfile[path.generic_string().size() + 1];
- strcpy(mdfile, path.generic_string().c_str()); //name and path of model file
- infile.open(mdfile);
+ boost::filesystem::ifstream file(path, std::ios::in);
// Exit if file opening failed
- if (!infile.is_open()){
+ if (!file.is_open()){
CORSIKA_LOG_ERROR("Failed opening WMM.COF.");
abort();
}
- infile >> epoch >> model_name >> release_date;
+ file >> epoch >> model_name >> release_date;
for (int i = 0; i < length; i++) {
- infile >> n[i] >> m[i] >> g[i] >> h[i] >> g_dot[i] >> h_dot[i];
+ file >> n[i] >> m[i] >> g[i] >> h[i] >> g_dot[i] >> h_dot[i];
// Time interpolation
g[i] = g[i] + (year - epoch) * g_dot[i];
h[i] = h[i] + (year - epoch) * h_dot[i];
}
- infile.close();
+ file.close();
double legendre, next_legendre, derivate_legendre;
double magneticfield[3] = {0, 0, 0};
for (int j = 0; j < length; j++) {
- // cout << assoc_legendre(n[0], m[0], sin(lat_sph)) << endl;
legendre = boost::math::legendre_p(n[j], m[j], sin(lat_sph));
next_legendre = boost::math::legendre_p(n[j] + 1, m[j], sin(lat_sph));
// Schmidt semi-normalization and Condon-Shortley phase term
if (m[j] > 0) {
- legendre *= sqrt(2 * factorial(n[j] - m[j]) / factorial(n[j] + m[j])) * pow(-1, m[j]);
- next_legendre *= sqrt(2 * factorial(n[j] + 1 - m[j]) / factorial(n[j] + 1 + m[j])) * pow(-1, m[j]);
+ legendre *= sqrt(2 * boost::math::factorial<double>(n[j] - m[j]) / boost::math::factorial<double>(n[j] + m[j])) * pow(-1, m[j]);
+ next_legendre *= sqrt(2 * boost::math::factorial<double>(n[j] + 1 - m[j]) / boost::math::factorial<double>(n[j] + 1 + m[j])) * pow(-1, m[j]);
}
derivate_legendre = (n[j] + 1) * tan(lat_sph) * legendre - sqrt(pow(n[j] + 1, 2) - pow(m[j], 2)) / cos(lat_sph) * next_legendre;
- magneticfield[0] += pow(constants::EarthRadius::Mean / r, n[j] + 2) * (g[j] * cos(m[j] * lon) + h[j] * sin(m[j] * lon)) * derivate_legendre;
- magneticfield[1] += pow(constants::EarthRadius::Mean / r, n[j] + 2) * m[j] * (g[j] * sin(m[j] * lon) - h[j] * cos(m[j] * lon)) * legendre;
- magneticfield[2] += (n[j] + 1) * pow(constants::EarthRadius::Mean / r, n[j] + 2) * (g[j] * cos(m[j] * lon) + h[j] * sin(m[j] * lon)) * legendre;
+ magneticfield[0] += pow(constants::EarthRadius::Geomagnetic_reference / r, n[j] + 2) * (g[j] * cos(m[j] * lon) + h[j] * sin(m[j] * lon)) * derivate_legendre;
+ magneticfield[1] += pow(constants::EarthRadius::Geomagnetic_reference / r, n[j] + 2) * m[j] * (g[j] * sin(m[j] * lon) - h[j] * cos(m[j] * lon)) * legendre;
+ magneticfield[2] += (n[j] + 1) * pow(constants::EarthRadius::Geomagnetic_reference / r, n[j] + 2) * (g[j] * cos(m[j] * lon) + h[j] * sin(m[j] * lon)) * legendre;
}
magneticfield[0] *= -1;
magneticfield[1] /= cos(lat_sph);
@@ -119,4 +101,4 @@ namespace corsika {
return MagneticFieldVector{Cs, magneticfield_geo[0] * 1_nT, magneticfield_geo[1] * 1_nT, magneticfield_geo[2] * -1_nT};
}
-}
\ No newline at end of file
+}
--
GitLab