Change flux unit

From [cm^-2 s^-1 sr^-2 eV^-1] to [m^-2 s^-1 sr^-2 eV^-1]
AVSLab · May 8, 2023 · e4ea709 · e4ea709
1 parent c1fb1f7
commit e4ea709
Show file tree

Hide file tree

Showing 3 changed files with 18 additions and 16 deletions.
diff --git a/src/architecture/msgPayloadDefC/PlasmaFluxMsgPayload.h b/src/architecture/msgPayloadDefC/PlasmaFluxMsgPayload.h
@@ -25,10 +25,9 @@
 
 //!@brief space weather ambient plasma flux message definition.
 typedef struct {
-    double meanElectronFlux[MAX_PLASMA_FLUX_SIZE];     //!< [cm^-2 s^-1 sr^-2 eV^-1] differential flux
-    double meanIonFlux[MAX_PLASMA_FLUX_SIZE];          //!< [cm^-2 s^-1 sr^-2 eV^-1] differential flux
+    double meanElectronFlux[MAX_PLASMA_FLUX_SIZE];     //!< [m^-2 s^-1 sr^-2 eV^-1] differential flux
+    double meanIonFlux[MAX_PLASMA_FLUX_SIZE];          //!< [m^-2 s^-1 sr^-2 eV^-1] differential flux
     double energies[MAX_PLASMA_FLUX_SIZE];             //!< [eV]
 }PlasmaFluxMsgPayload;
 
-
 #endif
diff --git a/src/simulation/environment/dentonFluxModel/_UnitTest/test_dentonFluxModel.py b/src/simulation/environment/dentonFluxModel/_UnitTest/test_dentonFluxModel.py
@@ -166,26 +166,29 @@ def dentonFluxModelTestFunction(show_plots, param1_Kp, param2_LT, param3_z, para
     trueIonFluxData = np.array([0.0] * messaging.MAX_PLASMA_FLUX_SIZE)
     with open(filepath, 'r') as file:
         rows = np.loadtxt(file, delimiter=",", unpack=False)
+        # true flux data provided by Denton is in Units of [cm^-2 s^-1 sr^-2 eV^-1], but DentonFluxModel converts it to
+        # [m^-2 s^-1 sr^-2 eV^-1]. Need to multiply by 1e4
         trueEnergyData[0:module.numOutputEnergies] = rows[0]
-        trueElectronFluxData[0:module.numOutputEnergies] = 10.**(rows[1])
-        trueIonFluxData[0:module.numOutputEnergies] = 10.**(rows[2])
+        trueElectronFluxData[0:module.numOutputEnergies] = 10.**(rows[1]) * 1e4
+        trueIonFluxData[0:module.numOutputEnergies] = 10.**(rows[2]) * 1e4
 
     # make sure module output data is correct
     ParamsString = ' for Kp-Index=' + param1_Kp + ', LT=' + str(param2_LT)
     testFailCount, testMessages = unitTestSupport.compareDoubleArray(
-        trueEnergyData, energyData, accuracy, ('electron and ion energies' + ParamsString),
+        trueEnergyData, energyData, accuracy * 1e4, ('particle energies' + ParamsString),
         testFailCount, testMessages)
     testFailCount, testMessages = unitTestSupport.compareDoubleArray(
-        trueElectronFluxData, electronFluxData, accuracy, ('electron and ion energies' + ParamsString),
+        trueElectronFluxData, electronFluxData, accuracy * 1e4, ('electron flux' + ParamsString),
         testFailCount, testMessages)
     testFailCount, testMessages = unitTestSupport.compareDoubleArray(
-        trueIonFluxData, ionFluxData, accuracy, ('electron and ion energies' + ParamsString),
+        trueIonFluxData, ionFluxData, accuracy * 1e4, ('ion flux' + ParamsString),
         testFailCount, testMessages)
 
     # print out success or failure message
     if testFailCount == 0:
         print("PASSED: " + module.ModelTag)
-        print("This test uses an accuracy value of " + str(accuracy) + " (true values don't have higher precision)")
+        print("This test uses an accuracy value of " + str(accuracy * 1e4) +
+              " (true values don't have higher precision)")
     else:
         print("FAILED " + module.ModelTag)
         print(testMessages)

diff --git a/src/simulation/environment/dentonFluxModel/dentonFluxModel.cpp b/src/simulation/environment/dentonFluxModel/dentonFluxModel.cpp
@@ -235,8 +235,8 @@ void DentonFluxModel::UpdateState(uint64_t CurrentSimNanos)
         flux12 = this->mean_e_flux[this->kpIndexCounter][eHigherIndex][localTimeFloor];
         flux13 = this->mean_e_flux[this->kpIndexCounter][eLowerIndex][localTimeCeil];
         flux14 = this->mean_e_flux[this->kpIndexCounter][eHigherIndex][localTimeCeil];
-        
-        // ELECTRON: Find flux
+
+        // ELECTRON: Find flux (differential flux in units of [cm^-2 s^-1 sr^-2 eV^-1])
         finalElec = bilinear(localTimeFloor, localTimeCeil, logEnElec[eLowerIndex], logEnElec[eHigherIndex],
                              logInputEnergy, flux11, flux12, flux13, flux14);
         finalElec = pow(10.0, finalElec);
@@ -246,15 +246,15 @@ void DentonFluxModel::UpdateState(uint64_t CurrentSimNanos)
         flux12 = this->mean_i_flux[this->kpIndexCounter][iHigherIndex][localTimeFloor];
         flux13 = this->mean_i_flux[this->kpIndexCounter][iLowerIndex][localTimeCeil];
         flux14 = this->mean_i_flux[this->kpIndexCounter][iHigherIndex][localTimeCeil];
-        
-        // ION: Find flux
+
+        // ION: Find flux (differential flux in units of [cm^-2 s^-1 sr^-2 eV^-1])
         finalIon = bilinear(localTimeFloor, localTimeCeil, logEnProt[iLowerIndex], logEnProt[iHigherIndex],
                             logInputEnergy, flux11, flux12, flux13, flux14);
         finalIon = pow(10.0, finalIon);
 
-        // Store the output message
-        fluxOutMsgBuffer.meanElectronFlux[i] = finalElec;
-        fluxOutMsgBuffer.meanIonFlux[i] = finalIon;
+        // Store the output message (differential flux in units of [m^-2 s^-1 sr^-2 eV^-1])
+        fluxOutMsgBuffer.meanElectronFlux[i] = finalElec * 1e4;
+        fluxOutMsgBuffer.meanIonFlux[i] = finalIon * 1e4;
         fluxOutMsgBuffer.energies[i] = inputEnergies[i];
     }