magneticFieldUtils.hpp

Go to the documentation of this file.

/****************************************************************************** 
* Copyright (c) ATTX INC 2025. All Rights Reserved.
*
* This software and associated documentation (the "Software") are the 
* proprietary and confidential information of ATTX INC. The Software is 
* furnished under a license agreement between ATTX and the user organization 
* and may be used or copied only in accordance with the terms of the agreement.
* Refer to 'license/attx_license.adoc' for standard license terms.
*
* EXPORT CONTROL NOTICE: THIS SOFTWARE MAY INCLUDE CONTENT CONTROLLED UNDER THE
* INTERNATIONAL TRAFFIC IN ARMS REGULATIONS (ITAR) OR THE EXPORT ADMINISTRATION 
* REGULATIONS (EAR99). No part of the Software may be used, reproduced, or 
* transmitted in any form or by any means, for any purpose, without the express 
* written permission of ATTX INC.
******************************************************************************/ 
#ifndef UTILS_MAGNETIC_FIELD_UTILS_H
#define UTILS_MAGNETIC_FIELD_UTILS_H
 
#include <string>
#include <cmath>
#include <fstream>
#include <sstream>
#include <mutex>
#include "utils/fileutils.h"
#include "core/clockwerkerrors.h"
 
namespace warptwin {
    // Constants for spherical harmonics
    const int MAX_HARMONIC_COEFFICIENTS = 12;
    const int NUMBER_OF_READ_COEFFICIENTS = (MAX_HARMONIC_COEFFICIENTS*(MAX_HARMONIC_COEFFICIENTS+1))/2 + MAX_HARMONIC_COEFFICIENTS;

    int readWMMCoefficientsFile(const std::string &filename,
                                double g[NUMBER_OF_READ_COEFFICIENTS],
                                double h[NUMBER_OF_READ_COEFFICIENTS],
                                double gdot[NUMBER_OF_READ_COEFFICIENTS],
                                double hdot[NUMBER_OF_READ_COEFFICIENTS],
                                double &epoch) {
 
        // First verify our file exists, then open
        if(!fileExist(filename)) {return ERROR_NOT_FOUND;}
        std::ifstream infile(filename.c_str());
        if (!infile.is_open()) {return ERROR_NOT_FOUND;}
 
        // Get the epoch from the header line
        std::string line;
        if (!std::getline(infile, line)) return ERROR_INVALID_RANGE; // File empty
        std::istringstream header_stream(line);
        header_stream >> epoch;
 
        // Loop through each line until the end-of-data line is reached
        int index = 0;
        while (index < NUMBER_OF_READ_COEFFICIENTS && std::getline(infile, line)) {
            // Stop if at the end-of-data line
            if (line.find("9999999999") == 0) {break;}
 
            // Read data from line
            int n, m;
            double g_val, h_val, gdot_val, hdot_val;
            std::istringstream iss(line);
            if (!(iss >> n >> m >> g_val >> h_val >> gdot_val >> hdot_val)) {
                // Malformed line, skipping
                continue;
            }
 
            // Save variables
            g[index] = g_val;
            h[index] = h_val;
            gdot[index] = gdot_val;
            hdot[index] = hdot_val;
            ++index;
        }
 
        // Close file to prevent data leaks
        infile.close();
 
        return NO_ERROR;
 
    }

    double muWMM(double input) {
        return std::sin(input);
    }

    double dmuWMM(double input) {
        return std::cos(input);
    }

    int schmidtSemiNormLegendreWMM(int N, double x, double *Pcup, double *dPcup) {
 
        // Check input values to confirm they are in valid range
        if (N < 1 or std::abs(x) > 1.0) {return ERROR_INVALID_RANGE;}
 
        // Static cache for normalization factors
        static std::vector<double> norm;
        static int cached_N = -1;
 
        // Only recompute if N changes
        if (N != cached_N) {
            const int NumTerms = (N + 1) * (N + 2) / 2;
            norm.assign(NumTerms + 1, 1.0);
 
            auto idx = [](int n, int m) { return n * (n + 1) / 2 + m; };
            norm[0] = 1.0;
            for (int n = 1; n <= N; ++n) {
                int base = idx(n, 0);
                int prev = idx(n - 1, 0);
                norm[base] = norm[prev] * (2.0 * n - 1) / n;
                for (int m = 1; m <= n; ++m) {
                    int i = idx(n, m);
                    int ip = idx(n, m - 1);
                    double numer = (n - m + 1) * (m == 1 ? 2.0 : 1.0);
                    double denom = n + m;
                    norm[i] = norm[ip] * std::sqrt(numer / denom);
                }
            }
            cached_N = N;
        }
 
        // Initialize the iterative algorithm
        const double z = std::sqrt(1.0 - x * x);
        Pcup[0] = 1.0;
        dPcup[0] = 0.0;
 
        // Lambda function for getting the index
        auto idx = [](int n, int m) { return n * (n + 1) / 2 + m; };
 
        // Compute Gauss-normalized associated Legendre functions
        for (int n = 1; n <= N; ++n) {
            for (int m = 0; m <= n; ++m) {
                const int i = idx(n, m);
                if (n == m) {
                    const int i1 = idx(n - 1, m - 1);
                    Pcup[i] = z * Pcup[i1];
                    dPcup[i] = z * dPcup[i1] + x * Pcup[i1];
                } else if (n == 1 && m == 0) {
                    const int i1 = idx(0, 0);
                    Pcup[i] = x * Pcup[i1];
                    dPcup[i] = x * dPcup[i1] - z * Pcup[i1];
                } else {
                    const int i1 = idx(n - 2, m);
                    const int i2 = idx(n - 1, m);
                    if (m > n - 2) {
                        Pcup[i] = x * Pcup[i2];
                        dPcup[i] = x * dPcup[i2] - z * Pcup[i2];
                    } else {
                        const double k = ((n - 1.0) * (n - 1.0) - m * m) / ((2.0 * n - 1) * (2.0 * n - 3));
                        Pcup[i] = x * Pcup[i2] - k * Pcup[i1];
                        dPcup[i] = x * dPcup[i2] - z * Pcup[i2] - k * dPcup[i1];
                    }
                }
            }
        }
 
        // Apply normalization and sign change for derivative
        for (int n = 1; n <= N; ++n) {
            for (int m = 0; m <= n; ++m) {
                int i = idx(n, m);
                double scale = norm[i];
                Pcup[i] *= scale;
                dPcup[i] = -dPcup[i] * scale;
            }
        }
 
        return NO_ERROR;
 
    }

    int schmidtSemiNormLegendreSingularityWMM(int N, double x, double *PcupS) {
 
        // Check input values to confirm they are in valid range
        if (N < 1 || std::abs(x) > 1.0) return ERROR_INVALID_RANGE;
 
        // Static cache for normalization factors
        static std::vector<double> normS;
        static int cached_N = -1;
 
        // Only recompute if N changes
        if (N != cached_N) {
            normS.assign(N + 1, 1.0);
            double schmidtQuasiNorm1 = 1.0;
            for (int n = 1; n <= N; ++n) {
                double schmidtQuasiNorm2 = schmidtQuasiNorm1 * (2.0 * n - 1) / n;
                double schmidtQuasiNorm3 = schmidtQuasiNorm2 * std::sqrt((double)(n * 2) / (n + 1));
                normS[n] = schmidtQuasiNorm3;
                schmidtQuasiNorm1 = schmidtQuasiNorm2;
            }
            cached_N = N;
        }
 
 
        // Deal with the singularities at the geographic poles, see section 1.4 in WMM technical report
        PcupS[0] = 1.0;
        for (int n = 1; n <= N; ++n) {
            /* Compute the ratio between the Gauss-normalized associated Legendre
            functions and the Schmidt quasi-normalized version. This is equivalent to
            sqrt((m==0?1:2)*(n-m)!/(n+m!))*(2n-1)!!/(n-m)!  */
            // Compute normalization factors
            if (n == 1) {
                PcupS[n] = PcupS[n - 1];
            } else {
                double k = ((n - 1) * (n - 1) - 1.0) / ((2.0 * n - 1) * (2.0 * n - 3));
                PcupS[n] = (x * PcupS[n - 1] - k * PcupS[n - 2]) * normS[n];
            }
        }
 
        return NO_ERROR;
 
    }
 
}
 
#endif