Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Refactor/bsk 321 denton flux module #322

Merged
merged 4 commits into from
May 30, 2023
Merged

Refactor/bsk 321 denton flux module #322

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
9322f32
Clean up code
JulianHammerl May 5, 2023
b4d0bd8
Change flux unit
JulianHammerl May 5, 2023
2f3df72
Use numpy assert capability for UnitTest
JulianHammerl May 30, 2023
3e20b24
Add release notes
JulianHammerl May 8, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 3 additions & 4 deletions src/architecture/msgPayloadDefC/PlasmaFluxMsgPayload.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,12 +23,11 @@
#define MAX_PLASMA_FLUX_SIZE 50


//!@brief GEO space weather flux message definition.
//!@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
15 changes: 9 additions & 6 deletions src/simulation/environment/dentonFluxModel/_UnitTest/test_dentonFluxModel.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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),
JulianHammerl marked this conversation as resolved.
Show resolved Hide resolved
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)
Expand Down
69 changes: 31 additions & 38 deletions src/simulation/environment/dentonFluxModel/dentonFluxModel.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,25 +20,16 @@
#include "simulation/environment/dentonFluxModel/dentonFluxModel.h"
#include "architecture/utilities/linearAlgebra.h"
#include "architecture/utilities/astroConstants.h"
JulianHammerl marked this conversation as resolved.
Show resolved Hide resolved

#include <iostream>
#include <cstring>
#include <fstream>
#include <cmath> // trig
#include <cmath>

/*! This is the constructor for the module class. It sets default variable
values and initializes the various parts of the model */

// Final Desired Constructor
DentonFluxModel::DentonFluxModel()
{
}
DentonFluxModel::DentonFluxModel() = default;

/*! Module Destructor */
DentonFluxModel::~DentonFluxModel()
{

}
DentonFluxModel::~DentonFluxModel() = default;

/*! This method is used to reset the module and checks that required input messages are connect.
@param CurrentSimNanos current simulation time in nano-seconds
Expand Down Expand Up @@ -82,7 +73,8 @@ void DentonFluxModel::Reset(uint64_t CurrentSimNanos)
// Check that the Kp index is a string of length 2
if (!(this->kpIndex.length() == 2))
{
bskLogger.bskLog(BSK_ERROR, "DentonFluxModel.kpIndex must be a string of length 2, such as '1-', '3o', '4+' etc.");
bskLogger.bskLog(BSK_ERROR,
"DentonFluxModel.kpIndex must be a string of length 2, such as '1-', '3o', '4+' etc.");
}
// Convert Kp index (such as '0o', '1-', '5+' etc.) to Kp index counter (int 0-27)
char kpMain = this->kpIndex[0]; // main Kp index, between 0 and 9
Expand All @@ -106,7 +98,8 @@ void DentonFluxModel::Reset(uint64_t CurrentSimNanos)
// Check that Kp index is between 0o and 9o (corresponding to Kp index counter 0-27)
if (this->kpIndexCounter < 0 || this->kpIndexCounter > MAX_NUM_KPS - 1)
{
bskLogger.bskLog(BSK_ERROR, "DentonFluxModel: Kp index must be between 0o and 9o. Indices 0- and 9+ do not exist.");
bskLogger.bskLog(BSK_ERROR,
"DentonFluxModel: Kp index must be between 0o and 9o. Indices 0- and 9+ do not exist.");
}

// convert energies to log10 values
Expand All @@ -119,7 +112,9 @@ void DentonFluxModel::Reset(uint64_t CurrentSimNanos)
// Define Energy Array
double step = (40000 - 1)/this->numOutputEnergies;

// start at 100eV (fluxes of smaller energies are unreliable due to contamination with secondary electrons and photoelectrons, according to Denton)
// start at 100eV
// (fluxes of smaller energies are unreliable due to contamination with secondary electrons and photoelectrons,
// according to Denton)
this->inputEnergies[0] = 100;
for (int i = 1; i < numOutputEnergies; i++)
{
Expand All @@ -132,7 +127,6 @@ void DentonFluxModel::Reset(uint64_t CurrentSimNanos)

}


/*! This is the main method that gets called every time the module is updated. Provide an appropriate description.
@param CurrentSimNanos current simulation time in nano-seconds
@return void
Expand Down Expand Up @@ -168,7 +162,8 @@ void DentonFluxModel::UpdateState(uint64_t CurrentSimNanos)
double tol = 4000e3; // tolerance how far spacecraft can be away from GEO
if (v2Norm(r_BE_N) < r_GEO - tol || v2Norm(r_BE_N) > r_GEO + tol || abs(r_BE_N[2]) > tol)
{
bskLogger.bskLog(BSK_WARNING, "DentonFluxModel: Spacecraft not in GEO regime. Denton Model not valid outside of GEO.");
bskLogger.bskLog(BSK_WARNING,
"DentonFluxModel: Spacecraft not in GEO regime. Denton Model not valid outside of GEO.");
}

// Find local lime from spacecraft and Earth state messages
Expand Down Expand Up @@ -201,9 +196,6 @@ void DentonFluxModel::UpdateState(uint64_t CurrentSimNanos)
eLowerIndex = k-1;
break;
}
else
{
}
}

// IONS: Find nearest neighbors in energy
Expand All @@ -227,9 +219,6 @@ void DentonFluxModel::UpdateState(uint64_t CurrentSimNanos)
iLowerIndex = k-1;
break;
}
else
{
}
}

int localTimeFloor = floor(this->localTime);
Expand All @@ -246,24 +235,26 @@ 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
finalElec = bilinear(localTimeFloor, localTimeCeil, logEnElec[eLowerIndex], logEnElec[eHigherIndex], logInputEnergy, flux11, flux12, flux13, flux14);

// 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);

// ION: Gather four nearest *MEAN* flux values
flux11 = this->mean_i_flux[this->kpIndexCounter][iLowerIndex][localTimeFloor];
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
finalIon = bilinear(localTimeFloor, localTimeCeil, logEnProt[iLowerIndex], logEnProt[iHigherIndex], logInputEnergy, flux11, flux12, flux13, flux14);

// 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];
}

Expand All @@ -278,7 +269,8 @@ void DentonFluxModel::UpdateState(uint64_t CurrentSimNanos)
*/
void DentonFluxModel::calcLocalTime(double r_SE_N[3], double r_BE_N[3])
{
// r_SE_N and r_BE_N are projected onto the equatorial plane to compute angle, thus only x and y components are used (z component is perpendicular to equator)
// r_SE_N and r_BE_N are projected onto the equatorial plane to compute angle,
// thus only x and y components are used (z component is perpendicular to equator)
double r_BE_N_hat[2]; /* unit vector from Earth to spacecraft */
double r_SE_N_hat[2]; /* unit vector from Earth to Sun */
v2Normalize(r_BE_N, r_BE_N_hat);
Expand All @@ -305,7 +297,8 @@ void DentonFluxModel::calcLocalTime(double r_SE_N[3], double r_BE_N[3])
/*! Bilinear interpolation method
@return void
*/
double DentonFluxModel::bilinear(int x1, int x2, double y1, double y2, double y, double f11, double f12, double f13, double f14)
double DentonFluxModel::bilinear(int x1, int x2, double y1, double y2, double y, double f11, double f12, double f13,
double f14)
{
// Define variables
double R1, R2, bilinear = 0.0;
Expand All @@ -325,13 +318,13 @@ double DentonFluxModel::bilinear(int x1, int x2, double y1, double y2, double y,

}


/*! Read in the Denton data file
@param fileName data file name
@param data data array pointer
@return void
*/
void DentonFluxModel::readDentonDataFile(std::string fileName, double data[MAX_NUM_KPS][MAX_NUM_ENERGIES][MAX_NUM_LOCAL_TIMES])
void DentonFluxModel::readDentonDataFile(std::string fileName,
double data[MAX_NUM_KPS][MAX_NUM_ENERGIES][MAX_NUM_LOCAL_TIMES])
{
double temp = 0.0;

Expand All @@ -341,7 +334,8 @@ void DentonFluxModel::readDentonDataFile(std::string fileName, double data[MAX_N
// Read data from file:
inputFile.open(this->dataPath + fileName);

// Read information into array: Data includes information about mean, standard deviation, median and percentiles (7 types of values in total). Only mean is relevant for this module
// Read information into array: Data includes information about mean, standard deviation,
// median and percentiles (7 types of values in total). Only mean is relevant for this module
if (inputFile.is_open()) {
for (int i = 0; i < MAX_NUM_KPS*MAX_NUM_VALUE_TYPES; i++)
{ for (int j = 0; j < MAX_NUM_ENERGIES; j++)
Expand All @@ -367,4 +361,3 @@ void DentonFluxModel::readDentonDataFile(std::string fileName, double data[MAX_N

return;
}

54 changes: 23 additions & 31 deletions src/simulation/environment/dentonFluxModel/dentonFluxModel.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,6 @@

*/


#ifndef DENTONFLUXMODEL_H
#define DENTONFLUXMODEL_H

Expand All @@ -28,56 +27,56 @@
#include "architecture/utilities/bskLogging.h"
#include "architecture/messaging/messaging.h"


#define MAX_NUM_KPS 28
#define MAX_NUM_ENERGIES 40
#define MAX_NUM_LOCAL_TIMES 24
#define MAX_NUM_VALUE_TYPES 7


/*! @brief This module provides the 10-year averaged GEO elecon and ion flux as discussed in the paper by Denton.
*/
class DentonFluxModel: public SysModel {
public:

// Constructor And Destructor
DentonFluxModel();
~DentonFluxModel();

// Methods
void Reset(uint64_t CurrentSimNanos);
void UpdateState(uint64_t CurrentSimNanos);
void Reset(uint64_t CurrentSimNanos) override;
void UpdateState(uint64_t CurrentSimNanos) override;

/* public variables */
int numOutputEnergies = -1; //!< number of energy bins used in the output message
std::string kpIndex = ""; //!< Kp index
std::string dataPath = ""; //!< -- String with the path to the Denton GEO data
std::string eDataFileName = "model_e_array_all.txt"; //!< file name of the electron data file
std::string iDataFileName = "model_i_array_all.txt"; //!< file name of the ion data file
int numOutputEnergies = -1; //!< number of energy bins used in the output message
std::string kpIndex = ""; //!< Kp index
std::string dataPath = ""; //!< -- String with the path to the Denton GEO data
std::string eDataFileName = "model_e_array_all.txt"; //!< file name of the electron data file
std::string iDataFileName = "model_i_array_all.txt"; //!< file name of the ion data file

ReadFunctor<SCStatesMsgPayload> scStateInMsg; //!< spacecraft state input message
ReadFunctor<SpicePlanetStateMsgPayload> earthStateInMsg; //!< Earth planet state input message
ReadFunctor<SpicePlanetStateMsgPayload> sunStateInMsg; //!< sun state input message
ReadFunctor<SCStatesMsgPayload> scStateInMsg; //!< spacecraft state input message
ReadFunctor<SpicePlanetStateMsgPayload> earthStateInMsg; //!< Earth planet state input message
ReadFunctor<SpicePlanetStateMsgPayload> sunStateInMsg; //!< sun state input message

Message<PlasmaFluxMsgPayload> fluxOutMsg; //!< output message with ion and electron fluxes
Message<PlasmaFluxMsgPayload> fluxOutMsg; //!< output message with ion and electron fluxes

BSKLogger bskLogger; //!< -- BSK Logging
BSKLogger bskLogger; //!< -- BSK Logging

private:
void calcLocalTime(double v1[3], double v2[3]);
double bilinear(int, int, double, double, double, double, double, double, double);
void readDentonDataFile(std::string fileName, double data[MAX_NUM_KPS][MAX_NUM_ENERGIES][MAX_NUM_LOCAL_TIMES]);

int kpIndexCounter; //!< Kp index counter (betweeen 0 and 27)
double localTime; /* spacecraft location time relative to sun heading at GEO */
double logEnElec[MAX_NUM_ENERGIES]; /* log of the electron energies */
double logEnProt[MAX_NUM_ENERGIES]; /* log of the proton energies */
double inputEnergies[MAX_NUM_ENERGIES]; /* input energies considered in this module */
int kpIndexCounter; //!< Kp index counter (betweeen 0 and 27)
double localTime; //!< spacecraft location time relative to sun heading at GEO
double logEnElec[MAX_NUM_ENERGIES]; //!< log of the electron energies
double logEnProt[MAX_NUM_ENERGIES]; //!< log of the proton energies
double inputEnergies[MAX_NUM_ENERGIES]; //!< input energies considered in this module

// Electron Flux:
//!< Electron Flux:
double mean_e_flux[MAX_NUM_KPS][MAX_NUM_ENERGIES][MAX_NUM_LOCAL_TIMES];

// Ion Flux:
//!< Ion Flux:
double mean_i_flux[MAX_NUM_KPS][MAX_NUM_ENERGIES][MAX_NUM_LOCAL_TIMES];

// Fill average centre energies, normalized by satellite
//!< Fill average centre energies, normalized by satellite
double enElec[40] = {1.034126, 1.346516, 1.817463, 2.399564,
3.161048, 4.153217, 5.539430, 7.464148,
9.836741, 12.543499, 16.062061, 20.876962,
Expand All @@ -99,13 +98,6 @@ class DentonFluxModel: public SysModel {
2170.886474, 2829.989013, 3691.509765, 4822.499023,
6300.260742, 8217.569335, 10726.390625, 14001.280273,
18276.244140, 23856.085937, 31140.962890, 40649.562500};


void calcLocalTime(double v1[3], double v2[3]); //!< calculate the local time
double bilinear(int, int, double, double, double, double, double, double, double);
void readDentonDataFile(std::string fileName, double data[MAX_NUM_KPS][MAX_NUM_ENERGIES][MAX_NUM_LOCAL_TIMES]);

};


#endif