Merge pull request #185 from AetherModel/fism

Merge Fism Branch into Develop - this is the last open branch, so development can be completely synched up at this point.

Author: aaronjridley

.gitignore

@@ -39,3 +39,7 @@ _deps
 #vs code settings
 .vscode
 .DS_Store
+
+# fism files
+srcPython/FISM*.nc
+srcPython/fism*.txt

doc/internals/fism.md

@@ -0,0 +1,42 @@
+# FISM-2
+
+FISM-2 can be used as a EUV model.
+
+Needs the FISM-2 files automatically made by `srcPython/fism.py`.
+
+FISM-2 contains the binned flux already, however Aether must still be provided with an
+EUV csv file containing the cross-sections.
+
+The different fism models available in fism.py contain different numbers of bins, so
+the euv file must be different.
+
+|  Model  |  number of bins  | euv file |
+|  :---  |  :-------------: | -------: |
+| HFG    | 23  | euv_solomon.csv |
+| Solomon | 23  | euv_solomon.csv |
+| NEUVAC  | 37/59   | euv.csv / euv_59.csv |
+| EUVAC  |  37  | euv.csv |
+
+The input format, when using fism data:
+
+    "Euv" : {
+        "doUse" : true,
+        "Model" : "fism",
+        "File" : "UA/inputs/euv_59.csv",
+        "fismFile": "fism2_file_59.txt",
+        "IncludePhotoElectrons" : true,
+        "HeatingEfficiency" : 0.05,
+        "dt" : 60.0
+    },
+
+To generate FISM-2 irradiances between two dates, one may call fism.py like so:
+
+srcPython/fism.py 20110319 20110321 -b neuvac
+
+Note that the optional argument '-b' defaults to the binning scheme of the NEUVAC model,
+which employs 59 bins. To use the 37 bins of the EUVAC model, the argument should be
+'euvac', while to use the 23 bins used by the HFG model, the argument should be 'solomon'.
+
+fism.py should always be run before Aether is run using the FISM-2 model. Even though the
+entire FISM-2 irradiances are stored in the repository, Aether will need a separate .csv
+file containing the temporal subset of FISM-2 irradiances in the desired binning scheme.

include/euv.h

@@ -94,6 +94,11 @@ class Euv {
   /// NEUVAC model intercept:
   std::vector<float> neuvac_int;
 
+  // To avoid having to start from 0 each iteration:
+  int fism_prev_index = 0;
+  // Declare this so it is not passed between function:
+  index_file_output_struct fismData;
+  
   // --------------------------------------------------------------------
   // Functions:
 
@@ -109,13 +114,26 @@ class Euv {
    **/
   bool euvac(Times time, Indices indices);
 
+ /**********************************************************************
+     \brief Compute the EUV spectrum given F107 and F107a
+     \param time The times within the model (dt is needed)
+     \param indices Need the F107 and F107a
+   **/
+ bool solomon_hfg(Times time, Indices indices);
+
   /**********************************************************************
-      \brief Compute the EUV spectrum given F107 and F107a
-      \param time The times within the model (dt is needed)
-      \param indices Need the F107 and F107a
-    **/
-  bool solomon_hfg(Times time, Indices indices);
+     \brief returns the FISM spectrum for a given time
+
+    Unlike the other EUV models ([N]EUVAC, Solomon, etc.), the spectrum
+    is read from a file (stored in fismData). This does the same thing
+    as get_index, however FISM is not stored in Indices since it can
+    have variable # of bins
+
+     \param time The times within the model (dt is needed)
+   **/
 
+ bool get_fism(Times time);
+  
   /**********************************************************************
      \brief Compute the EUV spectrum given F107 and F107a (new version)
      \param time The times within the model (dt is needed)
@@ -161,6 +179,17 @@ class Euv {
    **/
   bool read_file();
 
+  /**********************************************************************
+     \brief Read in the FISM  file
+
+     Read in the CSV file with FISM data. This can be made with
+      srcPython/fism.py. The data are read into a index_file_output_struct,
+      where each row is one time, and each col is a "variable". These should
+      match the number of bins in the provided EUV file.
+   **/
+ index_file_output_struct read_fism(std::string fism_filename);
+
+
   /**********************************************************************
      \brief Interprets the EUV CSV rows and returns the relevant row
 

include/inputs.h

@@ -228,6 +228,12 @@ class Inputs {
      \brief returns settings["
      \param
    **/
+  std::string get_euv_fismfile();
+  
+  /**********************************************************************
+     \brief returns settings["
+     \param 
+   **/
   bool get_euv_douse();
 
   /**********************************************************************

share/run/UA/inputs/euv_59_v2.csv

@@ -49,6 +49,9 @@ bool calc_euv(Planets planet,
       didWork = euv.neuvac(time, indices);
     else if (euvModel == "hfg")
       didWork = euv.solomon_hfg(time, indices);
+    else if (euvModel == "fism"){
+      didWork = euv.get_fism(time);
+    }
 
     if (didWork)
       euv.scale_from_1au(planet, time);

src/calc_euv.cpp

@@ -44,6 +44,11 @@ Euv::Euv() {
         }
       }
 
+      // Read in FISM data - does not need to be "slotted"
+      if (input.get_euv_model() == "fism")
+        fismData = read_fism(input.get_euv_fismfile());
+      // Read in NEUVAC data - also does not need to be "slotted"
+      
       // Slot the EUVAC model coefficients:
       if (input.get_euv_model() == "euvac") {
         IsOk = slot_euv("F74113", "", euvac_f74113);
@@ -167,6 +172,62 @@ bool Euv::read_file() {
   return DidWork;
 }
 
+// -------------------------------------------------------------------------------
+// Read in FISM data. FISM files are created with srcPython/fism.py,
+// and the data are read in to an index_file_output_struct.
+// Inside the struct, we have time & each of the "variables" correspond to a
+// FISM bin. This number of bins should match the number of bins in the EUV file
+// -------------------------------------------------------------------------------
+
+index_file_output_struct Euv::read_fism(std::string fism_filename) {
+
+  std::ifstream fismfstream;
+  fismfstream.open(fism_filename);
+  std::vector<std::vector<std::string>>  fism_file;
+  fism_file = read_csv(fismfstream);
+
+  index_file_output_struct fism_contents;
+
+  // one row per time
+  fism_contents.nTimes = fism_file.size();
+  // first six cols are the YYYY,MM,DD,HH,mm,ss (no ms)
+  // the rest are the binned fism data
+  fism_contents.nVars = fism_file[0].size() - 6;
+
+  // check that the user provided the correct EUV file
+  // The number of bins in euv file should match the number of fism bins ("nVars")
+  if (fism_contents.nVars != nWavelengths) {
+    report.error("Number of FISM wavelengths does not match the EUV file provided!");
+    report.error("Either change EUV file or check your FISM file is correct.");
+    IsOk = false;
+  }
+
+  std::vector<int> itime(7, 0);
+  std::vector<std::vector<float>> values; // holds all values
+  std::vector<float> values_tmp(fism_contents.nVars); // holds values in each row
+
+  for (int iLine = 0; iLine < fism_file.size(); iLine ++) {
+
+    itime[0] = stoi(fism_file[iLine][0]);
+    itime[1] = stoi(fism_file[iLine][1]);
+    itime[2] = stoi(fism_file[iLine][2]);
+    itime[3] = stoi(fism_file[iLine][3]);
+    itime[4] = stoi(fism_file[iLine][4]);
+    itime[5] = stoi(fism_file[iLine][5]);
+    itime[6] = 0; // 0 ms
+    fism_contents.times.push_back(time_int_to_real(itime));
+
+    for (int iVar = 0; iVar < fism_contents.nVars; iVar++)
+      values_tmp[iVar] = stof(fism_file[iLine][iVar + 6]);
+
+    values.push_back(values_tmp);
+  }
+
+  fism_contents.values = values;
+
+  return fism_contents;
+}
+
 // ---------------------------------------------------------------------------
 // Match rows in EUV file to different types of things, such as cross
 // sections and spectra
@@ -371,6 +432,52 @@ bool Euv::euvac(Times time,
   return didWork;
 }
 
+// --------------------------------------------------------------------------
+// From the FISM file, interpolate the nearest 2 data to the current time
+// --------------------------------------------------------------------------
+
+bool Euv::get_fism(Times time) {
+  // This is functionally similar to get_indices, however we do not store FISM in
+  // the Indices class since it has variable number of bins.
+  
+  std::string function = "Euv::get_fism";
+  static int iFunction = -1;
+  report.enter(function, iFunction);
+
+  double time_now = time.get_current();
+  bool didWork = true;
+
+  if (fism_prev_index == 0) {
+    // This is probably the first time we're "running" fism.
+    // Make sure the file covers the entire time range of the run.
+    double end_time = time.get_end();
+
+    if (time_now < fismData.times[0] && end_time > fismData.times[-1]) {
+      report.error("FISM data does not cover the entire time range!");
+      report.error("Please check that your FISM file is correct.");
+      didWork = false;
+    }
+  }
+
+  // Get the index prior to the current time
+  while (fismData.times[fism_prev_index + 1] <= time_now)
+    fism_prev_index ++;
+
+  // Determine time-interpolation weighting factor
+  precision_t dt_fism;
+  dt_fism = fismData.times[fism_prev_index + 1] - fismData.times[fism_prev_index];
+  precision_t x = (time_now - fismData.times[fism_prev_index]) / dt_fism;
+
+  // store the wavelength:
+  for (int iWave = 0; iWave < nWavelengths; iWave ++)
+    wavelengths_intensity_1au[iWave] =
+      (1.0 - x) * fismData.values[fism_prev_index][iWave]
+      + x * fismData.values[fism_prev_index + 1][iWave];
+
+  report.exit(function);
+  return didWork;
+}
+
 // --------------------------------------------------------------------------
 // Calculate EUVAC
 // --------------------------------------------------------------------------

src/euv.cpp

@@ -936,13 +936,22 @@ std::string Inputs::get_diffuse_auroral_model() {
 }
 
 // -----------------------------------------------------------------------
-// Return the EUV model used (EUVAC only option now)
+// Return the EUV model used (EUVAC, NEUVAC, FISM, etc.)
 // -----------------------------------------------------------------------
 
 std::string Inputs::get_euv_model() {
   return mklower(check_settings_str("Euv", "Model"));
 }
 
+
+// -----------------------------------------------------------------------
+// Return the FISM data file
+// -----------------------------------------------------------------------
+
+std::string Inputs::get_euv_fismfile(){
+  return get_setting_str("Euv", "fismFile");
+}
+
 // -----------------------------------------------------------------------
 // Return the heating efficiency of the neutrals for EUV
 // -----------------------------------------------------------------------