IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 61be3cbd authored by ralfulrich's avatar ralfulrich
Browse files

clang-format

parent 47c87589
No related branches found
No related tags found
1 merge request!404adding Geomagnetic Models
Hide whitespace changes
Inline Side-by-side
corsika/detail/media/WMM.inl 0 → 100644
+ 121
− 0
View file @ 61be3cbd
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
#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>
namespace corsika {
inline MagneticFieldVector get_wmm(const CoordinateSystemPtr Cs, const double year,
const LengthType altitude, const double latitude,
const double longitude) {
if (year < 2020 || year > 2025) {
CORSIKA_LOG_WARN("Year has to be between 2020 and 2025.");
}
if (altitude < -1_km || altitude > 850_km) {
CORSIKA_LOG_WARN("Altitude should be between -1_km and 850_km.");
}
if (latitude < -90 || latitude > 90) {
CORSIKA_LOG_ERROR("Latitude has to be between -90 and 90 degree.");
abort();
} else if (latitude < -89.992 || latitude > 89.992) {
CORSIKA_LOG_WARN("Latitude is close to the poles.");
}
if (longitude < -180 || longitude > 180) {
CORSIKA_LOG_WARN("Longitude should be between -180 and 180 degree.");
}
const double lat_geo = latitude * constants::pi / 180;
const double lon = longitude * constants::pi / 180;
// Transform into spherical coordinates
const double f = 1 / 298.257223563;
const double e_squared = f * (2 - f);
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);
double lat_sph = asin(z / r);
const int length = 90;
double epoch, g[length], h[length], g_dot[length], h_dot[length];
std::string model_name;
std::string release_date;
int n[length], m[length];
// Read in coefficients
boost::filesystem::path const path = corsika::corsika_data("GeoMag/WMM.COF");
boost::filesystem::ifstream file(path, std::ios::in);
// Exit if file opening failed
if (!file.is_open()) {
CORSIKA_LOG_ERROR("Failed opening WMM.COF.");
abort();
}
file >> epoch >> model_name >> release_date;
for (int i = 0; i < length; 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];
}
file.close();
double legendre, next_legendre, derivate_legendre;
double magneticfield[3] = {0, 0, 0};
for (int j = 0; j < length; j++) {
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 * 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::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);
magneticfield[2] *= -1;
// Transform back into geodetic coordinates
double magneticfield_geo[3];
magneticfield_geo[0] = magneticfield[0] * cos(lat_sph - lat_geo) -
magneticfield[2] * sin(lat_sph - lat_geo);
magneticfield_geo[1] = magneticfield[1];
magneticfield_geo[2] = magneticfield[0] * sin(lat_sph - lat_geo) +
magneticfield[2] * cos(lat_sph - lat_geo);
return MagneticFieldVector{Cs, magneticfield_geo[0] * 1_nT,
magneticfield_geo[1] * 1_nT, magneticfield_geo[2] * -1_nT};
}
} // namespace corsika
\ No newline at end of file
corsika/media/WMM.hpp
+ 6
− 4
View file @ 61be3cbd
...@@ -25,14 +25,16 @@ namespace corsika { ...@@ -25,14 +25,16 @@ namespace corsika {
* @param latitude Latitude of the location to evaluate the field at, * @param latitude Latitude of the location to evaluate the field at,
in degrees between -90 and 90 (negative for southern hemisphere). in degrees between -90 and 90 (negative for southern hemisphere).
* @param longitute Longitude of the location to evaluate the field at, * @param longitute Longitude of the location to evaluate the field at,
in degrees between -180 and 180 (negative for western hemisphere). in degrees between -180 and 180 (negative for western
hemisphere).
* *
* @returns The magnetic field vector in nT. * @returns The magnetic field vector in nT.
* *
*/ */
inline MagneticFieldVector get_wmm(const CoordinateSystemPtr Cs, const double year, inline MagneticFieldVector get_wmm(const CoordinateSystemPtr Cs, const double year,
const LengthType altitude, const double latitude, const double longitude); const LengthType altitude, const double latitude,
const double longitude);
} // namespace corsika } // namespace corsika
#include <corsika/detail/media/WMM.inl> #include <corsika/detail/media/WMM.inl>
tests/media/testMagneticField.cpp
+ 58
− 68
View file @ 61be3cbd
...@@ -30,72 +30,62 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") { ...@@ -30,72 +30,62 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
CoordinateSystemPtr const& gCS = get_root_CoordinateSystem(); CoordinateSystemPtr const& gCS = get_root_CoordinateSystem();
Point const gOrigin(gCS, {0_m, 0_m, 0_m}); Point const gOrigin(gCS, {0_m, 0_m, 0_m});
// setup our interface types SECTION("UniformMagneticField interface") {
using IModelInterface = IMagneticFieldModel<IMediumModel>;
using AtmModel = UniformMagneticField<HomogeneousMedium<IModelInterface>>; // setup our interface types
using IModelInterface = IMagneticFieldModel<IMediumModel>;
// the composition we use for the homogenous medium using AtmModel = UniformMagneticField<HomogeneousMedium<IModelInterface>>;
NuclearComposition const protonComposition({Code::Proton}, {1.});
// the composition we use for the homogenous medium
// create a magnetic field vector NuclearComposition const protonComposition({Code::Proton}, {1.});
Vector B0(gCS, 0_T, 0_T, 0_T);
// create a magnetic field vector
// the constant density Vector B0(gCS, 0_T, 0_T, 0_T);
const auto density{19.2_g / cube(1_cm)};
// the constant density
// create our atmospheric model const auto density{19.2_g / cube(1_cm)};
AtmModel medium(B0, density, protonComposition);
// create our atmospheric model
// and test at several locations AtmModel medium(B0, density, protonComposition);
CHECK(B0.getComponents(gCS) ==
medium.getMagneticField(Point(gCS, -10_m, 4_m, 35_km)).getComponents(gCS)); // and test at several locations
CHECK( CHECK(B0.getComponents(gCS) ==
B0.getComponents(gCS) == medium.getMagneticField(Point(gCS, -10_m, 4_m, 35_km)).getComponents(gCS));
medium.getMagneticField(Point(gCS, 1000_km, -1000_km, 1000_km)).getComponents(gCS)); CHECK(B0.getComponents(gCS) ==
CHECK(B0.getComponents(gCS) == medium.getMagneticField(Point(gCS, 1000_km, -1000_km, 1000_km))
medium.getMagneticField(Point(gCS, 0_m, 0_m, 0_m)).getComponents(gCS)); .getComponents(gCS));
CHECK(B0.getComponents(gCS) ==
// create a new magnetic field vector medium.getMagneticField(Point(gCS, 0_m, 0_m, 0_m)).getComponents(gCS));
Vector B1(gCS, 23_T, 57_T, -4_T);
// create a new magnetic field vector
// and update this atmospheric model Vector B1(gCS, 23_T, 57_T, -4_T);
medium.setMagneticField(B1);
// and update this atmospheric model
// and test at several locations medium.setMagneticField(B1);
CHECK(B1.getComponents(gCS) ==
medium.getMagneticField(Point(gCS, -10_m, 4_m, 35_km)).getComponents(gCS)); // and test at several locations
CHECK( CHECK(B1.getComponents(gCS) ==
B1.getComponents(gCS) == medium.getMagneticField(Point(gCS, -10_m, 4_m, 35_km)).getComponents(gCS));
medium.getMagneticField(Point(gCS, 1000_km, -1000_km, 1000_km)).getComponents(gCS)); CHECK(B1.getComponents(gCS) ==
CHECK(B1.getComponents(gCS) == medium.getMagneticField(Point(gCS, 1000_km, -1000_km, 1000_km))
medium.getMagneticField(Point(gCS, 0_m, 0_m, 0_m)).getComponents(gCS)); .getComponents(gCS));
CHECK(B1.getComponents(gCS) ==
// check the density and nuclear composition medium.getMagneticField(Point(gCS, 0_m, 0_m, 0_m)).getComponents(gCS));
CHECK(density == medium.getMassDensity(Point(gCS, 0_m, 0_m, 0_m))); }
medium.getNuclearComposition();
SECTION("WorldMagneticModel") {
// create a line of length 1 m // create earth magnetic field vector
Line const line(gOrigin, Vector<SpeedType::dimension_type>( MagneticFieldVector Earth_B_1 = get_wmm(gCS, 2022.5, 100_km, -80, -120);
gCS, {1_m / second, 0_m / second, 0_m / second})); MagneticFieldVector Earth_B_2 = get_wmm(gCS, 2022.5, 0_km, 0, 120);
MagneticFieldVector Earth_B_3 = get_wmm(gCS, 2020, 100_km, 80, 0);
// the end time of our line CHECK(Earth_B_1.getX(gCS) / 1_nT == Approx(5815).margin(0.03));
auto const tEnd = 1_s; CHECK(Earth_B_1.getY(gCS) / 1_nT == Approx(14803).margin(0.03));
CHECK(Earth_B_1.getZ(gCS) / 1_nT == Approx(49755.3).margin(0.03));
// and the associated trajectory CHECK(Earth_B_2.getX(gCS) / 1_nT == Approx(39684.7).margin(0.03));
setup::Trajectory const track = CHECK(Earth_B_2.getY(gCS) / 1_nT == Approx(-42.2).margin(0.03));
setup::testing::make_track<setup::Trajectory>(line, tEnd); CHECK(Earth_B_2.getZ(gCS) / 1_nT == Approx(10809.5).margin(0.03));
CHECK(Earth_B_3.getX(gCS) / 1_nT == Approx(6261.8).margin(0.03));
// create earth magnetic field vector CHECK(Earth_B_3.getY(gCS) / 1_nT == Approx(-185.5).margin(0.03));
MagneticFieldVector Earth_B_1 = get_wmm(gCS, 2022.5, 100_km, -80, -120); CHECK(Earth_B_3.getZ(gCS) / 1_nT == Approx(-52429.1).margin(0.03));
MagneticFieldVector Earth_B_2 = get_wmm(gCS, 2022.5, 0_km, 0, 120); }
MagneticFieldVector Earth_B_3 = get_wmm(gCS, 2020, 100_km, 80, 0);
CHECK(Earth_B_1.getX(gCS) / 1_nT == Approx(5815).margin(0.03));
CHECK(Earth_B_1.getY(gCS) / 1_nT == Approx(14803).margin(0.03));
CHECK(Earth_B_1.getZ(gCS) / 1_nT == Approx(49755.3).margin(0.03));
CHECK(Earth_B_2.getX(gCS) / 1_nT == Approx(39684.7).margin(0.03));
CHECK(Earth_B_2.getY(gCS) / 1_nT == Approx(-42.2).margin(0.03));
CHECK(Earth_B_2.getZ(gCS) / 1_nT == Approx(10809.5).margin(0.03));
CHECK(Earth_B_3.getX(gCS) / 1_nT == Approx(6261.8).margin(0.03));
CHECK(Earth_B_3.getY(gCS) / 1_nT == Approx(-185.5).margin(0.03));
CHECK(Earth_B_3.getZ(gCS) / 1_nT == Approx(-52429.1).margin(0.03));
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum