Geodesy.WMM2025.magneticFieldEnu

Parameters, inputs, outputs, defaults and equations from the source file.

function magneticFieldEnu
    "Evaluate WMM2025 at WGS84 ellipsoid altitude and return ENU tesla"
    input Real latitude_deg(unit = "deg");
    input Real longitude_deg(unit = "deg");
    input Real altitudeEllipsoid_m(unit = "m");
    input Real decimalYear = 2026.0;
    output Real magneticFieldWorldEnu_T[3](each unit = "T");
  protected
    Real latitudeGeodetic_rad;
    Real longitude_rad;
    Real sinLatitudeGeodetic;
    Real cosLatitudeGeodetic;
    Real primeVerticalRadius_m;
    Real cylindricalRadius_m;
    Real z_m;
    Real radius_m;
    Real latitudeSpherical_rad;
    Real sinLatitudeSpherical;
    Real cosLatitudeSpherical;
    Real deltaYear;
    Real g[13, 13];
    Real h[13, 13];
    Real P[13, 13];
    Real dP[13, 13];
    Real sinMultipleLongitude[13];
    Real cosMultipleLongitude[13];
    Real basePower;
    Real radialPower;
    Real f1;
    Real f2;
    Real northSpherical_nT;
    Real east_nT;
    Real downSpherical_nT;
    Real north_nT;
    Real down_nT;
    Real cosDeltaLatitude;
    Real sinDeltaLatitude;
    Integer coefficientRow;
  algorithm
    latitudeGeodetic_rad := latitude_deg * pi / 180.0;
    longitude_rad := longitude_deg * pi / 180.0;
    sinLatitudeGeodetic := sin(latitudeGeodetic_rad);
    cosLatitudeGeodetic := cos(latitudeGeodetic_rad);
    primeVerticalRadius_m := Wgs84A_m / sqrt(1.0
      - Wgs84EccentricitySquared * sinLatitudeGeodetic
        * sinLatitudeGeodetic);
    cylindricalRadius_m := (primeVerticalRadius_m + altitudeEllipsoid_m)
      * cosLatitudeGeodetic;
    z_m := (primeVerticalRadius_m * (1.0 - Wgs84EccentricitySquared)
      + altitudeEllipsoid_m) * sinLatitudeGeodetic;
    radius_m := sqrt(cylindricalRadius_m * cylindricalRadius_m + z_m * z_m);
    latitudeSpherical_rad := asin(z_m / radius_m);
    sinLatitudeSpherical := sin(latitudeSpherical_rad);
    cosLatitudeSpherical := cos(latitudeSpherical_rad);
    deltaYear := decimalYear - 2025.0;

    g := zeros(13, 13);
    h := zeros(13, 13);
    coefficientRow := 0;
    for n in 1:MaximumDegree loop
      for m in 0:n loop
        coefficientRow := coefficientRow + 1;
        g[n + 1, m + 1] := coefficients[coefficientRow, 1]
          + deltaYear * coefficients[coefficientRow, 3];
        h[n + 1, m + 1] := coefficients[coefficientRow, 2]
          + deltaYear * coefficients[coefficientRow, 4];
      end for;
    end for;

    // Schmidt semi-normalized associated Legendre functions and latitude
    // derivatives. This fixed degree-12 recurrence is the Modelica port of
    // the NOAA evaluation used by the reference C implementation.
    P := zeros(13, 13);
    dP := zeros(13, 13);
    P[1, 1] := 1.0;
    P[2, 1] := sinLatitudeSpherical;
    dP[2, 1] := cosLatitudeSpherical;
    for n in 2:MaximumDegree loop
      f1 := (2.0 * n - 1.0) / n;
      f2 := (n - 1.0) / n;
      P[n + 1, 1] := f1 * sinLatitudeSpherical * P[n, 1]
        - f2 * P[n - 1, 1];
      dP[n + 1, 1] := n * (P[n, 1]
        - sinLatitudeSpherical * P[n + 1, 1])
          / cosLatitudeSpherical;
    end for;
    for m in 1:MaximumDegree loop
      // Schmidt quasi-normalization has a distinct m=1 seed. Applying the
      // general diagonal recurrence at m=1 incorrectly scales every
      // non-zonal term by 1/sqrt(2).
      P[m + 1, m + 1] := if m == 1 then cosLatitudeSpherical
        else P[m, m] * cosLatitudeSpherical
          * sqrt((2.0 * m - 1.0) / (2.0 * m));
      dP[m + 1, m + 1] := -m * sinLatitudeSpherical
        * P[m + 1, m + 1] / cosLatitudeSpherical;
      if m < MaximumDegree then
        P[m + 2, m + 1] := sinLatitudeSpherical * sqrt(2.0 * m + 1.0)
          * P[m + 1, m + 1];
        dP[m + 2, m + 1] := (sqrt(2.0 * m + 1.0)
          * P[m + 1, m + 1] - sinLatitudeSpherical * (m + 1.0)
            * P[m + 2, m + 1]) / cosLatitudeSpherical;
      end if;
      for n in m + 2:MaximumDegree loop
        f1 := (2.0 * n - 1.0) / sqrt((n + m) * (n - m));
        f2 := sqrt((n - m - 1.0) * (n + m - 1.0)
          / ((n + m) * (n - m)));
        P[n + 1, m + 1] := sinLatitudeSpherical * f1
          * P[n, m + 1] - f2 * P[n - 1, m + 1];
        dP[n + 1, m + 1] := (sqrt((n + m) * (n - m))
          * P[n, m + 1] - n * sinLatitudeSpherical
            * P[n + 1, m + 1]) / cosLatitudeSpherical;
      end for;
    end for;

    sinMultipleLongitude := zeros(13);
    cosMultipleLongitude := zeros(13);
    cosMultipleLongitude[1] := 1.0;
    for m in 1:MaximumDegree loop
      sinMultipleLongitude[m + 1] := sinMultipleLongitude[m]
          * cos(longitude_rad) + cosMultipleLongitude[m]
          * sin(longitude_rad);
      cosMultipleLongitude[m + 1] := cosMultipleLongitude[m]
          * cos(longitude_rad) - sinMultipleLongitude[m]
          * sin(longitude_rad);
    end for;

    northSpherical_nT := 0.0;
    east_nT := 0.0;
    downSpherical_nT := 0.0;
    basePower := GeomagneticReferenceRadius_m / radius_m;
    radialPower := basePower * basePower;
    for n in 1:MaximumDegree loop
      radialPower := radialPower * basePower;
      for m in 0:n loop
        northSpherical_nT := northSpherical_nT - radialPower
          * (g[n + 1, m + 1] * cosMultipleLongitude[m + 1]
            + h[n + 1, m + 1] * sinMultipleLongitude[m + 1])
          * dP[n + 1, m + 1];
        east_nT := east_nT + radialPower / cosLatitudeSpherical * m
          * (g[n + 1, m + 1] * sinMultipleLongitude[m + 1]
            - h[n + 1, m + 1] * cosMultipleLongitude[m + 1])
          * P[n + 1, m + 1];
        downSpherical_nT := downSpherical_nT
          - (n + 1.0) * radialPower
          * (g[n + 1, m + 1] * cosMultipleLongitude[m + 1]
            + h[n + 1, m + 1] * sinMultipleLongitude[m + 1])
          * P[n + 1, m + 1];
      end for;
    end for;
    cosDeltaLatitude := cos(latitudeSpherical_rad - latitudeGeodetic_rad);
    sinDeltaLatitude := sin(latitudeSpherical_rad - latitudeGeodetic_rad);
    north_nT := northSpherical_nT * cosDeltaLatitude
      - downSpherical_nT * sinDeltaLatitude;
    down_nT := northSpherical_nT * sinDeltaLatitude
      + downSpherical_nT * cosDeltaLatitude;
    magneticFieldWorldEnu_T := {east_nT, north_nT, -down_nT} * 1.0e-9;
  end magneticFieldEnu;