electrodynamics.h

// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members)
// Full license can be found in License.md
 
#ifndef INCLUDE_ELECTRODYNAMICS_H_
#define INCLUDE_ELECTRODYNAMICS_H_
 
/**************************************************************
 * This is the electrodynamics class for Aether.
 
   To use this, you have to:
   1. initialize it, which reads in the electrodynamics file (if used).
   2. set the time to use.
   3. set any indices to you need for the electrodynmics.
   4. set the magnetic latitudes (2d) you need
   5. set the magnetic local times (2d) you need
   6. call get_potential or get_eflux or ... to get the electrodynamics
      on the grid that was specified.
   2-6 can all be called as many times as you want.
 *
 **************************************************************/
 
#include <vector>
 
#include "aether.h"
 
class Electrodynamics {
 
 public:
 
  /**************************************************************
     \brief Initialize electrodynamics variables and routines
 
     This does the following:
     - initialize all variables to missing values
     - read in file if it exists
 
     \param input Need to pass Input class, so code can get info
                  about how user has configured things. Inside
                  input, the function uses .efield_model, .auroral_model,
                  .electrodynamics_file.
 
     \param report Need to pass Report class, so reporting can occur  
   **/
  Electrodynamics(Inputs input, Report &report);
 
  /**************************************************************
     \brief update the potential and aurora
 
     \param planet need info about SZA and stuff to get MLTs
     \param gGrid need grid
     \param time need current time
     \param ions Going to set the potential and aurora
 
     \param report Need to pass Report class, so reporting can occur  
   **/
  
  int update(Planets planet,
         Grid gGrid,
         Times time,
         Ions &ions,
         Report &report);
 
  /**************************************************************
     \brief used in main.cpp to ensure electrodynamics times and
      aether input times match up. Returns false if times are
      misaligned, true if they are aligned
 
     \param inputStartTime input file starting time
 
     \param inputEndTime input file ending time
   **/
 
  bool check_times(double inputStartTime, double inputEndTime);
  
  /**************************************************************
     \brief used in advance.cpp to get potential, eflux, avee
 
     \param magLat magnetic latitude
 
     \param magLocalTime magnetic local time
 
     \param report reporting
   **/
 
  std::tuple<arma_cube,
         arma_mat,
         arma_mat> get_electrodynamics(arma_cube magLat,
                       arma_cube magLocalTime,
                       Report &report);
 
  /**************************************************************
     \brief Gets interpolation indices
 
     Performs 2d interpolation over search vector to get indices
 
     \param vals the 2d array that needs indices
     \param search The vector of values to interpolate over  
   **/
 
  arma_mat get_interpolation_indices(arma_mat vals, arma_vec search);
 
  /**************************************************************
     \brief Sets time needed for electrodynamics
 
     Internally, if there is a file read, this function:
     - finds interpolation indices in time
     - interpolates the primaary quantities (pot, avee, eflux, etc) to
       the given time, creating internal states of these quantities. This
       is done for all relavant grids.
 
     \param time the time requested.
     \param report Need to pass Report class, so reporting can occur  
   **/
  
  void set_time(double time, Report &report);
 
  /**************************************************************
     \brief Sets the current grid to request data on
 
     Internally, if there is a file read, this function finds the
     interpolation indices in space for each relevant grid
 
     \param lats a 2D matrix of magnetic latitudes to interpolate to
     \param mlts a 2D matrix of magnetic local times to interpolate to
     \param report Need to pass Report class, so reporting can occur  
   **/
  void set_grid(arma_mat lats, arma_mat mlts, Report &report);
 
  /**************************************************************
     \brief Set the IMF Bx for internal usage
     \param value Value to assign to IMF Bx (nT)
   **/
  void set_imf_bx(precision_t value);
 
  /**************************************************************
     \brief Set the IMF By for internal usage
     \param value Value to assign to IMF By (nT)
   **/
  void set_imf_by(precision_t value);
 
  /**************************************************************
     \brief  Set the IMF Bz for internal usage
     \param value Value to assign to IMF Bz (nT)
   **/
  void set_imf_bz(precision_t value);
 
  /**************************************************************
     \brief Set the Solar Wind Velocity for internal usage
     \param value Value to assign to Solar wind velocity (km/s)
   **/
  void set_sw_v(precision_t value);
 
  /**************************************************************
     \brief Set the Solar Wind Density for internal usage
     \param value Value to assign to Solar Wind Density (/cc)
   **/
  void set_sw_n(precision_t value);
 
  /**************************************************************
     \brief Set the Hemispheric Power for internal usage
     \param value Value to assign to Hemispheric Power (GW)
     \param
   **/
  void set_hp(precision_t value);
 
  /**************************************************************
     \brief Set the AU index for internal usage
     \param value Value to assign to Auroral Upper Index (nT)
   **/
  void set_au(precision_t value);
 
  /**************************************************************
     \brief Set the AL index for internal usage
     \param value Value to assign to Auroral Lower Index (nT)
   **/
  void set_al(precision_t value);
 
  /**************************************************************
     \brief Set the AE Index for internal usage
     \param value Value to assign to Auroral Electrojet Index (nT)
   **/
  void set_ae(precision_t value);
 
  /**************************************************************
     \brief Set the Kp Index for internal usage
     \param value Value to assign to Kp index
   **/
  void set_kp(precision_t value);
  
  /**************************************************************
     \brief Get 2D electric potential on specified grid
 
     This function returns the electric potential on the requested
     grid at the requested time (with the requested indices, if
     applicable)
 
     This function does the following:
     - creates an empty potential matrix ok, I see to return
     - Loops through the grids in priority order calling set_values
       with the potentials in the grids
 
     \param report Need to pass Report class, so reporting can occur  
   **/
  arma_cube get_potential(arma_cube magLat,
              arma_cube magLocalTime,
              Report &report);
 
  /**************************************************************
     \brief Get 2D electron energy flux on specified grid
 
     This function returns the electron energy flux on the requested
     grid at the requested time (with the requested indices, if
     applicable)
 
     This function does the following:
     - creates an empty eflux matrix to return
     - Loops through the grids in priority order calling set_values
       with the eflux in the grids
 
     \param report Need to pass Report class, so reporting can occur  
   **/
  arma_mat get_eflux(arma_cube magLat, arma_cube magLocalTime, Report &report);
 
  /**************************************************************
     \brief Get 2D electron average energy on specified grid
 
     This function returns the electron average energy on the requested
     grid at the requested time (with the requested indices, if
     applicable)
 
     This function does the following:
     - creates an empty avee matrix to return
     - Loops through the grids in priority order calling set_values
       with the avee in the grids
 
     \param report Need to pass Report class, so reporting can occur  
   **/
  arma_mat get_avee(arma_cube magLat, arma_cube magLocalTime, Report &report);
 
  /**************************************************************
     \brief Get 2D ion energy flux on specified grid
 
     This function returns the ion energy flux on the requested
     grid at the requested time (with the requested indices, if
     applicable)
 
     This function does the following:
     - creates an empty ion eflux matrix to return
     - Loops through the grids in priority order calling set_values
       with the ion eflux in the grids
 
     \param report Need to pass Report class, so reporting can occur  
   **/
  arma_mat get_ion_eflux(Report &report);
 
  /**************************************************************
     \brief Get 2D ion average energy on specified grid
 
     This function returns the ion average energy on the requested
     grid at the requested time (with the requested indices, if
     applicable)
 
     This function does the following:
     - creates an empty ion avee matrix to return
     - Loops through the grids in priority order calling set_values
       with the ion avee in the grids
 
     \param report Need to pass Report class, so reporting can occur  
   **/
  arma_mat get_ion_avee(Report &report);
  
 private:
 
  int iTimeInterpolationMethod;
  const int iPrevious_ = 1;
  const int iNext_ = 2;
  // Use the closest value:
  const int iClosest_ = 3; 
  const int iInterp_ = 4;
 
  std::string input_file;
 
  bool HaveElectrodynamics;
 
 
  double time_needed;
 
  arma_mat lats_needed;
 
  arma_mat mlts_needed;
  
  precision_t imf_bx_needed;
  precision_t imf_by_needed;
  precision_t imf_bz_needed;
  precision_t sw_v_needed;
  precision_t sw_n_needed;
  precision_t hp_needed;
  precision_t au_needed;
  precision_t al_needed;
  precision_t ae_needed;
  precision_t kp_needed;
 
  int iUseGridBasedModel;
 
  std::string efield_model_to_use;
 
  std::string auroral_model_to_use;
  
  precision_t time_index;
 
  /**************************************************************
   * input_electrodynamics_struct is a structure that contains
   * information about the electrodynamics on a magnetic latitude /
   * local time grid. It is assumed that the quantities are on a
   * regular grid, in that the lat and mlt can be described with 1D
   * arrays and the quantities can be described with 2D arrays.
   * Time then makes these matrices into vectors. The quantities
   * stored are energy flux, average energy, electric potential,
   * ion energy flux, ion average energy.
   **/
  struct input_electrodynamics_struct {
 
    int nLats;
    int nMlts;
 
    arma_vec mlats;
    arma_vec mlts;
 
    std::vector<double> times;
 
    std::vector<arma_mat> potential;
 
    arma_mat potential_current;
    
    std::vector<arma_mat> energy_flux;
    arma_mat energy_flux_current;
    
    std::vector<arma_mat> average_energy;
    arma_mat average_energy_current;
    
    std::vector<arma_mat> ion_energy_flux;
    arma_mat ion_energy_flux_current;
    
    std::vector<arma_mat> ion_average_energy;
    arma_mat ion_average_energy_current;
 
    int DoesIncludeIonPrecip;
    
    int grid_priority;
 
    arma_mat lats_indices;
    arma_mat mlts_indices;
    
  };
  
  std::vector<input_electrodynamics_struct> input_electrodynamics;
  
  std::vector<int> grid_order;
  
  int nElectrodynamicsGrids;
  
  /**************************************************************
     \brief Reads a netcdf file that has the electrodynamics specification
 
     Reads a netcdf file that contains at a minimum:
     - Potential
     - Auroral energy flux
     - Auroral average energy
     May contain:
     - Ion energy flux
     - Ion average energy
 
     These are on a magnetic lat/mlt grid as a function of time, and
     should be put into the input_electrodynamics structure.
 
     \param filename
     \param report Need to pass Report class, so reporting can occur
   **/
  void read_netcdf_electrodynamics_file(std::string filename,
                    Report &report);
 
  /**************************************************************
     \brief Takes the pot/eflux/avee/etc and interpolates to the grid
 
     This function takes values from the input_electrodynamics
     structure and interpolates the values onto the user specified
     grid. The tricky bit is that there can be multiple overlapping
     grids.  So, this needs to be called for each of the existing
     grids, so that the values are overwritten. To keep it
     "functional", we pass in the last round of values and those are
     moved into the output values and then the overlapping region is
     overwritten (e.g., in the get_potential function, the 
     grids need to be cycled through calling get_values with the
     potential on that grid and the interpolation indices for the grid.
 
     \param values_current the pot/eflux/avee/etc from
     input_electrodynamics grid
 
     \param lats_indices the interpolation indices for the current
     grid latitudes
 
     \param mlts_indices the interpolation indices for the current
     grid mlts
 
     \param values_old the output of this function for the last grid
  **/  
  
  arma_mat get_values(arma_mat matToInterpolateOn, int rows, int cols);
};
 
#endif // INCLUDE_ELECTRODYNAMICS_H_