CDPP / AMDA_Kernel

Commit 38bebd9ceb82823a1c280f35d95664be1b691c49

Authored by Hacene SI HADJ MOHAND
1 parent 3699433c
Exists in master and in 83 other branches 10038, 10098, 5390, 5770, 6728, 7616, 7811, 7875, 8034_7628, 8038, 8041, 8045_8044, 8295, 8459, 8528, 8928, 9016, 9016Bis, 9017, 9017-new, 9049, 9084, 9259, 9326, 9347, 9383, 9719, 9719.new, AllLegends, CDF_ISTP, CDF_ISTP_HOPEFULLY_ALL_WORKING, FER_10400, FER_10412, FER_10663, FER_10765, FER_11000, FER_11100, FER_11242, FER_11277, FER_11500, FER_11537, FER_11588, FER_11651, FER_19836, FER_5344, FER_6031, FER_6598, FER_6966, FER_7318, FER_7618, FER_7857, FER_8455, FER_9629, FER_CR41, FER_CR42, FER_CR43, FER_TFCat, Histo1D2D, JUNO_JEDI, MAZ_10936, MAZ_6018, MAZ_9016_New, MAZ_FIX_9016, SpeasyGet, amdadev, amdadev_test, cu2022, das2, density_maps, develop, evol_9042, expression-parser-bre, feature/9259, frames-transformation, pickle, rm_6912, rm_7926, rm_8050, rm_8932, shue, shue97, test-merge, ubuntu

lib Jupeter field added

Showing 10 changed files with 1194 additions and 1 deletions   Show diff stats
Inline Side-by-side
CMakeLists.txt
  Diff comments   View file @38bebd9
... ... @@ -145,6 +145,7 @@ add_subdirectory(src/ParamGetImpl/ConstantInterface)
145 145 add_subdirectory(src/InternLib)
146 146 add_subdirectory(src/ExternLib/Deriv)
147 147 add_subdirectory(src/ExternLib/Morschhauser)
  148 +add_subdirectory(src/ExternLib/Jupeter_JRM09_CAN_81)
148 149 add_subdirectory(src/ExternLib/Ram_Presure)
149 150 add_subdirectory(src/ExternLib/TimeShifted)
150 151 add_subdirectory(src/ExternLib/GetClbInfo)
... ...
src/ExternLib/Jupeter_JRM09_CAN_81/AMDAPlugin.cc 0 → 100644
  Diff comments   View file @38bebd9
... ... @@ -0,0 +1,38 @@
  1 +/*
  2 + * To change this license header, choose License Headers in Project Properties.
  3 + * To change this template file, choose Tools | Templates
  4 + * and open the template in the editor.
  5 + */
  6 +
  7 +/*
  8 + * File: AMDAPlugin.cc
  9 + * Author: hacene
  10 + *
  11 + * Created on June 5, 2020, 2:31 PM
  12 + */
  13 +
  14 +#include "Jupeter_JRM09_CAN81_MagProcess.hh"
  15 +#include "ServicesServer.hh"
  16 +#include "PluginManager.hh"
  17 +
  18 +using namespace AMDA::Parameters;
  19 +
  20 +/**
  21 + Retrieve the Plugin version we're going to expect
  22 +*/
  23 +extern "C" const char* getPluginVersion()
  24 +{
  25 + return "(Version)";
  26 +}
  27 +
  28 +/**
  29 + Tells us to register our functionality to an engine kernel
  30 +*/
  31 +extern "C" void registerPlugin(AMDA::Plugins::PluginManager & pm)
  32 +{
  33 + ProcessFactory factJupeter_JRM09_CAN81_MagProcess = boost::factory<Jupeter_JRM09_CAN81_MagProcess*>();
  34 + ServicesServer::getInstance()->addProcessFactory("jupeter_JRM09_CAN81_bmag", factJupeter_JRM09_CAN81_MagProcess);
  35 + ServicesServer::getInstance()->linkProcessWithPlugin("jupeter_JRM09_CAN81_bmag", pm.getCurrentPluginPath());
  36 +}
  37 +
  38 +
... ...
src/ExternLib/Jupeter_JRM09_CAN_81/CMakeLists.txt 0 → 100644
  Diff comments   View file @38bebd9
... ... @@ -0,0 +1,31 @@
  1 +PROJECT(Jupeter_JRM09_CAN_81)
  2 +
  3 +set(LIBRARY_OUTPUT_PATH ${PLUGIN_OUTPUT_PATH})
  4 +
  5 +include_directories(
  6 + ${CMAKE_HOME_DIRECTORY}/src/InternLib/
  7 + ${CMAKE_HOME_DIRECTORY}/src/Common/
  8 + ${CMAKE_HOME_DIRECTORY}/src/Parameters/
  9 + ${CMAKE_HOME_DIRECTORY}/src/Plugins/
  10 + ${CMAKE_HOME_DIRECTORY}/src/helpers/
  11 + ${LOG4CXX_INCLUDE_DIR}
  12 + ${Boost_INCLUDE_DIR}
  13 +)
  14 +
  15 +#Library configuration
  16 +file(
  17 + GLOB_RECURSE
  18 + source_files
  19 + ./*
  20 +)
  21 +
  22 +ADD_LIBRARY( Jupeter_JRM09_CAN_81Process SHARED ${source_files} )
  23 +
  24 +target_link_libraries(
  25 + MorschhauserProcess
  26 + ${LOG4CXX_LIBRARIES}
  27 + Parameters
  28 + InternLib
  29 + TimeTableCatalog
  30 + SpiceKernel
  31 +)
0 32 \ No newline at end of file
... ...
src/ExternLib/Jupeter_JRM09_CAN_81/Jupeter_JRM09_CAN81.hh 0 → 100644
  Diff comments   View file @38bebd9
... ... @@ -0,0 +1,312 @@
  1 +/*
  2 + * To change this license header, choose License Headers in Project Properties.
  3 + * To change this template file, choose Tools | Templates
  4 + * and open the template in the editor.
  5 + */
  6 +
  7 +/*
  8 + * File: getJupeterMag.hh
  9 + * Author: hacene
  10 + *
  11 + * Created on May 29, 2020, 2:50 PM
  12 + */
  13 +
  14 +#ifndef GETJUPETERMAG_HH
  15 +#define GETJUPETERMAG_HH
  16 +
  17 +#include "Jupeter_JRM09_CAN81_constants.hh"
  18 +#include "Jupeter_JRM09_CAN81_tools.hh"
  19 +
  20 +#include "Parameter.hh"
  21 +#include "ParamData.hh"
  22 +#include "DataTypeMath.hh"
  23 +#include "Operation.hh"
  24 +#include "SpiceKernelMgr.hh"
  25 +
  26 +namespace AMDA {
  27 + namespace Parameters{
  28 + using namespace GETJUPETERMAG_CONSTANTS;
  29 +
  30 + /**
  31 + * @class GetJupeterMagCommon
  32 + * @brief It is responsible to compute Jupeter Magnetic Field along an orbit. Abstract class
  33 + * @details This class implement the interface Operation. Input vector are [R, Z]
  34 + */
  35 +
  36 + template<typename DataType, class TOutputParamData>
  37 + class Jupeter_JRM09_CAN81 : public Operation {
  38 + public:
  39 +
  40 + /**
  41 + * @brief Constructor.
  42 + * @details Create the ParamData type of the input ParamData.
  43 + */
  44 + Jupeter_JRM09_CAN81(Process& pProcess, ParamDataSpec<std::vector<DataType> >&paramInput, can) : Operation(pProcess),
  45 + _paramInput(paramInput),
  46 + _paramOutput(new TOutputParamData()), _can(can) {
  47 + _paramDataOutput = _paramOutput;
  48 + }
  49 +
  50 + virtual ~Jupeter_JRM09_CAN81() {
  51 + }
  52 +
  53 + /**
  54 + * @overload Operation::write(ParamDataIndexInfo &pParamDataIndexInfo)
  55 + */
  56 + void write(ParamDataIndexInfo &pParamDataIndexInfo) {
  57 + for (unsigned int _index = pParamDataIndexInfo._startIndex;
  58 + _index < pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess;
  59 + ++_index) {
  60 + std::vector<DataType> lVal = _paramInput.getDataList()[_index];
  61 + double crtTime = _paramInput.getTime(_index);
  62 +
  63 + std::vector<DataType> input;
  64 + input.resize(3);
  65 +
  66 + short sign = (lVal[0] > 0) ? 1 : 0;
  67 +
  68 + std::vector<DataType> outputB;
  69 + outputB.resize(3);
  70 +
  71 + std::vector<DataType>outputBxyz
  72 + outputBxyz.resize(3);
  73 +
  74 + DataType outputBm;
  75 +
  76 + // GetJupiterMag00(std::vector<DataType> inputSphCoord, bool CAN, std::vector<DataType>& outputBrtp, DataType& outputBm, std::vector<DataType>& outputBxyz)
  77 + GetJupiterMag00(input, _can, outputBrtp, outputBm, outputBxyz);
  78 + _paramOutput->pushTime(crtTime);
  79 + pushResult(outputBrtp, outputBm, outputBxyz);
  80 + }
  81 + }
  82 + protected:
  83 +
  84 + virtual void pushResult(std::vector<DataType> fieldSph, DataType magnitude_, std::vector<DataType> fieldCoor) = 0;
  85 + private:
  86 + bool _can;
  87 +
  88 + int can81(std::vector<DataType> inputCoor, std::vector<DataType> outpField) {
  89 + /**
  90 + * @brief Jupiter's external magnetic field model:
  91 + * @details ; The Connerney+1981_JGR model: Modeling the Jovian Current Sheet and Inner Magnetophere.
  92 + *The CAN model should be only used near the Jupiter (<~30 RJ).
  93 + *Input: rho and z in cylindrical coords, z is along the magnetic axis.
  94 + *Output: Brho and Bz (nT) also in cylindrical coords. Note that the model is axis-symmetric.
  95 + */
  96 +
  97 + // [constants]
  98 + double D = 2.5; // Rj
  99 + double a = 5.0; // Rj
  100 + double b = 50.0; // 70.0 ; Rj ; 50.0 in Connerney1981
  101 + double mu0_I0 = 450.0; //; 175.0*2 ; 225.0*2 in Connerney1981
  102 + double minval = 1.0e-15;
  103 +
  104 + double rho = inputCoor[0];
  105 + double z = inputCoor[1];
  106 + double F1, F2, F3, Brhe, tmp, Bz, Brho, Brho2, Bz2;
  107 +
  108 + // [Approximate formulas given in Connerney+1981_The magnetic field in Jupiter.]
  109 + // [Cylindrical coordinates: in nT]
  110 +
  111 + if (rho < a) {
  112 + F1 = std::sqrt((z - D)*(z - D)* + a * a);
  113 + F2 = std::sqrt((z + D)*(z + D) + a * a);
  114 + F3 = std::sqrt(z * z + a * a);
  115 +
  116 + Brho = 0.5 * rho * (1.0 / F1 - 1.0 / F2);
  117 + tmp = (z - D) / (F1 * F1 * F1) - (z + D) / (F2 * F2 * F2);
  118 + Bz = 2.0 * D / F3 - 0.25 * rho * rho*tmp;
  119 + } else {
  120 + F1 = std::sqrt((z - D)*(z - D)* + rho * rho);
  121 + F2 = std::sqrt((z + D)*(z + D) + rho * rho);
  122 + F3 = std::sqrt(z * z + rho * rho);
  123 +
  124 + Brho = (F1 - F2 + 2 * D) / rho;
  125 + if (std::abs(z) > D && z < 0) Brho = (F1 - F2 - 2 * D) / rho;
  126 + if (std::abs(z) < D) Brho = (F1 - F2 + 2 * z) / rho;
  127 + Brho -= 0.25 * a * a * rho * (1.0 / (F1 * F1 * F1) - 1.0 / (F2 * F2 * F2));
  128 + tmp = (z - D) / (F1 * F1 * F1) - (z + D) / (F2 * F2 * F2);
  129 + Bz = 2.0 * D / F3 - 0.25 * a * a*tmp;
  130 + }
  131 + F1 = std::sqrt((z - D)*(z - D) + b * b);
  132 + F2 = std::sqrt((z + D)*(z + D) + b * b);
  133 + F3 = std::sqrt(z * z + b * b);
  134 + Brho2 = 0.5 * rho * (1 / F1 - 1 / F2);
  135 + Bz2 = 2 * D / F3 - 0.25 * rho * rho * ((z - D) / (F1 * F1 * F1) - (z + D) / (F2 * F2 * F2));
  136 +
  137 + Brho -= Brho2;
  138 + Bz -= Bz2;
  139 +
  140 + Brho *= 0.5 * mu0_I0;
  141 + Bz *= 0.5 * mu0_I0;
  142 + outpField[0] = Brho;
  143 + outpField[1] = Brho;
  144 +
  145 + return 0;
  146 + };
  147 +
  148 + /**
  149 + * GetJupiterMag00.pro
  150 + =================================================
  151 + Description:
  152 + This program is to calculate the Jovian magnetic field of JRM09 (Connerney et al., 2018) and
  153 + CAN (Connerney et al., 1981) model.
  154 +
  155 + All of the relevant subroutines are included in this single file, including the Schmdt coefficients
  156 + needed in JRM09 model.
  157 +
  158 + (1) Jupiter's Magnetic field Model: JRM09 (+CAN), CAN.
  159 + (2) Refer to: Connerney+2018_GRL JRM09 model
  160 + Connerney+1981 CAN model
  161 + (3) All the coordinates used here are Right-Hand-Side, because of SPICE convention.
  162 + (4) The positions can be in Arrays, i.e., multi-positions at one time, which can be more efficient.
  163 +
  164 + Input:
  165 + r (Rj), theta (deg), phi (deg), spherical coordinates in System III RHS (i.e., IAU_Jupiter), better use double!
  166 +
  167 + optional:
  168 + /CAN include the CAN model or not (true = yes)
  169 + Output:
  170 + Br (nT), Bt (nT), Bp (nT), Bmag (nT) in System III RHS.
  171 + Bxyz = Bxyz optional output Bxyz in System III RHS
  172 +
  173 +Example:
  174 + GetJupiterMag00, [20.0d0,20.00000001d0], [90.0d0,90.000001d0], [120.0d0, 120.0d0], $
  175 + Br, Bt, Bp, Bmag1, Bxyz = Bxyz, /can
  176 + GetJupiterMag00, 20.0d0, 90.000001d0, 120.0d0, Br, Bt, Bp, Bmag2, Bxyz = Bxyz, /can
  177 +
  178 + History:
  179 + (1) writen by Yuxian Wang 2020-01-16 02:26:49
  180 + ywang@irap.omp.eu + yxwang@spaceweather.ac.cn
  181 + =================================================
  182 + */
  183 +
  184 + int GetJupiterMag00(std::vector<DataType> inputSphCoord, bool CAN, std::vector<DataType>& outputBrtp, DataType& outputBm, std::vector<DataType>& outputBxyz) {
  185 + // [position in: ['Rj, deg, deg'] to ['Rj, rad,rad']]
  186 + const long double PI = std::atan(1)*4.0;
  187 + long double r, theta, phi, x, radialfactor;
  188 + r = inputCoord[0];
  189 + theta = inputCoord[1] * PI / 180;
  190 + phi = inputCoord[2] * PI / 180;
  191 +
  192 + // [try to avoid invoke restore one more time.]
  193 + x = std::cos(theta);
  194 + int maxN = 10 // JRM09 degrees
  195 + std::vector<std::vector<long double>> Plm;
  196 + std::vector<std::vector<long double>> dPlm;
  197 +
  198 + //int AllLegendre(T1 x, int lma, int mmax, std::vector<std::vector<T2>> Plm,
  199 + //std::vector<std::vector<T2>> dPlm, bool isSineTheta, int SineTheta, bool DoDerivative , int normalize = 1)
  200 + int res = AMDA::GETJUPETERMAG_TOOLS::Jupeter_JRM09_CAN81_tools::AllLegendre(x, maxN, maxN, Plm, dPlm, true, 1, true);
  201 + // [Calculate Components of B]
  202 +
  203 + for (int l = 1; l < maxN; l++) {
  204 + radialfactor = std::pow(r, (-l + 2.0));
  205 + for (int m = 0; m < l; l++) {
  206 + //-------------find radial b field component---------------
  207 + long double gsmdtlm = GSMDT[l][m];
  208 + long double hsmdtlm = HSMDT[l][m];
  209 + long double phifactor = gsmdtlm * std::cos(m * phi) + hsmdtlm * std::sin(m * phi);
  210 + long double Pnm = Plm[l][m];
  211 + outputBrtp[0] += (l + 1) * radialfactor * phifactor * Pnm;
  212 +
  213 + //--------------find theta b field component---------------
  214 + long double dpnmdtheta = dPlm[l][m]; // dpdtheta(L, M, theta, /inrad)
  215 + outputBrtp[1] -= radialfactor * phifactor * dpnmdtheta;
  216 +
  217 + //-------------------------------------------------------
  218 + outputBrtp[2] += m * radialfactor * (gsmdtlm * sin(m * phi) - hsmdtlm * cos(m * phi)) * pnm;
  219 + }
  220 + }
  221 + // --------correct bphi and don't divide by 0!------------
  222 + long double sintheta = sin(outputBrtp[1]);
  223 + if (sintheta != 0) {
  224 + outputBrtp[2] /= sintheta;
  225 + } else {
  226 + outputBrtp[2] = 0.0;
  227 + }
  228 + // [Add CAN model]
  229 + if (CAN) {
  230 + // [IAU -> JSM xyz]
  231 + std::vector<DataType> IAUxyz = AMDA::GETJUPETERMAG_TOOLS::Jupeter_JRM09_CAN81_tools::Sph2CarP(inputCoord);
  232 + std::vector<DataType>JSMxyz = AMDA::GETJUPETERMAG_TOOLS::Jupeter_JRM09_CAN81_tools::Coord_SIII2JSM(IAUxyz);
  233 + // [JSM: xyz -> rho z]
  234 + long double rho = sqrt(JSMxyz[0] * JSMxyz[0] + JSMxyz[1] * JSMxyz[1]);
  235 + long double zzz = JSMxyz[2];
  236 + std::vector<long double> coorIn{rho, zzz};
  237 + std::vector<long double> Bcan; // [rho z, point]
  238 + int resCan81 = CAN81(coorIn, Bcan);
  239 + // [B_JSM: rho z -> xyz]
  240 + phi = std::atan(JSMxyz[1] / JSMxyz[0]);
  241 + if (JSMxyz[0] < 0)
  242 + phi += PI;
  243 + if (phi < 0)
  244 + phi += 2.0 * PI;
  245 + std::vector<long double> BJSMxyz{Bcan[0] * std::cos(phi), Bcan[0] * std::sin(phi), Bcan[1]};
  246 + // [JSM to System III]
  247 + std::vector<long double> BIAUxyz = AMDA::GETJUPETERMAG_TOOLS::Jupeter_JRM09_CAN81_tools::Coord_JSM2SIII(BJSMxyz);
  248 + std::vector<long double> BIAUrtp = AMDA::GETJUPETERMAG_TOOLS::Jupeter_JRM09_CAN81_tools::Car2SphV(IAUxyz, BIAUxyz);
  249 + outputBrtp[0] += BIAUrtp[0];
  250 + outputBrtp[1] += BIAUrtp[1];
  251 + outputBrtp[2] += BIAUrtp[2];
  252 + }
  253 + // [Bx, By, Bz: Bxyz in System III RHS]
  254 + Bmag = std::sqrt(outputBrtp[0] * outputBrtp[0] + outputBrtp[1] * outputBrtp[1] + outputBrtp[2] * outputBrtp[2]);
  255 + Bxyz = Sph2CarV(inputSphCoord, outputBrtp);
  256 + return 0;
  257 + };
  258 +
  259 +
  260 + };
  261 + template <typename DataType>
  262 + class Jupeter_JRM09_CAN81BMag:public Jupeter_JRM09_CAN81<DataType, ParamDataSpec<DataType>> {
  263 + public:
  264 + /**
  265 + * @brief Constructor.
  266 + * @details Create the ParamData type of the input ParamData.
  267 + */
  268 + Jupeter_JRM09_CAN81BMag(Process& pProcess, ParamDataSpec<std::vector<DataType> >& paramInput, bool can) :
  269 + Jupeter_JRM09_CAN81<DataType, ParamDataSpec<DataType>>(pProcess, paramInput, can) {
  270 + }
  271 + protected:
  272 + virtual void pushResult(std::vector<DataType> /*fieldSph_*/, DataType magnitude_, std::vector<DataType> /*fieldCart_*/) {
  273 + Jupeter_JRM09_CAN81BMag<DataType, ParamDataSpec>>::_paramOutput->getDataList().push_back(magnitude_);
  274 + }
  275 + };
  276 +
  277 + template <typename DataType>
  278 + class Jupeter_JRM09_CAN81BSphr:public Jupeter_JRM09_CAN81<DataType, ParamDataSpec<DataType>> {
  279 + public:
  280 + /**
  281 + * @brief Constructor.
  282 + * @details Create the ParamData type of the input ParamData.
  283 + */
  284 + Jupeter_JRM09_CAN81Sphr(Process& pProcess, ParamDataSpec<std::vector<DataType> >& paramInput, bool can) :
  285 + Jupeter_JRM09_CAN81<DataType, ParamDataSpec<DataType>>(pProcess, paramInput, can) {
  286 + }
  287 + protected:
  288 + virtual void pushResult(std::vector<DataType> fieldSph_, DataType/*magnitude_*/, std::vector<DataType> /*fieldCart_*/) {
  289 + Jupeter_JRM09_CAN81Sphr<DataType, ParamDataSpec>>::_paramOutput->getDataList().push_back(fieldSph);
  290 + }
  291 + };
  292 +
  293 + template <typename DataType>
  294 + class Jupeter_JRM09_CAN81BCart:public Jupeter_JRM09_CAN81<DataType, ParamDataSpec<DataType>> {
  295 + public:
  296 + /**
  297 + * @brief Constructor.
  298 + * @details Create the ParamData type of the input ParamData.
  299 + */
  300 + Jupeter_JRM09_CAN81Cart(Process& pProcess, ParamDataSpec<std::vector<DataType> >& paramInput, bool can) :
  301 + Jupeter_JRM09_CAN81<DataType, ParamDataSpec<DataType>>(pProcess, paramInput, can) {
  302 + }
  303 + protected:
  304 + virtual void pushResult(std::vector<DataType> /*fieldSph_*/, DataType/*magnitude_*/, std::vector<DataType>& fieldCart_) {
  305 + Jupeter_JRM09_CAN81Cart<DataType, ParamDataSpec>>::_paramOutput->getDataList().push_back(fieldSph);
  306 + }
  307 + };
  308 + }
  309 +} // end AMDA
  310 +
  311 +#endif /* GETJUPETERMAG_HH */
  312 +
... ...
src/ExternLib/Jupeter_JRM09_CAN_81/Jupeter_JRM09_CAN81Creator.hh 0 → 100644
  Diff comments   View file @38bebd9
... ... @@ -0,0 +1,209 @@
  1 +
  2 +/*
  3 + * To change this license header, choose License Headers in Project Properties.
  4 + * To change this template file, choose Tools | Templates
  5 + * and open the template in the editor.
  6 + */
  7 +
  8 +
  9 +/*
  10 + * File: Jupeter_JRM09_CAN81Creator.hh
  11 + * Author: hacene
  12 + *
  13 + * Created on June 5, 2020, 11:01 AM
  14 + */
  15 +
  16 +#ifndef JUPETER_JRM09_CAN81CREATOR_HH
  17 +#define JUPETER_JRM09_CAN81CREATOR_HH
  18 +
  19 +#include "ParamData.hh"
  20 +#include "VisitorOfParamData.hh"
  21 +#include "Jupeter_JRM09_CAN81.hh"
  22 +
  23 +namespace AMDA {
  24 +namespace Parameters {
  25 +
  26 +/**
  27 + @class Jupeter_RM09_CAN81Creator
  28 + *@brief Creator of the Operation used to compute Jupeter_RM09_CAN81
  29 + */
  30 +
  31 +class Jupeter_RM09_CAN81Creator : public VisitorOfParamData {
  32 +public:
  33 +
  34 + Jupeter_RM09_CAN81Creator(Process& pProcess, ParamData& paramInput, std::string output, bool can)
  35 + : _process(pProcess), _paramData(paramInput), _operation(NULL), _output(output) {
  36 +
  37 + _paramData.accept(*this);
  38 + }
  39 +
  40 + /**
  41 + * @overload VisitorOfParamData::visit(ParamDataScalaireShort *)
  42 + */
  43 + void visit(ParamDataScalaireShort *) {
  44 + BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("ParamDataScalaireShort data not supported"));
  45 + }
  46 +
  47 + /**
  48 + * @overload VisitorOfParamData::visit(ParamDataScalaireFloat *)
  49 + */
  50 + void visit(ParamDataScalaireFloat *) {
  51 + BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("ParamDataScalaireFloat data not supported"));
  52 + }
  53 +
  54 + /**
  55 + * @overload VisitorOfParamData::visit(ParamDataScalaireDouble *)
  56 + */
  57 + void visit(ParamDataScalaireDouble *) {
  58 + BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("ParamDataScalaireDouble data not supported"));
  59 + }
  60 +
  61 + /**
  62 + * @overload VisitorOfParamData::visit(ParamDataScalaireLongDouble *)
  63 + */
  64 + void visit(ParamDataScalaireLongDouble *) {
  65 + BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("ParamDataScalaireLongDouble data not supported"));
  66 + }
  67 +
  68 + /**
  69 + * @overload VisitorOfParamData::visit(ParamDataScalaireInt *)
  70 + */
  71 + void visit(ParamDataScalaireInt *) {
  72 + BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("ParamDataScalaireInt data not supported"));
  73 + }
  74 +
  75 + /**
  76 + * @overload VisitorOfParamData::visit(ParamDataLogicalData *)
  77 + */
  78 + void visit(ParamDataLogicalData *) {
  79 + BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("ParamDataLogicalData data not supported"));
  80 + }
  81 +
  82 + /**
  83 + * @overload VisitorOfParamData::visit(ParamDataTab1DShort *)
  84 + */
  85 + void visit(ParamDataTab1DShort *) {
  86 + if (_output == "sphr")
  87 + _operation = new Jupeter_RM09_CAN81Sphr<short>(_process, dynamic_cast<ParamDataTab1DShort &> (_paramData));
  88 + else if (_output == "mag")
  89 + _operation = new Jupeter_RM09_CAN81BMag<short>(_process, dynamic_cast<ParamDataTab1DShort &> (_paramData));
  90 + else
  91 + _operation = new Jupeter_RM09_CAN81Cart<short>(_process, dynamic_cast<ParamDataTab1DShort &> (_paramData));
  92 + }
  93 +
  94 + /**
  95 + * @overload VisitorOfParamData::visit(ParamDataTab1DFloat *)
  96 + */
  97 + void visit(ParamDataTab1DFloat *) {
  98 + if (_output == "sphr")
  99 + _operation = new Jupeter_RM09_CAN81Sphr<short>(_process, dynamic_cast<ParamDataTab1DShort &> (_paramData));
  100 + else if (_output == "mag")
  101 + _operation = new Jupeter_RM09_CAN81BMag<short>(_process, dynamic_cast<ParamDataTab1DShort &> (_paramData));
  102 + else
  103 + _operation = new Jupeter_RM09_CAN81Cart<short>(_process, dynamic_cast<ParamDataTab1DShort &> (_paramData));
  104 + }
  105 +
  106 + /**
  107 + * @overload VisitorOfParamData::visit(ParamDataTab1DDouble *)
  108 + */
  109 + void visit(ParamDataTab1DDouble *) {
  110 + if (_output == "sphr")
  111 + _operation = new Jupeter_RM09_CAN81Sphr<short>(_process, dynamic_cast<ParamDataTab1DShort &> (_paramData));
  112 + else if (_output == "mag")
  113 + _operation = new Jupeter_RM09_CAN81BMag<short>(_process, dynamic_cast<ParamDataTab1DShort &> (_paramData));
  114 + else
  115 + _operation = new Jupeter_RM09_CAN81Cart<short>(_process, dynamic_cast<ParamDataTab1DShort &> (_paramData));
  116 + }
  117 +
  118 + /**
  119 + * @overload VisitorOfParamData::visit(ParamDataTab1DLongDouble *)
  120 + */
  121 + void visit(ParamDataTab1DLongDouble *) {
  122 + if (_output == "sphr")
  123 + _operation = new Jupeter_RM09_CAN81Sphr<short>(_process, dynamic_cast<ParamDataTab1DShort &> (_paramData));
  124 + else if (_output == "mag")
  125 + _operation = new Jupeter_RM09_CAN81BMag<short>(_process, dynamic_cast<ParamDataTab1DShort &> (_paramData));
  126 + else
  127 + _operation = new Jupeter_RM09_CAN81Cart<short>(_process, dynamic_cast<ParamDataTab1DShort &> (_paramData));
  128 + }
  129 +
  130 + /**
  131 + * @overload VisitorOfParamData::visit(ParamDataTab1DInt *)
  132 + */
  133 + void visit(ParamDataTab1DInt *) {
  134 + if (_output == "sphr")
  135 + _operation = new Jupeter_RM09_CAN81Sphr<short>(_process, dynamic_cast<ParamDataTab1DShort &> (_paramData));
  136 + else if (_output == "mag")
  137 + _operation = new Jupeter_RM09_CAN81BMag<short>(_process, dynamic_cast<ParamDataTab1DShort &> (_paramData));
  138 + else
  139 + _operation = new Jupeter_RM09_CAN81Cart<short>(_process, dynamic_cast<ParamDataTab1DShort &> (_paramData));
  140 + }
  141 +
  142 + /**
  143 + * @overload VisitorOfParamData::visit(ParamDataTab1DLogicalData *)
  144 + */
  145 + void visit(ParamDataTab1DLogicalData *) {
  146 + BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("ParamDataTab1DLogicalData data not supported"));
  147 + }
  148 +
  149 + /**
  150 + * @overload VisitorOfParamData::visit(ParamDataTab2DShort *)
  151 + */
  152 + void visit(ParamDataTab2DShort *) {
  153 + BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("ParamDataTab2DShort data not supported"));
  154 + }
  155 +
  156 + /**
  157 + * @overload VisitorOfParamData::visit(ParamDataTab2DFloat *)
  158 + */
  159 + void visit(ParamDataTab2DFloat *) {
  160 + BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("ParamDataTab2DFloat data not supported"));
  161 + }
  162 +
  163 + /**
  164 + * @overload VisitorOfParamData::visit(ParamDataTab2DDouble *)
  165 + */
  166 + void visit(ParamDataTab2DDouble *) {
  167 + BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("ParamDataTab2DDouble data not supported"));
  168 + }
  169 +
  170 + /**
  171 + * @overload VisitorOfParamData::visit(ParamDataTab2DLongDouble *)
  172 + */
  173 + void visit(ParamDataTab2DLongDouble *) {
  174 + BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("ParamDataTab2DLongDouble data not supported"));
  175 + }
  176 +
  177 + /**
  178 + * @overload VisitorOfParamData::visit(ParamDataTab2DInt *)
  179 + */
  180 + void visit(ParamDataTab2DInt *) {
  181 + BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("ParamDataTab2DInt data not supported"));
  182 + }
  183 +
  184 + /**
  185 + * @overload VisitorOfParamData::visit(ParamDataTab2DLogicalData *)
  186 + */
  187 + void visit(ParamDataTab2DLogicalData *) {
  188 + BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_ERROR_UNKNOWN) << AMDA::ex_msg("ParamDataTab2DLogicalData data not supported"));
  189 + }
  190 +
  191 + /**
  192 + * @brief get the TimeShifted parameterized operation.
  193 + */
  194 + Operation * getOperation() const {
  195 + return _operation;
  196 + }
  197 +
  198 +private:
  199 + Process &_process;
  200 + ParamData &_paramData;
  201 + Operation *_operation;
  202 + std::string _output;
  203 + bool _can;
  204 +};
  205 +
  206 +
  207 +}
  208 +#endif /* JUPETER_JRM09_CAN81CREATOR_HH */
  209 +
... ...
src/ExternLib/Jupeter_JRM09_CAN_81/Jupeter_JRM09_CAN81_MagProcess.cc 0 → 100644
  Diff comments   View file @38bebd9
... ... @@ -0,0 +1,51 @@
  1 +/*
  2 + * To change this license header, choose License Headers in Project Properties.
  3 + * To change this template file, choose Tools | Templates
  4 + * and open the template in the editor.
  5 + */
  6 +
  7 +/*
  8 + * File: Jupeter_JRM09_CAN81_MagProcess.cc
  9 + * Author: hacene
  10 + *
  11 + * Created on June 5, 2020, 10:30 AM
  12 + */
  13 +#include <stdlib.h>
  14 +#include <string>
  15 +
  16 +#include "Operation.hh"
  17 +#include "ParameterManager.hh"
  18 +#include "MorschhauserBmagProcess.hh"
  19 +#include "MorschhauserCreator.hh"
  20 +#include "ParameterCreatorFromExpression.hh"
  21 +#include "Jupeter_JRM09_CAN81_MagProcess.hh"
  22 +
  23 +using namespace std;
  24 +using namespace boost;
  25 +using namespace log4cxx;
  26 +namespace AMDA {
  27 +namespace Parameters {
  28 +
  29 +Jupeter_JRM09_CAN81_MagProcess::Jupeter_JRM09_CAN81_MagProcess(Parameter &parameter, bool can) : SingleParamProcess_CRTP(parameter), _can(can){
  30 +}
  31 +
  32 +Jupeter_JRM09_CAN81_MagProcess::Jupeter_JRM09_CAN81_MagProcess(const Jupeter_JRM09_CAN81_MagProcess& orig) {
  33 +}
  34 +
  35 +Jupeter_JRM09_CAN81_MagProcess::~Jupeter_JRM09_CAN81_MagProcess() {
  36 +}
  37 +
  38 + TimeStamp Jupeter_JRM09_CAN81_MagProcess::init() {
  39 + TimeStamp timeStamp = _parameterInput->init( this, _timeIntervalList);
  40 + Parameter::InfoList lInfoList = _parameterInput->getInfoList();
  41 + _parameter.getInfoList().insert(lInfoList.begin(), lInfoList.end());
  42 + _paramInput = _parameterInput->getParamData(this).get();
  43 + Jupeter_RM09_CAN81Creator lJupeter_RM09_CAN81Creator(*this, *_paramInput, "sphr", _can);
  44 + _operation = lMorschhauserCreator.getOperation();
  45 + _paramData = ParamDataSPtr(_operation->getParamOutput());
  46 + _paramData->setMinSampling(_paramInput->getMinSampling());
  47 +
  48 + return timeStamp;
  49 + }
  50 +}
  51 +}
... ...
src/ExternLib/Jupeter_JRM09_CAN_81/Jupeter_JRM09_CAN81_MagProcess.hh 0 → 100644
  Diff comments   View file @38bebd9
... ... @@ -0,0 +1,40 @@
  1 +/*
  2 + * To change this license header, choose License Headers in Project Properties.
  3 + * To change this template file, choose Tools | Templates
  4 + * and open the template in the editor.
  5 + */
  6 +
  7 +/*
  8 + * File: Jupeter_JRM09_CAN81_MagProcess.hh
  9 + * Author: AKKA IS
  10 + *
  11 + * Created on June 5, 2020, 10:30 AM
  12 + */
  13 +
  14 +#ifndef JUPETER_JRM09_CAN81_MAGPROCESS_HH
  15 +#define JUPETER_JRM09_CAN81_MAGPROCESS_HH
  16 +#include "SingleParamProcess.hh"
  17 +
  18 +
  19 +
  20 +namespace AMDA {
  21 +namespace Parameters {
  22 +class Jupeter_JRM09_CAN81_MagProcess : public AMDA::Parameters::SingleParamProcess_CRTP<Jupeter_JRM09_CAN81_MagProcess> {
  23 +public:
  24 + Jupeter_JRM09_CAN81_MagProcess(Parameter &parameter, bool can);
  25 + Jupeter_JRM09_CAN81_MagProcess(const Jupeter_JRM09_CAN81_MagProcess& orig, Parameter &pParameter, bool can);
  26 + virtual ~Jupeter_JRM09_CAN81_MagProcess();
  27 + /**
  28 + * @overload DataWriter::init()
  29 + */
  30 + TimeStamp init();
  31 +private:
  32 + bool _can;
  33 +
  34 +};
  35 +
  36 +}
  37 +}
  38 +
  39 +#endif /* JUPETER_JRM09_CAN81_MAGPROCESS_HH */
  40 +
... ...
src/ExternLib/Jupeter_JRM09_CAN_81/Jupeter_JRM09_CAN81_constants.hh 0 → 100644
  Diff comments   View file @38bebd9
... ... @@ -0,0 +1,63 @@
  1 +/*
  2 + * To change this license header, choose License Headers in Project Properties.
  3 + * To change this template file, choose Tools | Templates
  4 + * and open the template in the editor.
  5 + */
  6 +
  7 +/*
  8 + * File: GetJupeterMag_constants.hh
  9 + * Author: AKKA IS
  10 + * Created on May 29, 2020, 1:47 PM
  11 + */
  12 +
  13 +#ifndef GETJUPETERMAG_CONSTANTS_HH
  14 +#define GETJUPETERMAG_CONSTANTS_HH
  15 +
  16 +#include "Jupeter_JRM09_CAN81_tools.hh"
  17 +
  18 +namespace GETJUPETERMAG_CONSTANTS {
  19 +
  20 +/**
  21 + * The coefficients are from Connerney et al. (2018).
  22 + */
  23 +std::vector<std::vector <const double>> GSMDT = GetJupeterMag_tools::transpose({
  24 + {0.0000000000000000, 410244.7000000000, 11670.40000000000, 4018.600000000000, -34645.40000000000, -18023.60000000000, -20819.60000000000, 598.4000000000000, 10059.20000000000, 9671.799999999999, -2299.500000000000},
  25 + {0.0000000000000000, -71498.30000000000, -56835.80000000000, -37791.10000000000, -8247.600000000000, 4683.900000000000, 9992.900000000000, 4665.900000000000, 1934.400000000000, -3046.200000000000, 2009.700000000000},
  26 + {0.0000000000000000, 0.0000000000000000, 48689.50000000000, 15926.30000000000, -2406.100000000000, 16160.00000000000, 11791.80000000000, -6495.700000000000, -6702.900000000000, 260.9000000000000, 2127.800000000000},
  27 + {0.0000000000000000, 0.0000000000000000, 0.0000000000000000, -2710.500000000000, -11083.80000000000, -16402.00000000000, -12574.70000000000, -2516.500000000000, 153.7000000000000, 2071.300000000000, 3498.300000000000},
  28 + {0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, -17837.20000000000, -2600.700000000000, 2669.700000000000, -6448.500000000000, -4124.200000000000, 3329.600000000000, 2967.600000000000},
  29 + {0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, -3660.700000000000, 1113.200000000000, 1855.300000000000, -867.2000000000001, -2523.100000000000, 16.30000000000000},
  30 + {0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 7584.900000000000, -2892.900000000000, -3740.600000000000, 1787.100000000000, 1806.500000000000},
  31 + {0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 2968.000000000000, -732.4000000000000, -1148.200000000000, -46.50000000000000},
  32 + {0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, -2433.200000000000, 1276.500000000000, 2897.800000000000},
  33 + {0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, -1976.800000000000, 574.5000000000000},
  34 + {0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 1298.900000000000},
  35 +});
  36 +
  37 +std::vector<std::vector <const double>> HSMDT = GetJupeterMag_tools::transpose({
  38 + {0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000},
  39 + {0.0000000000000000, 21330.50000000000, -42027.30000000000, -32957.30000000000, 31994.50000000000, 45347.90000000000, 14533.10000000000, -7626.300000000000, -2409.700000000000, -8467.400000000000, -4692.600000000000},
  40 + {0.0000000000000000, 0.0000000000000000, 19353.20000000000, 42084.50000000000, 27811.20000000000, -749.0000000000000, -10592.90000000000, -10948.40000000000, -11614.60000000000, -1383.800000000000, 4445.800000000000},
  41 + {0.0000000000000000, 0.0000000000000000, 0.0000000000000000, -27544.20000000000, -926.1000000000000, 6268.500000000000, 568.6000000000000, 2633.300000000000, 9287.000000000000, 5697.700000000000, -2378.600000000000},
  42 + {0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 367.1000000000000, 10859.60000000000, 12871.70000000000, 5394.200000000000, -911.9000000000000, -2056.300000000000, -2204.300000000000},
  43 + {0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 9608.400000000000, -4147.800000000000, -6050.800000000000, 2754.500000000000, 3081.500000000000, 164.1000000000000},
  44 + {0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 3604.400000000000, -1526.000000000000, -2446.100000000000, -721.2000000000001, -1361.600000000000},
  45 + {0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, -5684.200000000000, 1207.300000000000, 1352.500000000000, -2031.500000000000},
  46 + {0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, -2887.300000000000, -210.1000000000000, 1411.800000000000},
  47 + {0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 1567.600000000000, -714.3000000000000},
  48 + {0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 1676.500000000000},
  49 +});
  50 +std::vector<std::vector <const double>> RZY = GetJupeterMag_tools::transpose({
  51 + {-0.91419964974697621, 0.36936050381396512, -0.16676875948905945},
  52 + {-0.37460649013519287, -0.92718392610549927, 0.00000000000000000},
  53 + {-0.15462531317480988, 0.062472659656396701, 0.98599600791931152}
  54 +});
  55 +std::vector<std::vector <const double>> RYZ = GetJupeterMag_tools::transpose({
  56 + {-0.91419964974697621, -0.37460649013519287, -0.15462531317480988},
  57 + {0.36936050381396512, -0.92718392610549927, 0.062472659656396701},
  58 + {-0.16676875948905945, 0.00000000000000000, 0.98599600791931152}
  59 +});
  60 +}
  61 +
  62 +#endif /* GETJUPETERMAG_CONSTANTS_HH */
  63 +
... ...
src/ExternLib/Jupeter_JRM09_CAN_81/Jupeter_JRM09_CAN81_tools.hh 0 → 100644
  Diff comments   View file @38bebd9
... ... @@ -0,0 +1,448 @@
  1 +/*
  2 + * To change this license header, choose License Headers in Project Properties.
  3 + * To change this template file, choose Tools | Templates
  4 + * and open the template in the editor.
  5 + */
  6 +
  7 +/*
  8 + * File: tools.hh
  9 + * Author: AKKA IS
  10 + * Created on June 1, 2020, 10:42 AM
  11 + */
  12 +
  13 +#ifndef TOOLS_HH
  14 +#define TOOLS_HH
  15 +
  16 +#include <vector>
  17 +#include <cmath>
  18 +#include <math.h>
  19 +namespace AMDA{
  20 +namespace GETJUPETERMAG_TOOLS {
  21 + using namespace GETJUPETERMAG_CONSTANTS;
  22 +
  23 + class Jupeter_JRM09_CAN81_tools {
  24 + public:
  25 +
  26 + template<typename T1, typename T2>
  27 + std::vector<T1> operator ^(std::vector<T1> a, std::vector<T2> b) {
  28 + if (a.size() != b.size() || b.size() != 3) {
  29 + BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_OPER_NOT_ALLOWED) << AMDA::ex_msg("Vector product on two vector of different size"));
  30 + }
  31 + std::vector<T1> r;
  32 + r[0] = a[1] * b[2] - a[2] * b[1];
  33 + r[1] = a[2] * b[0] - a[0] * b[2];
  34 + r[2] = a[0] * b[1] - a[1] * b[0];
  35 + return r;
  36 + }
  37 +
  38 + template<typename T1>
  39 + std::vector<std::vector<T1>> transpose(std::vector<std::vector<T1>> a) {
  40 + // a[rows][cols]
  41 + int rows = a[0].size();
  42 + int cols = a.size();
  43 + std::vector<std::vector < T1>> a_trans;
  44 + for (int i = 0; i < rows; i++)
  45 + for (int j = 0; j < cols; j++)
  46 + a_trans [j][i] = a[i][j];
  47 + return a_trans;
  48 + }
  49 +
  50 + template<typename T1>
  51 + std::vector<std::vector<T1>> transpose(T1 a[][]) {
  52 + // a[rows][cols]
  53 + int rows = sizeof a/ sizeof a[0]; //
  54 + int cols = sizeof a[0] / sizeof(int);
  55 + std::vector<std::vector < T1>> a_trans;
  56 + for (int i = 0; i < rows; i++)
  57 + for (int j = 0; j < cols; j++)
  58 + a_trans [j][i] = a[i][j];
  59 + return a_trans;
  60 + }
  61 +
  62 + /**
  63 + Transform from Jupiter Magnetic Coordinate (x, y, z) to System III coordinate (x, y, z)
  64 +More details about the Magnetic coordinate:
  65 +Jupiter Solar Magnetospheric (JSM) coordinates Coordinates used during the Galileo mission,
  66 + which used a dipole tilt of 9.6 degrees and lamda_III = 202 degrees based on the O4 model
  67 +of Connerney (1981). Therefore this coordinates are special for CAN current model.
  68 +More details about the System III coordinate:
  69 +Here refers to right-hand System III (S3RH).
  70 + * */
  71 + template<typename T>
  72 + std::vector<T> Coord_JSM2SIII(std::vector<T> jsm) {
  73 + return RZY ^jsm;
  74 + }
  75 +
  76 + /**
  77 + (INVERSE TRANSFORMATION)
  78 + Transform from Jupiter Magnetic Coordinate (x, y, z) to System III coordinate (x, y, z)
  79 + More details about the Magnetic coordinate:
  80 + Jupiter Solar Magnetospheric (JSM) coordinates Coordinates used during the Galileo mission,
  81 + which used a dipole tilt of 9.6 degrees and lamda_III = 202 degrees based on the O4 model
  82 + of Connerney (1981). Therefore this coordinates are special for CAN current model.
  83 + More details about the System III coordinate:
  84 + Here refers to right-hand System III (S3RH).
  85 + Written by Yuxian
  86 + * */
  87 + template<typename T>
  88 + std::vector<T> Coord_SIII2JSM(std::vector<T> siii) {
  89 + return RYZ^siii;
  90 + }
  91 +
  92 + /**
  93 + Cartesian to Spherical coordinates for positions.
  94 +Note: there are differences between points transformation and that of vectors.
  95 +@Input: x, y, z in cartesian coords.
  96 +@Output: r, theta, phi in degrees
  97 + */
  98 + template<typename T>
  99 + std::vector<T> Car2SphP(std::vector<T> cart) {
  100 + T x = cart[0];
  101 + T y = cart[1];
  102 + T z = cart[2];
  103 + T r = std::sqrt(x * x + y * y + z * z);
  104 + if (r == 0)
  105 + r = 1e-31;
  106 + if (x == 0)
  107 + x = 1e-31;
  108 + T theta = std::acos(z / r)*180.0 / PI;
  109 + T phi = std::acos(y / x)*180.0 / PI;
  110 + if (x < 0)
  111 + phi += 180.0;
  112 + if (phi < 0)
  113 + phi += 360.0;
  114 + std::vector<T> result{r, theta, phi};
  115 + return result;
  116 + }
  117 +
  118 + /**
  119 + Spherical to Cartesian coordinates for positions
  120 + Note: there are differences between points transformation and that of vectors
  121 + @Output: x, y, z in cartesian coords. [xyz, npoints]
  122 + @Input: r, theta, phi in degrees
  123 + */
  124 + template<typename T>
  125 + std::vector<T> Car2SphP(std::vector<T> sph) {
  126 + T r = sph[0];
  127 + T theta = sph[1];
  128 + T phi = sph[2];
  129 +
  130 + return std::vector<T> vect{r * std::sin(theta) * std::cos(phi), r * std::sin(theta) * std::sin(phi), r * std::cos(theta)};
  131 + }
  132 +
  133 + /**
  134 + Cartesian to Spherical coordinates for vectors
  135 + Note: there are differences between points transformation and that of vectors.
  136 + @input: x, y, z in cartesian coords., Bx, By, Bz
  137 + @output: Br, Bt, Bp
  138 + */
  139 + template<typename T1, typename T2>
  140 + std::vector<T1> * Car2SphV(std::vector<T1> posi, std::vector<T1> b) {
  141 + T1 x = posi[0];
  142 + T1 y = posi[1];
  143 + T1 z = posi[2];
  144 + T1 Bx = b[0];
  145 + T1 By = b[1];
  146 + T1 Bz = b[2];
  147 +
  148 + T1 tmp = Car2SphP(posi);
  149 + T1 theta = tmp[1] * PI / 180.0;
  150 + T1 phi = tmp[2] * PI / 180.0;
  151 +
  152 + T1 sint = std::sin(theta);
  153 + T1 sinp = std::sin(phi);
  154 + T1 cost = std::cos(theta);
  155 + T1 cosp = std::cos(phi);
  156 +
  157 + T1 Br = Bx * sint * cosp + By * sint * sinp + Bz*cost;
  158 + T1 Bt = Bx * cost * cosp + By * cost * sinp - Bz*sint;
  159 + T Bp = -Bx * sinp + By*cosp;
  160 +
  161 + return std::vector<T1> vect{Br, Bt, Bp};
  162 + }
  163 +
  164 + /**
  165 + Spherical to Cartesian coordinates for vectors
  166 + Note: there are differences between points transformation and that of vectors.
  167 + Function:
  168 +@Output: x, y, z in cartesian coords.
  169 +@Input: posi[r, theta, phi] in degrees, Br, Bt, Bp
  170 + * */
  171 + template<typename T1, typename T2>
  172 + std::vector<T1> * Sph2CarV(std::vector<T1> posi, std::vector<T1> b) {
  173 + T1 theta = posi[1] * PI / 180.0;
  174 + T1 phi = posi[2] * PI / 180.0;
  175 +
  176 + T1 sint = std::sin(theta);
  177 + T1 sinp = std::sin(phi);
  178 + T1 cost = std::cos(theta);
  179 + T1 cosp = std::cos(phi);
  180 +
  181 + T1 Bx = Br * sint * cosp + Bt * cost * cosp - Bp*sinp;
  182 + T1 By = Br * sint * sinp + Bt * cost * sinp + Bp*cosp;
  183 + T1 Bz = Br * cost - Bt * sint;
  184 +
  185 + return std::vector<T1> vect{Bx, By, Bz};
  186 + }
  187 +
  188 + /**
  189 +;Do a "safe" division of numerator by denominator, without a divide by zero.
  190 + ;Both arguments are assumed to have same number of elements.
  191 + ;Where both the numerator and denomiator are zero, the result is normally returned as 1.,
  192 + ; unless the keyword VALUEBOTH0 is set to a different number, such as zero.
  193 + ;Otherwise, the result is returned as near infinity (1.e36), with same sign
  194 + ;as numerator. Since some calculations, such as trigometric functions like
  195 + ;COS(pi/2) result in very small numbers ~4.37e-8 instead of exactly zero, then any
  196 + ;number in the numerator that is smaller than 1.e-8 will be treated like zero,
  197 + ;This value is slightly larger than cos(!PI/2)
  198 + * */
  199 + template<typename T1, typename T2>
  200 + std::vector<T1> * Safe0Division(std::vector<T1> numerator, std::vector<T2> den, std::vector<int> zeropts = {}, int numzeros = 0, T1 valueboth0 = 1.0) {
  201 + std::vector<T1> res;
  202 + if (zeropts.zise() == 0) {
  203 + for (int i = 0; i < den.size(); i++) {
  204 + if (den == 0)
  205 + zeropts.push_back(i);
  206 + }
  207 + }
  208 + if (numzeros == 0)
  209 + numzeros = zeropts.size();
  210 + if (numzeros == 0) {
  211 + for (int i = 0; i < den.size(); i++) {
  212 + res = numerator / den;
  213 + }
  214 + return res;
  215 + } else {
  216 + std::vector<T2> temp = den;
  217 + for (int i : numzeros)
  218 + temp = 1.0;
  219 + for (int i = 0; i < den.size(); i++) {
  220 + res = numerator / temp;
  221 + }
  222 + // result[zeropts]=1. ;one version used 0 before 1/13/2011
  223 + for (int i : numzeros)
  224 + res = valueboth0;
  225 + // value of result for where both numerator and denominator are zero.
  226 + T1 nearzero = 1.e-8;
  227 + T1 inf = 1.e36;
  228 + for (int i : numzeros) {
  229 + if (numerator < -nearzero) {
  230 + res = -inf;
  231 + } else if (numerator > nearzero) {
  232 + res = inf;
  233 + }
  234 + }
  235 + return res;
  236 + }
  237 + }
  238 +
  239 + template<typename T1, typename T2, typename T3>
  240 + T1 Safe0Division(T1 numerator, T2 den, T3 valueboth0 = 1.0) {
  241 + T1 nearzero = 1.e-8;
  242 + T1 inf = 1.e36;
  243 + T1 res;
  244 + if (den != 0)
  245 + return (T1) numerator / den;
  246 + if (numerator < -nearzero) {
  247 + res = -inf;
  248 + } else if (numerator > nearzero) {
  249 + res = inf;
  250 + } else {
  251 + res = valueboth0;
  252 + }
  253 + return res;
  254 + }
  255 +
  256 + /**
  257 + PRO AllLegendre,x,lmax,mmax,Plm,dPlm,SineTheta=SineTheta
  258 + Routine for calculating Associated Legendre Polynomials.
  259 + Returns results in an array for ALL L and M values, and also provides first derivatives.
  260 + Uses recurence formulas. A normalization function is applied.
  261 + At the time it was developed, the recurence formula that was choosen was to avoid copyright conflicts with
  262 + the Numerical Receipes Book, particularly in the FORTRAN vesion that was distributed with the electric potential models.
  263 + compute Associate Legendre Function P_l^m(x) in Double Precision
  264 + returns an array of dimension (n,lmax+1,mmax+1) , where n is size of x
  265 + if X is out of range ( abs(x)>1 ) then value is returns as if x=1.
  266 + Plm is the returned array of Associated Legendre functions.
  267 + dPlm contains the first derivatives. Omit this parameters if the derivatives are not wanted.
  268 + If doing the derivatives and X contains cos(theta),
  269 + then you should set the SINETHETA keyword parameter to non-zero,
  270 + so that the result contains the product of -sin(theta) and dP(cos(theta))/dtheta.
  271 + If any of the angles theta are negative (so that sine(theta) would be negative),
  272 + then you should provide an array with the Sine of the angle theta in the SINETHETA keyword parameter, or just an array of the signs.
  273 + The default behavior is to apply a Schmidt quasi-normalized factor.
  274 + Set the keyword NORMALIZE to 0 to skip the normalization.
  275 + */
  276 + template<typename T1, typename T2>
  277 + int AllLegendre(T1 x, int lma, int mmax, std::vector<std::vector<T2>> Plm,
  278 + std::vector<std::vector<T2>> dPlm, bool isSineTheta, int SineTheta, bool DoDerivative , int normalize = 1) {
  279 + if (lmax < 1 or mmax < 0 || mmax > lmax) {
  280 + BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_OPER_NOT_ALLOWED) << AMDA::ex_msg("Bad arguments to Legendre calculation"));
  281 + }
  282 + // do a range check, as xx must be in the range of -1 to +1
  283 + for (auto val : x) {
  284 + if (val < -1.0)
  285 + val = -1.0;
  286 + if (val > 1.0)
  287 + val = 1.0;
  288 + }
  289 + // First calculate the P_lm polynomials.
  290 + long double P[lmax + 1][mmax + 1][3];
  291 + //P[*,*,*,0] terms multiply by 1 (0 power)
  292 + //P[*,*,*,1] terms contain multiples of sqrt(1-x^2) (or 1/2 power)
  293 + //P[*,*,*,2] terms contain multiples of (1-x^2), later used in the derivatives
  294 + P[0][0][0] = 1.0;
  295 + P[1][0][0] = x;
  296 + //uses this formula for m=0:
  297 + //{P_l}(x) = \frac{{(2l - 1)x{P_{l - 1}}(x) - (l - 1){P_{l - 2}}(x)}}{l}
  298 + if (lmax > 1) {
  299 + for (int j = 2; j < lmax; j++) {
  300 + P[j][0][0] = ((2. * j - 1) * x * P[j - 1][0][0] - (j - 1) * P[j - 2][0][0]) / j;
  301 + }
  302 + }
  303 + long double fact = 1.0 - std::pow(x, 2);
  304 + long double sfact = std::sqrt(fact);
  305 + /**
  306 +
  307 + Following uses this recursion formula:
  308 + {P_l}^m(x) = P_{l-2}^m(x) - (2l-1)\sqrt{1-{x^2}} P_{l-1}^{m-1}(x)
  309 + Most likely this was originally (in 1997) from Abramowitz et al., or Gradshteyn & Ryzhik
  310 + This formula follows from #6 in Wikipedia at:
  311 + https://en.wikipedia.org/wiki/Associated_Legendre_polynomials#Recurrence_formula
  312 + \sqrt{1 - {x^2}} P_l^m(x) = \frac{1}{{2l + 1}}[ - P_{l + 1}^{m + 1} + P_{l - 1}^{m + 1}(x
  313 + */
  314 + if (mmax > 0) {
  315 + for (int m = 1; m < mmax; m++) {
  316 + for (int l = m; l < lmax; l++) {
  317 + int l2 = (l - 2 > 0) ? l - 2 : 0;
  318 + P[l][m][0] = P[l2][m][0];
  319 + P[l][m][1] = P[l2][m][1] - (2 * l - 1) * P[l - 1][m - 1][0] - (2 * l - 1) * P[l - 1][m - 1][2] * fact;
  320 + P[l][m][2] = P[l2][m][2] - (2 * l - 1) * P[l - 1][m - 1][1]; // sqrt(1-x^2) multiplied by itself becomes (1-X^2)
  321 +
  322 + //The P[*,l,m,0] part holds numbers without (1-x^2)
  323 + //The P[*,l,m,1] part holds multiples of sqrt(1-x^2), kept separate for the use in finding 1st derivatives,
  324 + //The P[*,l,m,2] terms contain multiples of (1-x^2)
  325 +
  326 + }
  327 + }
  328 + }
  329 + if (Doderivative) {
  330 + //Keep track of terms containing powers of (1.-x^2)^(1/2) while the recursion formulas are applied to the derivatinves,
  331 + //then evaluate the result at the end.
  332 + //This process is needed because some exponents will cancel out, with an ending power of 0, producing 1 in the end.
  333 + //If recursion formulas were computed at each of the intermediated steps,
  334 + //then a divide by zero error may occur in (1.-x^2)^(-1/2) terms where |x| is 1.
  335 + long double dP[n][lmax + 1][mmax + 1][3];
  336 + /**
  337 + *
  338 + * dP[*,*,*,0] terms have no (1-x^2) multiples
  339 + *dP[*,*,*,1] terms contain multiples of 1./sqrt(1-x^2) (or -1/2 power)
  340 + */
  341 + for (int i = 0; i < n; i++) {
  342 + dP[1][0][0] = 1.0;
  343 + }
  344 + // formula for m=0 is {dP_l}^m(x) = ( (2l-1)P_{l-1}^m(x) + (2l-1)xdP_{l-1}^m(x)- (l-1)dP_{l-2}^m(x) )/l
  345 +
  346 + if (lmax > 1) {
  347 + for (int i = 0; i < n; i++) {
  348 + for (int l = 2; l < lmax; l++) {
  349 +
  350 + l2m1 = 2. * l - 1;
  351 + dP[l][0][0] = (l2m1 * P[l - 1][0][0] + l2m1 * x * dP[l - 1][0][0]-(l - 1) - dP[l - 2][0][0]) / l;
  352 + }
  353 + }
  354 + }
  355 + /**
  356 + formula for m > 0 is
  357 + \frac{dP_l^m(x)}{dx}= \frac{ -lxP_l^m(x) + (l+m) P_(l-1)^m(x) } { (1-x^2) }
  358 + This formula follows from #11 in Wikipedia at:
  359 + https://en.wikipedia.org/wiki/Associated_Legendre_polynomials#Recurrence_formula
  360 + * */
  361 + bool all_zero = true;
  362 + if (mmax > 0) {
  363 + for (int m = 1; m < mmax; m++) {
  364 + for (int l = m; l < lmax; l++) {
  365 + int lm1 = l - 1;
  366 + dP[l][m][0] = -l * x * P[l][m][2] + (l + m) * P[l - 1][m][2]; // A (1-x^2) in the numerator cancels the division by (1-x^2)
  367 + dP[l][m][1] = -l * x * P[l][m][1] + (l + m) * P[l - 1][m][1];
  368 + if (dP[l][m][1] != 0)
  369 + all_zero = false;
  370 + /**
  371 + P[*,*,*,0] terms multiply by 1 (0 power)
  372 + P[*,*,*,1] terms contain multiples of sqrt(1-x^2) (or 1/2 power)
  373 + P[*,*,*,2] terms contain multiples of (1-x^2) , which when divided by (1-x^2) becomes 1
  374 + dP[*,*,*,0] terms have no (1-x^2) multiples
  375 + dP[*,*,*,1] terms contain multiples of 1./sqrt(1-x^2) (or -1/2 power)
  376 + */
  377 + }
  378 + }
  379 + }
  380 + long double signs = -1.0;
  381 + for (int l = 0; l < lmax; l++) {
  382 + mlimit = (mmax < l) ? mmax : l;
  383 + // if there is any chance that the angle theta could have negative values,
  384 + // then the keyword SineTheta needs to contain the actual value of sine(theta), or just the signs, in an array of the same size as x
  385 + if (isSineTheta) {
  386 + // Set SineTheta if X is Cosine(Theta), so need to multiply derivative by -sine(theta)
  387 + // Since sin^2(theta) +cos^2(theta)=1, then sin(theta)=sqrt(1-x^2) , which is already in the variabel sfact.
  388 + dPlm[l][0] = signs *sfact*dP[l][0][0];
  389 + if (mlimit > 0) {
  390 + for (m = 1; m < mlimt; m++) {
  391 + dPlm[l][m] = signs * sfact* dP[l][m][0];
  392 + if (!all_zero) //this next part is done only if there are good data in dP[*,L,M,1]
  393 + dPlm[l][m] += signs * dP[l][m][1]; //multiplying by sfact, sine(theta), cancels out the division by sfact!
  394 + }
  395 + }
  396 + } else {
  397 + //derivative is requested, but the keyword sinetheta was not set
  398 + // dP[*,*,*,0] terms have no (1-x^2) multiples
  399 + //dP[*,*,*,1] terms contain multiples of 1./sqrt(1-x^2) (or -1/2 power)
  400 + dPlm[l][0] = dp[l][0][0];
  401 + if (mlimit > 0) {
  402 + for (int m = 0; m < mlimit; m++) {
  403 + if (!all_zero) {
  404 + if (sfact == 0) {
  405 + dPlm[l][m] += dP[l][m][1] / sfact;
  406 + } else {
  407 + dPlm[l][m] += Safe0Division(dP[l][m][1], sfact, 0.);
  408 + }
  409 + }
  410 +
  411 + }
  412 + }
  413 +
  414 + }
  415 + }
  416 + } // if Doderivative
  417 + for (int l = 0; l < lmax; l++) {
  418 + int mlimit = (mmax < l) ? mmax : l;
  419 + for (int m = 0; m < mlimit; m++) {
  420 + long double nf;
  421 + if (m == 0 || normalize == 0) {
  422 + nf = 1.0;
  423 + } else {
  424 + long double asign = (m % 2 == 0) ? 1 : -1; // If m is odd, then multiply by -1
  425 + // nf=asign*SQRT( 2. *factorial(l-m)/factorial(l+m) ) ;apply the Schmidt quasi-normalized function.
  426 + nf = asign * std::sqrt(2.0 * exp(lgamma(1 - m + 1)) - lgamma(1 + m + 1)); // more accurate way to get (l-m)!/(l+m)!
  427 + }
  428 + Plm[l][m] = nf * (P[l][m][0] + P[l][m][1] * sfact);
  429 + // P[*,*,*,0] terms multiply by 1 (0 power)
  430 + // P[*,*,*,1] terms contain multiples of sqrt(1-x^2) (or 1/2 power), in variable sfact
  431 + // P[*,*,*,2] terms contain multiples of (1-x^2) , in variable fact
  432 + if (DoDerivative)
  433 + dPlm[l][m] *= nf;
  434 + }
  435 + }
  436 +
  437 + return 0;
  438 +
  439 + }
  440 +
  441 + private:
  442 + const double PI = std::atan(1)*4;
  443 +
  444 + };
  445 +} // namespace
  446 +}// namespace amda
  447 +#endif /* TOOLS_HH */
  448 +
... ...
src/Parameters/DataTypeMath.hh
  Diff comments   View file @38bebd9
... ... @@ -322,7 +322,7 @@ Type1 alfvenVelocity(Type1 density, Type2 mag) {
322 322 * @brief FAIRFIELD 1970 model of neutral sheet position
323 323 * @details
324 324 * z_sm-std::sqrt(pow(11.0, 2.0) - pow(11.0, 2.0) / pow(15.0, 2.0) * pow(y_sm, 2.0)) * std::sin(x_ss);
325   - * where y_sm and z_sm are the solar magnetospheric coordinates of the spacecraft and x_ss is geomagnetic latitude of the sun
  325 + * where y_sm and z_sm are the solar magnetospheric coordinates of the spacecraft and x_ss is geomagnetic latitude of the sun (tilt angle)
326 326 */
327 327 template <typename Type1, typename Type2, typename Type3>
328 328 Type1 fairfield70(Type1 y_sm, Type2 z_sm,Type3 x_ss) {
... ...