Merge pull request #1264 from AVSLab/feature/new-wmm-coef

Update WMM coef. file to 2025 version

Author: schaubh

docs/source/Support/bskReleaseNotes.rst

@@ -36,6 +36,7 @@ Version |release|
 - Fixed range of grammatical typos in the Basilisk documentation
 - Ensure the HTML documentation figures from the example scenario look more consistent
 - Avoid a memory leak when saving off figures in the CI Unit tests to build the online documentation
+- Updated the World Magnetic Model (WMM) Coefficients file to the 2025 version.
 
 
 Version 2.9.0 (Jan. 28 2026)

examples/scenarioMagneticFieldWMM.py

@@ -50,7 +50,7 @@
 vector components with respect to the inertial frame.
 
 As this :ref:`MagneticFieldWMM` model is specific to Earth, there are
-no parameters to set of tune.  Rather, the ``WMM.COF`` WMM coefficient
+no parameters to set of tune.  Rather, the ``WMM2025.COF`` WMM coefficient
 file is set via the ``configureWMMFile()`` method and loaded on reset.
 
 The default planet's position vector is assumed to be the inertial
@@ -127,6 +127,7 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
+
 # The path to the location of Basilisk
 # Used to get the location of supporting data.
 from Basilisk import __path__
@@ -138,13 +139,19 @@
 # import simulation related support
 from Basilisk.simulation import spacecraft
 from Basilisk.simulation import magneticFieldWMM
+
 # general support file with common unit test functions
 # import general simulation support files
-from Basilisk.utilities import (SimulationBaseClass, macros, orbitalMotion,
-                                simIncludeGravBody, unitTestSupport)
+from Basilisk.utilities import (
+    SimulationBaseClass,
+    macros,
+    orbitalMotion,
+    simIncludeGravBody,
+    unitTestSupport,
+)
 from Basilisk.utilities.supportDataTools.dataFetcher import get_path, DataFile
 
-#attempt to import vizard
+# attempt to import vizard
 from Basilisk.utilities import vizSupport
 
 
@@ -158,7 +165,6 @@ def run(show_plots, orbitCase):
 
     """
 
-
     # Create simulation variable names
     simTaskName = "simTask"
     simProcessName = "simProcess"
@@ -172,7 +178,7 @@ def run(show_plots, orbitCase):
     dynProcess = scSim.CreateNewProcess(simProcessName)
 
     # create the dynamics task and specify the integration update time
-    simulationTimeStep = macros.sec2nano(60.)
+    simulationTimeStep = macros.sec2nano(60.0)
     dynProcess.addTask(scSim.CreateNewTask(simTaskName, simulationTimeStep))
 
     #
@@ -189,7 +195,7 @@ def run(show_plots, orbitCase):
     # setup Gravity Body
     gravFactory = simIncludeGravBody.gravBodyFactory()
     planet = gravFactory.createEarth()
-    planet.isCentralBody = True          # ensure this is the central gravitational body
+    planet.isCentralBody = True  # ensure this is the central gravitational body
     mu = planet.mu
     req = planet.radEquator
 
@@ -203,15 +209,19 @@ def run(show_plots, orbitCase):
     magModule.ModelTag = "WMM"
 
     # set the minReach and maxReach values if on an elliptic orbit
-    if orbitCase == 'elliptical':
-        magModule.envMinReach = 10000*1000.
-        magModule.envMaxReach = 20000*1000.
+    if orbitCase == "elliptical":
+        magModule.envMinReach = 10000 * 1000.0
+        magModule.envMaxReach = 20000 * 1000.0
 
     # set epoch date/time message
-    epochMsg = unitTestSupport.timeStringToGregorianUTCMsg('2019 June 27, 10:23:0.0 (UTC)')
+    epochMsg = unitTestSupport.timeStringToGregorianUTCMsg(
+        "2019 June 27, 10:23:0.0 (UTC)"
+    )
 
     # add spacecraft to the magnetic field module so it can read the sc position messages
-    magModule.addSpacecraftToModel(scObject.scStateOutMsg)  # this command can be repeated if multiple
+    magModule.addSpacecraftToModel(
+        scObject.scStateOutMsg
+    )  # this command can be repeated if multiple
 
     # add the magnetic field module to the simulation task stack
     scSim.AddModelToTask(simTaskName, magModule)
@@ -221,12 +231,12 @@ def run(show_plots, orbitCase):
     #
     # setup the orbit using classical orbit elements
     oe = orbitalMotion.ClassicElements()
-    rPeriapses = req*1.1     # meters
-    if orbitCase == 'circular':
+    rPeriapses = req * 1.1  # meters
+    if orbitCase == "circular":
         oe.a = rPeriapses
         oe.e = 0.0000
-    elif orbitCase == 'elliptical':
-        rApoapses = req*3.5
+    elif orbitCase == "elliptical":
+        rApoapses = req * 3.5
         oe.a = (rPeriapses + rApoapses) / 2.0
         oe.e = 1.0 - rPeriapses / oe.a
     else:
@@ -249,8 +259,8 @@ def run(show_plots, orbitCase):
 
     # set the simulation time
     n = np.sqrt(mu / oe.a / oe.a / oe.a)
-    P = 2. * np.pi / n
-    simulationTime = macros.sec2nano(1. * P)
+    P = 2.0 * np.pi / n
+    simulationTime = macros.sec2nano(1.0 * P)
 
     # connect messages
     magModule.epochInMsg.subscribeTo(epochMsg)
@@ -259,17 +269,22 @@ def run(show_plots, orbitCase):
     #   Setup data logging before the simulation is initialized
     #
     numDataPoints = 100
-    samplingTime = unitTestSupport.samplingTime(simulationTime, simulationTimeStep, numDataPoints)
+    samplingTime = unitTestSupport.samplingTime(
+        simulationTime, simulationTimeStep, numDataPoints
+    )
     dataLog = scObject.scStateOutMsg.recorder(samplingTime)
     magLog = magModule.envOutMsgs[0].recorder(samplingTime)
     scSim.AddModelToTask(simTaskName, dataLog)
     scSim.AddModelToTask(simTaskName, magLog)
 
     # if this scenario is to interface with the BSK Viz, uncomment the following line
     if vizSupport.vizFound:
-        viz = vizSupport.enableUnityVisualization(scSim, simTaskName, scObject,
-                                                  # saveFile=fileName,
-                                                  )
+        viz = vizSupport.enableUnityVisualization(
+            scSim,
+            simTaskName,
+            scObject,
+            # saveFile=fileName,
+        )
         viz.epochInMsg.subscribeTo(epochMsg)
 
         viz.settings.show24hrClock = 1
@@ -301,41 +316,56 @@ def run(show_plots, orbitCase):
     plt.figure(1)
     fig = plt.gcf()
     ax = fig.gca()
-    ax.ticklabel_format(useOffset=False, style='sci')
-    ax.get_yaxis().set_major_formatter(plt.FuncFormatter(lambda x, loc: "{:,}".format(int(x))))
+    ax.ticklabel_format(useOffset=False, style="sci")
+    ax.get_yaxis().set_major_formatter(
+        plt.FuncFormatter(lambda x, loc: "{:,}".format(int(x)))
+    )
     rData = []
     for idx in range(0, len(posData)):
         rMag = np.linalg.norm(posData[idx])
-        rData.append(rMag / 1000.)
-    plt.plot(timeAxis / P, rData, color='#aa0000')
-    if orbitCase == 'elliptical':
-        plt.plot(timeAxis / P, [magModule.envMinReach/1000.]*len(rData), color='#007700', dashes=[5, 5, 5, 5])
-        plt.plot(timeAxis / P, [magModule.envMaxReach / 1000.] * len(rData),
-                 color='#007700', dashes=[5, 5, 5, 5])
-
-    plt.xlabel('Time [orbits]')
-    plt.ylabel('Radius [km]')
-    plt.ylim(min(rData)*0.9, max(rData)*1.1)
+        rData.append(rMag / 1000.0)
+    plt.plot(timeAxis / P, rData, color="#aa0000")
+    if orbitCase == "elliptical":
+        plt.plot(
+            timeAxis / P,
+            [magModule.envMinReach / 1000.0] * len(rData),
+            color="#007700",
+            dashes=[5, 5, 5, 5],
+        )
+        plt.plot(
+            timeAxis / P,
+            [magModule.envMaxReach / 1000.0] * len(rData),
+            color="#007700",
+            dashes=[5, 5, 5, 5],
+        )
+
+    plt.xlabel("Time [orbits]")
+    plt.ylabel("Radius [km]")
+    plt.ylim(min(rData) * 0.9, max(rData) * 1.1)
     figureList = {}
     pltName = fileName + "1" + orbitCase
     figureList[pltName] = plt.figure(1)
 
     plt.figure(2)
     fig = plt.gcf()
     ax = fig.gca()
-    ax.ticklabel_format(useOffset=False, style='sci')
-    ax.get_yaxis().set_major_formatter(plt.FuncFormatter(lambda x, loc: "{:,}".format(int(x))))
+    ax.ticklabel_format(useOffset=False, style="sci")
+    ax.get_yaxis().set_major_formatter(
+        plt.FuncFormatter(lambda x, loc: "{:,}".format(int(x)))
+    )
     for idx in range(3):
-        plt.plot(timeAxis / P, magData[:, idx] *1e9,
-                 color=unitTestSupport.getLineColor(idx, 3),
-                 label=r'$B\_N_{' + str(idx) + '}$')
-    plt.legend(loc='lower right')
-    plt.xlabel('Time [orbits]')
-    plt.ylabel('Magnetic Field [nT]')
+        plt.plot(
+            timeAxis / P,
+            magData[:, idx] * 1e9,
+            color=unitTestSupport.getLineColor(idx, 3),
+            label=r"$B\_N_{" + str(idx) + "}$",
+        )
+    plt.legend(loc="lower right")
+    plt.xlabel("Time [orbits]")
+    plt.ylabel("Magnetic Field [nT]")
     pltName = fileName + "2" + orbitCase
     figureList[pltName] = plt.figure(2)
 
-
     if show_plots:
         plt.show()
 
@@ -345,13 +375,12 @@ def run(show_plots, orbitCase):
     return figureList
 
 
-
 #
 # This statement below ensures that the unit test scrip can be run as a
 # stand-along python script
 #
 if __name__ == "__main__":
     run(
-        True,          # show_plots
-        'elliptical',  # orbit Case (circular, elliptical)
+        True,  # show_plots
+        "elliptical",  # orbit Case (circular, elliptical)
     )