diff --git a/edu/examples/arm.cpp b/edu/examples/arm.cpp index 12b19f9c..cd272b6b 100644 --- a/edu/examples/arm.cpp +++ b/edu/examples/arm.cpp @@ -10,7 +10,7 @@ is each to use these 3d arrys. Here is how I compiled the code (clearly, I installed it in /usr/local): - g++ arm.cpp -o arm.exe -std=c++11 -O2 -larmadillo -I/usr/local/include -L/usr/local/lib + g++ arm.cpp -o arm.exe -std=c++11 -O2 -larmadillo -I/usr/local/include */ diff --git a/include/advance.h b/include/advance.h index 5504fcf1..19f7dcc4 100644 --- a/include/advance.h +++ b/include/advance.h @@ -1,8 +1,8 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_ADVANCE_H_ -#define AETHER_INCLUDE_ADVANCE_H_ +#ifndef INCLUDE_ADVANCE_H_ +#define INCLUDE_ADVANCE_H_ #include "../include/times.h" #include "../include/inputs.h" @@ -12,15 +12,15 @@ #include "../include/planets.h" #include "../include/ions.h" -int advance( Planets &planet, - Grid &gGrid, - Times &time, - Euv &euv, - Neutrals &neutrals, - Ions &ions, - Chemistry &chemistry, - Indices &indices, - Inputs &args, - Report &report); +int advance(Planets &planet, + Grid &gGrid, + Times &time, + Euv &euv, + Neutrals &neutrals, + Ions &ions, + Chemistry &chemistry, + Indices &indices, + Inputs &args, + Report &report); -#endif // AETHER_INCLUDE_ADVANCE_H_ +#endif // INCLUDE_ADVANCE_H_ diff --git a/include/bfield.h b/include/bfield.h index a6dd95a7..3a3c99e1 100644 --- a/include/bfield.h +++ b/include/bfield.h @@ -1,8 +1,12 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_BFIELD_H_ -#define AETHER_INCLUDE_BFIELD_H_ +#ifndef INCLUDE_BFIELD_H_ +#define INCLUDE_BFIELD_H_ + +#include "../include/planets.h" +#include "../include/inputs.h" +#include "../include/report.h" struct bfield_info_type { float b[3]; @@ -11,25 +15,24 @@ struct bfield_info_type { }; bfield_info_type get_bfield(float lon, - float lat, - float alt, - Planets planet, - Inputs input, - Report &report); + float lat, + float alt, + Planets planet, + Inputs input, + Report &report); bfield_info_type get_dipole(float lon, - float lat, - float alt, - Planets planet, - Inputs input, - Report &report); + float lat, + float alt, + Planets planet, + Inputs input, + Report &report); bfield_info_type get_dipole(float lon, - float lat, - float alt, - Planets planet, - Inputs input, - Report &report); - -#endif // AETHER_INCLUDE_EUV_H_ + float lat, + float alt, + Planets planet, + Inputs input, + Report &report); +#endif // INCLUDE_EUV_H_ diff --git a/include/calc_chemistry.h b/include/calc_chemistry.h index 41f39deb..fa136989 100644 --- a/include/calc_chemistry.h +++ b/include/calc_chemistry.h @@ -1,61 +1,49 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_CALC_CHEMISTRY_H_ -#define AETHER_INCLUDE_CALC_CHEMISTRY_H_ +#ifndef INCLUDE_CALC_CHEMISTRY_H_ +#define INCLUDE_CALC_CHEMISTRY_H_ + +#include +#include +#include "../include/ions.h" +#include "../include/neutrals.h" +#include "../include/times.h" +#include "../include/grid.h" +#include "../include/report.h" class Chemistry { public: Chemistry(Inputs args, Report report); - - struct reaction_type { + struct reaction_type { // Reactions: // loss1 + loss2 + loss3 -> source1 + source2 + source3 - + std::vector sources_names; std::vector losses_names; - + std::vector sources_ids; std::vector sources_IsNeutrals; - + std::vector losses_ids; std::vector losses_IsNeutrals; float rate; - }; std::vector reactions; - - struct sources_and_losses_type { + struct sources_and_losses_type { float neutral_sources[nSpecies]; float neutral_losses[nSpecies]; float ion_sources[nIons]; float ion_losses[nIons]; - float heat_neutrals; float heat_ions; float heat_electrons; - }; - void calc_chemistry(Neutrals &neutrals, - Ions &ions, - Times time, - Grid grid, - Report report); - - void calc_chemical_sources(float neutral_density, - float ion_density, - float Tn, - float Ti, - float Te, - sources_and_losses_type &sources_and_losses; - Report report); - - -#endif // AETHER_INCLUDE_CALC_CHEMISTRY_H_ +#endif // INCLUDE_CALC_CHEMISTRY_H_ diff --git a/include/calc_euv.h b/include/calc_euv.h index de0c6baf..0e8681e1 100644 --- a/include/calc_euv.h +++ b/include/calc_euv.h @@ -1,8 +1,8 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_CALC_EUV_H_ -#define AETHER_INCLUDE_CALC_EUV_H_ +#ifndef INCLUDE_CALC_EUV_H_ +#define INCLUDE_CALC_EUV_H_ #include #include @@ -20,14 +20,14 @@ // // ------------------------------------------------------------------------- -int calc_euv( Planets planet, - Grid grid, - Times time, - Euv &euv, - Neutrals &neutrals, - Ions &ions, - Indices indices, - Inputs args, - Report &report); +int calc_euv(Planets planet, + Grid grid, + Times time, + Euv &euv, + Neutrals &neutrals, + Ions &ions, + Indices indices, + Inputs args, + Report &report); -#endif // AETHER_INCLUDE_CALC_EUV_H_ +#endif // INCLUDE_CALC_EUV_H_ diff --git a/include/chemistry.h b/include/chemistry.h index 4654a897..258ad189 100644 --- a/include/chemistry.h +++ b/include/chemistry.h @@ -1,8 +1,8 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_CHEMISTRY_H_ -#define AETHER_INCLUDE_CHEMISTRY_H_ +#ifndef INCLUDE_CHEMISTRY_H_ +#define INCLUDE_CHEMISTRY_H_ #include #include @@ -13,55 +13,47 @@ #include "../include/inputs.h" #include "../include/report.h" - class Chemistry { public: struct reaction_type { - // Reactions: // loss1 + loss2 + loss3 -> source1 + source2 + source3 - std::vector sources_names; std::vector losses_names; - + std::vector sources_ids; std::vector sources_IsNeutral; - + std::vector losses_ids; std::vector losses_IsNeutral; int nSources; int nLosses; - + float energy; float rate; float branching_ratio; - }; std::vector reactions; - long nReactions; - + int64_t nReactions; + Chemistry(Neutrals neutrals, - Ions ions, - Inputs args, - Report &report); + Ions ions, + Inputs args, + Report &report); - void calc_chemistry(Neutrals &neutrals, - Ions &ions, - Times time, - Grid grid, - Report &report); - - void calc_chemical_sources(float *neutral_density, - float *ion_density, - float Tn, - float Ti, - float Te, - Report &report); + Ions &ions, + Times time, + Grid grid, + Report &report); + + void calc_chemical_sources(Neutrals &neutrals, + Ions &ions, + Report &report); private: @@ -71,36 +63,32 @@ class Chemistry { float neutral_losses[nSpecies]; float ion_sources[nIons]; float ion_losses[nIons]; - + float heat_neutrals; float heat_ions; float heat_electrons; - }; sources_and_losses_type sources_and_losses; - + int read_chemistry_file(Neutrals neutrals, - Ions ions, - Inputs args, - Report &report); - + Ions ions, + Inputs args, + Report &report); + reaction_type interpret_reaction_line(Neutrals neutrals, - Ions ions, - std::vector line, - Report &report); + Ions ions, + std::vector line, + Report &report); void find_species_id(std::string name, - Neutrals neutrals, - Ions ions, - int &id_, - int &IsNeutral, - Report &report); + Neutrals neutrals, + Ions ions, + int &id_, + int &IsNeutral, + Report &report); void display_reaction(reaction_type reaction); - - - }; -#endif // AETHER_INCLUDE_CHEMISTRY_H_ +#endif // INCLUDE_CHEMISTRY_H_ diff --git a/include/constants.h b/include/constants.h index 35185f9e..070c6ebf 100644 --- a/include/constants.h +++ b/include/constants.h @@ -1,36 +1,36 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_CONSTANTS_H_ -#define AETHER_INCLUDE_CONSTANTS_H_ +#ifndef INCLUDE_CONSTANTS_H_ +#define INCLUDE_CONSTANTS_H_ #include // Header file to define a bunch of physical constants -const float avogadros_number = 6.02214086e23; // -const float boltzmanns_constant = 1.38064852e-23; // m2 kg /s2 /K -const float mass_proton = 1.6726219e-27; // kg +const float avogadros_number = 6.02214086e23; // per mole +const float boltzmanns_constant = 1.38064852e-23; // m2 kg /s2 /K +const float mass_proton = 1.6726219e-27; // kg const float amu = mass_proton; -const float mass_electron = 9.1094e-31; // Kg -const float planck_constant = 6.6261e-34; // Js -const float element_charge = 1.6022e-19; // C (J/eV) -const float speed_light = 2.9979e8; // m/s -const float univ_gas_constant = avogadros_number*boltzmanns_constant; // J/(moleK) -const float r_gas = univ_gas_constant*1.0E+07; // erg/(moleK) -const float pi = 3.141592653589793; +const float mass_electron = 9.1094e-31; // Kg +const float planck_constant = 6.6261e-34; // Js +const float element_charge = 1.6022e-19; // C (J/eV) +const float speed_light = 2.9979e8; // m/s +const float univ_gas_constant = avogadros_number*boltzmanns_constant; +const float r_gas = univ_gas_constant*1.0E+07; // erg/(moleK) +const float pi = 3.141592653589793; const float twopi = 2*pi; const float dtor = pi/180.0; const float rtod = 180.0/pi; -const float p00 = 1.0e5; // Pa -const float T00 = 273.0; // K -const float SpecificHeatVolume = 3./2.; +const float p00 = 1.0e5; // Pa +const float T00 = 273.0; // K +const float SpecificHeatVolume = 3./2.; const float SpecificHeatPressure = SpecificHeatVolume + 1.0; const float gamma_const = SpecificHeatPressure/SpecificHeatVolume; // A bunch of constants for converting time between seconds and arrays: -const std::vector days_of_month {31,28,31,30,31,30,31,31,30,31,30,31}; -const int reference_year = 1965; // I think that this has to be the year after a leap year. +const std::vector days_of_month {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}; +const int reference_year = 1965; const double julian_day_of_reference = 2438762.0; const double j2000 = 2451545.0; const double seconds_per_year = 31536000.0; @@ -41,15 +41,15 @@ const double seconds_per_minute = 60.0; // Conversion const float mtokm = 1.0/1000.0; const float kmtom = 1000.0; -const float pcm3topm3 = 1e6; // /cm3 to /m3 -const float pcm2topm2 = 1e4; // /cm2 to /m2 -const float pcmtopm = 100.0; // /cm to /m -const float atom = 1.0e-10; // angstrom to m +const float pcm3topm3 = 1e6; // /cm3 to /m3 +const float pcm2topm2 = 1e4; // /cm2 to /m2 +const float pcmtopm = 100.0; // /cm to /m +const float atom = 1.0e-10; // angstrom to m // Stellar constants: -const float au_to_m = 1.495978707e11; // m -const float gravitational_constant = 6.67408e-11; // m3/kg/s2 +const float au_to_m = 1.495978707e11; // m +const float gravitational_constant = 6.67408e-11; // m3/kg/s2 // This is for our sun: -const float solar_mass = 1.989e30; // kg +const float solar_mass = 1.989e30; // kg -#endif // AETHER_INCLUDE_CONSTANTS_H_ +#endif // INCLUDE_CONSTANTS_H_ diff --git a/include/earth.h b/include/earth.h index 9b2ef474..5f3603db 100644 --- a/include/earth.h +++ b/include/earth.h @@ -1,3 +1,5 @@ +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) +// Full license can be found in License.md #define iO_ 0 #define iO2_ 1 @@ -22,4 +24,3 @@ #define iOP_2P_ 7 #define iElec_ 8 #define nIons 8 - diff --git a/include/euv.h b/include/euv.h index ca5f4896..992f7c79 100644 --- a/include/euv.h +++ b/include/euv.h @@ -1,8 +1,8 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_EUV_H_ -#define AETHER_INCLUDE_EUV_H_ +#ifndef INCLUDE_EUV_H_ +#define INCLUDE_EUV_H_ #include #include @@ -23,14 +23,12 @@ class Euv { int nLines; struct waveinfotype { - std::string name; std::string to; std::string type; std::string units; std::string note; std::vector values; - }; std::vector waveinfo; @@ -61,9 +59,9 @@ class Euv { // ------------------------------------------------------------------------- int scale_from_1au(Planets planet, Times time); - + private: - + // -------------------------------------------------------------------------- // Read in the EUV file that describes all of the wavelengths and // cross sections @@ -76,10 +74,9 @@ class Euv { // -------------------------------------------------------------------------- int slot_euv(std::string item, - std::string item2, - std::vector &values, - Report report); - + std::string item2, + std::vector &values, + Report report); }; -#endif // AETHER_INCLUDE_EUV_H_ +#endif // INCLUDE_EUV_H_ diff --git a/include/file_input.h b/include/file_input.h index 96a331f1..3c2d0f6f 100644 --- a/include/file_input.h +++ b/include/file_input.h @@ -1,8 +1,8 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_FILE_INPUT_H_ -#define AETHER_INCLUDE_FILE_INPUT_H_ +#ifndef INCLUDE_FILE_INPUT_H_ +#define INCLUDE_FILE_INPUT_H_ #include #include @@ -17,4 +17,4 @@ std::string strip_spaces(std::string instring); std::string read_string(std::ifstream &file_ptr, std::string hash); int read_int(std::ifstream &file_ptr, std::string hash); -#endif // AETHER_INCLUDE_FILE_INPUT_H_ +#endif // INCLUDE_FILE_INPUT_H_ diff --git a/include/fill_grid.h b/include/fill_grid.h index 936070bd..8bed9723 100644 --- a/include/fill_grid.h +++ b/include/fill_grid.h @@ -1,8 +1,8 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_INIT_FILL_GRID_H_ -#define AETHER_INCLUDE_INIT_FILL_GRID_H_ +#ifndef INCLUDE_FILL_GRID_H_ +#define INCLUDE_FILL_GRID_H_ #include "inputs.h" #include "planets.h" @@ -10,7 +10,5 @@ #include "times.h" void fill_grid_radius(Grid &gGrid, Planets planet, Inputs input); -// void fill_grid_sza(Grid &gGrid, Planets planet, Times time, Inputs input); -// void fill_grid(Grid &gGrid, Planets planet, Inputs input); -#endif // AETHER_INCLUDE_INIT_FILL_GRID_H_ +#endif // INCLUDE_FILL_GRID_H_ diff --git a/include/grid.h b/include/grid.h index c3f53c29..6254ce40 100644 --- a/include/grid.h +++ b/include/grid.h @@ -1,14 +1,20 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_GRID_H_ -#define AETHER_INCLUDE_GRID_H_ +#ifndef INCLUDE_GRID_H_ +#define INCLUDE_GRID_H_ + +// The armadillo library is to allow the use of 3d cubes and other +// array types, with array math built in. This eliminates loops! +#include #include "inputs.h" #include "sizes.h" #include "planets.h" #include "times.h" +using namespace arma; + // We need a naming convention for the variables that are defined on // the grid. These could then match the formulas that are used to find // the given points on the grid. @@ -18,13 +24,20 @@ 1 - indication of whether variable is scalar (s) or vector (v) 2 - physical dimensions: 3 - 3d (e.g., lon, lat, alt or x, y, z) - - + c - added this to indicate that it is an armadillo "cube" + In theory, there should be a 2 and 1, but I have not created any examples + 3 - Whether there are ghost cells or not + g - ghostcells + n - no ghost cells. I have not created any of these yet. + 4 - Whether the variable is defined at the center or edge + c - center + e - edge For example: _s3gc : scalar variable, 3d, include ghost cells, cell centers - _31ne : - + _vcne : vector variable, cube (3d armadillo), no ghost cells, cell edges + _scgc : most common type, scalar, cube, ghost cells, cell centers + _vcgc : vector, cube, ghost, centers: vector */ @@ -35,33 +48,36 @@ For example: // Scalars, 3D, Include Ghostcells, Cell Centers: #define ijk_geo_s3gc(i,j,k) \ - ((i)*long(nGeoLatsG)*long(nGeoAltsG) + \ - (j)*long(nGeoAltsG) + \ + ((i)*int64_t(nGeoLatsG)*int64_t(nGeoAltsG) + \ + (j)*int64_t(nGeoAltsG) + \ (k)) #define ijk_mag_s3gc(i,j,k) \ - ((i)*long(nMagLatsG)*long(nGeoAltsG) + \ - (j)*long(nMagAltsG) + \ + ((i)*int64_t(nMagLatsG)*int64_t(nGeoAltsG) + \ + (j)*int64_t(nMagAltsG) + \ (k)) // Scalars, 3D, Include Ghostcells, Cell Edges (Altitude): #define ijk_geo_s3ge3(i,j,k) \ - ((i)*long(nGeoLatsG)*long(nGeoAltsG+1) + \ - (j)*long(nGeoAltsG+1) + \ + ((i)*int64_t(nGeoLatsG)*int64_t(nGeoAltsG+1) + \ + (j)*int64_t(nGeoAltsG+1) + \ (k)); // Vectors, 3D, Include Ghostcells, Cell Centers: #define ijkl_geo_v3gc(i,j,k,l) \ - ((i)*long(nGeoLatsG)*long(nGeoAltsG)*long(3) + \ - (j)*long(nGeoAltsG)*long(3) + \ - (k)*long(3) + \ + ((i)*int64_t(nGeoLatsG)*int64_t(nGeoAltsG)*int64_t(3) + \ + (j)*int64_t(nGeoAltsG)*int64_t(3) + \ + (k)*int64_t(3) + \ (l)) // Vectors, 3D, Include Ghostcells, Cell Centers: #define ijkl_mag_v3gc(i,j,k,l) \ - ((i)*long(nMagLatsG)*long(nMagAltsG)*long(3) + \ - (j)*long(nMagAltsG)*long(3) + \ - (k)*long(3) + \ + ((i)*int64_t(nMagLatsG)*int64_t(nMagAltsG)*int64_t(3) + \ + (j)*int64_t(nMagAltsG)*int64_t(3) + \ + (k)*int64_t(3) + \ (l)) +// ---------------------------------------------------------------------------- +// Grid class +// ---------------------------------------------------------------------------- class Grid { @@ -70,18 +86,29 @@ class Grid { int get_IsGeoGrid(); void set_IsGeoGrid(int value); - long get_nPointsInGrid(); - - // These define the geographic grid: - float *geoLon_s3gc, *geoX_s3gc; - float *geoLat_s3gc, *geoY_s3gc; - float *geoAlt_s3gc, *geoZ_s3gc; - + int64_t get_nPointsInGrid(); + + int64_t get_nX(); + int64_t get_nY(); + int64_t get_nZ(); + + int64_t get_nLons(); + int64_t get_nLats(); + int64_t get_nAlts(); + + int64_t get_nGCs(); + + // Armidillo Cube Versions: + fcube geoLon_scgc, geoX_scgc; + fcube geoLat_scgc, geoY_scgc; + fcube geoAlt_scgc, geoZ_scgc; + // These define the magnetic grid: - float *magLon_s3gc, *magX_s3gc; - float *magLat_s3gc, *magY_s3gc; - float *magAlt_s3gc, *magZ_s3gc; - float *magLocalTime_s3gc; + // Armidillo Cube Versions: + fcube magLon_scgc, magX_scgc; + fcube magLat_scgc, magY_scgc; + fcube magAlt_scgc, magZ_scgc; + fcube magLocalTime_scgc; std::string altitude_name = "Altitude"; std::string altitude_unit = "meters"; @@ -91,25 +118,25 @@ class Grid { std::string latitude_name = "Latitude"; std::string latitude_unit = "radians"; - + // These are derived variables from the grid: - float *radius_s3gc; - float *radius_sq_s3gc; - float *radius_inv_sq_s3gc; - float *gravity_s3gc; - float *sza_s3gc, *cos_sza_s3gc; - - float *bfield_v3gc; - float *bfield_mag_s3gc; - - float *geoX_edges_s3ge, *geoX_cell_width_s3gc; - - float *alt_edges; - float *alt_cell_width; - float *dalt_center_s3gc; - float *dalt_lower_s3gc; - - Grid(int nX, int nY, int nZ); + + // Switch to armadillo variables (float cubes): + fcube radius_scgc; + fcube radius2_scgc; + fcube radius2i_scgc; + fcube gravity_scgc; + + fcube sza_scgc; + fcube cos_sza_scgc; + + fcube dalt_center_scgc; + fcube dalt_lower_scgc; + + std::vector bfield_vcgc; + fcube bfield_mag_scgc; + + Grid(int nX_in, int nY_in, int nZ_in, int nGCs_in); void calc_sza(Planets planet, Times time, Report &report); void fill_grid(Planets planet, Report &report); @@ -121,6 +148,12 @@ class Grid { int IsGeoGrid; + int64_t nX, nLons; + int64_t nY, nLats; + int64_t nZ, nAlts; + + int nGCs; // number of ghostcells + }; -#endif // AETHER_INCLUDE_GRID_H_ +#endif // INCLUDE_GRID_H_ diff --git a/include/indices.h b/include/indices.h index 990571c1..d0ddf244 100644 --- a/include/indices.h +++ b/include/indices.h @@ -1,14 +1,12 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_INDICES_H_ -#define AETHER_INCLUDE_INDICES_H_ +#ifndef INCLUDE_INDICES_H_ +#define INCLUDE_INDICES_H_ #include #include - - class Indices { public: @@ -23,10 +21,8 @@ class Indices { private: struct ind_time_pair { - float index; double time; - }; std::vector f107; @@ -37,11 +33,10 @@ class Indices { // This is the general method for setting indices: int set_index(std::vector time, - std::vector indexarray, - std::vector &index); + std::vector indexarray, + std::vector &index); float get_index(double time, std::vector index); - }; -#endif // AETHER_INCLUDE_INDICES_H_ +#endif // INCLUDE_INDICES_H_ diff --git a/include/init_geo_grid.h b/include/init_geo_grid.h index d7deee26..f0077253 100644 --- a/include/init_geo_grid.h +++ b/include/init_geo_grid.h @@ -1,8 +1,8 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_INIT_GEO_GRID_H_ -#define AETHER_INCLUDE_INIT_GEO_GRID_H_ +#ifndef INCLUDE_INIT_GEO_GRID_H_ +#define INCLUDE_INIT_GEO_GRID_H_ #include "inputs.h" #include "planets.h" @@ -10,5 +10,5 @@ void init_geo_grid(Grid &gGrid, Planets planet, Inputs input); -#endif // AETHER_INCLUDE_INIT_GEO_GRID_H_ +#endif // INCLUDE_INIT_GEO_GRID_H_ diff --git a/include/inputs.h b/include/inputs.h index 99a81fb7..e4764c9f 100644 --- a/include/inputs.h +++ b/include/inputs.h @@ -1,8 +1,8 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_INPUTS_H_ -#define AETHER_INCLUDE_INPUTS_H_ +#ifndef INCLUDE_INPUTS_H_ +#define INCLUDE_INPUTS_H_ #include #include @@ -30,25 +30,23 @@ class Inputs { std::string get_planetary_file(); std::string get_planet_species_file(); std::string get_bfield_type(); - + // ------------------------------ // Grid inputs: struct grid_input_struct { - std::string alt_file; int IsUniformAlt; float alt_min; float dalt; - float lat_min; float lat_max; float lon_min; float lon_max; }; - grid_input_struct get_grid_inputs(); - + grid_input_struct get_grid_inputs(); + int iVerbose; private: @@ -63,16 +61,15 @@ class Inputs { std::string planet_species_file = ""; std::string bfield = "none"; - + grid_input_struct grid_input; - + float euv_heating_eff_neutrals; float euv_heating_eff_electrons; std::vector dt_output; std::vector type_output; float dt_euv; - }; -#endif // AETHER_INCLUDE_INPUTS_H_ +#endif // INCLUDE_INPUTS_H_ diff --git a/include/ions.h b/include/ions.h index 7b7ebd62..72693f1c 100644 --- a/include/ions.h +++ b/include/ions.h @@ -1,53 +1,68 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_IONS_H_ -#define AETHER_INCLUDE_IONS_H_ +#ifndef INCLUDE_IONS_H_ +#define INCLUDE_IONS_H_ + +#include +#include + +// The armadillo library is to allow the use of 3d cubes and other +// array types, with array math built in. This eliminates loops! +#include #include "inputs.h" #include "report.h" #include "grid.h" +using namespace arma; + class Ions { public: - struct species_chars { + // This struct contains all of the information needed for a single + // species of ion. We will then have a vector of these species. + struct species_chars { std::string cName; float mass; int charge; - + int DoAdvect; - - float *density_s3gc; - float *par_velocity_v3gc; - float *perp_velocity_v3gc; - float *temperature_s3gc; - // Sources and Losses: - float *ionization_s3gc; - + fcube density_scgc; + fcube par_velocity_vcgc; + fcube perp_velocity_vcgc; + + fcube temperature_scgc; + + // Sources and Losses: + + fcube ionization_scgc; + + fcube sources_scgc; + fcube losses_scgc; }; // bulk quantities (states): - float *density_s3gc; - float *velocity_v3gc; - float *exb_v3gc; - float *ion_temperature_s3gc; - float *electron_temperature_s3gc; + fcube density_scgc; + fcube ion_temperature_scgc; + fcube electron_temperature_scgc; + + // This is the vector that will contain all of the different species: std::vector species; - + // ------------------------------ // Functions: - - Ions(Inputs input, Report report); - species_chars create_species(); + + Ions(Grid grid, Inputs input, Report report); + species_chars create_species(Grid grid); int read_planet_file(Inputs input, Report report); - void fill_electrons(Grid grid, Report &report); + void fill_electrons(Report &report); }; -#endif // AETHER_INCLUDE_NEUTRALS_H_ +#endif // INCLUDE_IONS_H_ diff --git a/include/neutrals.h b/include/neutrals.h index 2ee95de0..a25b4a79 100644 --- a/include/neutrals.h +++ b/include/neutrals.h @@ -1,8 +1,15 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_NEUTRALS_H_ -#define AETHER_INCLUDE_NEUTRALS_H_ +#ifndef INCLUDE_NEUTRALS_H_ +#define INCLUDE_NEUTRALS_H_ + +#include +#include + +// The armadillo library is to allow the use of 3d cubes and other +// array types, with array math built in. This eliminates loops! +#include #include "grid.h" #include "euv.h" @@ -11,20 +18,23 @@ #include "inputs.h" #include "report.h" +using namespace arma; + class Neutrals { public: + // This struct contains all of the information needed for a single + // species of neutrals. We will then have a vector of these species. + struct species_chars { - std::string cName; float mass; float vibe; int DoAdvect; - - float *density_s3gc; - float *velocity_v3gc; + + fcube density_scgc; std::vector diff0; std::vector diff_exp; @@ -39,31 +49,34 @@ class Neutrals { std::vector iEuvIonId_; // Some derived quantities: - float *chapman_s3gc; + fcube chapman_scgc; + fcube scale_height_scgc; // Sources and Losses: - float *ionization_s3gc; + fcube ionization_scgc; - // If we want a flat lower BC: + fcube sources_scgc; + fcube losses_scgc; + + // If we want a fixed lower BC: float lower_bc_density; - }; - + // bulk quantities (states): - float *density_s3gc; - float *velocity_v3gc; - float *temperature_s3gc; - - float *rho_s3gc; - float *mean_major_mass_s3gc; - float *pressure_s3gc; - float *sound_s3gc; - + + fcube density_scgc; + fcube temperature_scgc; + + fcube rho_scgc; + fcube mean_major_mass_scgc; + fcube pressure_scgc; + fcube sound_scgc; + // For heating/cooling: - float *Cv_s3gc; - float *gamma_s3gc; - float *kappa_s3gc; + fcube Cv_scgc; + fcube gamma_scgc; + fcube kappa_scgc; std::vector neutrals; @@ -71,36 +84,33 @@ class Neutrals { // Source terms: - float *heating_euv_s3gc; - float *conduction_s3gc; + fcube conduction_scgc; + fcube heating_euv_scgc; float heating_efficiency; - + // This is an initial temperature profile, read in through the // planet.in file: float *initial_temperatures, *initial_altitudes; - int nInitial_temps=0; + int nInitial_temps = 0; // names and units - std::string density_unit="(/m3)"; - std::string density_name="Neutral Bulk Density"; + std::string density_unit = "(/m3)"; + std::string density_name = "Neutral Bulk Density"; - std::string velocity_unit="(m/s)"; + std::string velocity_unit = "(m/s)"; std::vector velocity_name; - std::string temperature_unit="(K)"; - std::string temperature_name="Temperature"; - + std::string temperature_unit = "(K)"; + std::string temperature_name = "Temperature"; + // ------------------------------ // Functions: - + Neutrals(Grid grid, Inputs input, Report report); - species_chars create_species(); + species_chars create_species(Grid grid); int read_planet_file(Inputs input, Report report); int initial_conditions(Grid grid, Inputs input, Report report); - float calc_scale_height(int iSpecies, - long index, - Grid grid); int pair_euv(Euv euv, Ions ions, Report report); void calc_mass_density(Report &report); void calc_specific_heat(Report &report); @@ -108,10 +118,8 @@ class Neutrals { void calc_ionization_heating(Euv euv, Ions &ions, Report &report); void calc_conduction(Grid grid, Times time, Report &report); void add_sources(Times time, Report &report); - + void set_bcs(Report &report); }; - - -#endif // AETHER_INCLUDE_NEUTRALS_H_ +#endif // INCLUDE_NEUTRALS_H_ diff --git a/include/output.h b/include/output.h index 109438f7..ac1c855c 100644 --- a/include/output.h +++ b/include/output.h @@ -1,8 +1,8 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_OUTPUT_H_ -#define AETHER_INCLUDE_OUTPUT_H_ +#ifndef INCLUDE_OUTPUT_H_ +#define INCLUDE_OUTPUT_H_ #include "../include/neutrals.h" #include "../include/grid.h" @@ -12,12 +12,11 @@ #include "../include/report.h" int output(Neutrals neutrals, - Ions ions, - Grid grid, - Times time, - Planets planet, - Inputs args, - Report &report); + Ions ions, + Grid grid, + Times time, + Planets planet, + Inputs args, + Report &report); - -#endif // AETHER_INCLUDE_OUTPUT_H_ +#endif // INCLUDE_OUTPUT_H_ diff --git a/include/planets.h b/include/planets.h index 49844e9e..f57819b7 100644 --- a/include/planets.h +++ b/include/planets.h @@ -1,4 +1,4 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #ifndef AETHER_INCLUDE_PLANETS_H_ @@ -54,9 +54,10 @@ class Planets { float rates_nodelongitude; float planet_tilt; + // rotation_period, omega, and length of day are all related to each other float rotation_period; - float omega; // rotation frequency - float length_of_day; // this and rotation period are related to each other + float omega; + float length_of_day; double length_of_year; float longitude_jb2000; float longitude_offset; @@ -81,7 +82,6 @@ class Planets { float ls; double update_time; - }; std::vector planets; @@ -89,7 +89,6 @@ class Planets { planet_chars planet; int read_file(Inputs args, Report report); - }; -#endif // AETHER_INCLUDE_PLANETS_H_ +#endif // INCLUDE_PLANETS_H_ diff --git a/include/read_f107_file.h b/include/read_f107_file.h index 4c8d6cfe..aeb538ab 100644 --- a/include/read_f107_file.h +++ b/include/read_f107_file.h @@ -1,14 +1,14 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_READ_F107_FILE_H_ -#define AETHER_INCLUDE_READ_F107_FILE_H_ +#ifndef INCLUDE_READ_F107_FILE_H_ +#define INCLUDE_READ_F107_FILE_H_ #include #include int read_f107_file(std::string f107_file, - std::vector &time, - std::vector &f107); + std::vector &time, + std::vector &f107); -#endif // AETHER_INCLUDE_READ_F107_FILE_H_ +#endif // INCLUDE_READ_F107_FILE_H_ diff --git a/include/report.h b/include/report.h index 7c537ddc..d589d502 100644 --- a/include/report.h +++ b/include/report.h @@ -1,8 +1,8 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_REPORT_H_ -#define AETHER_INCLUDE_REPORT_H_ +#ifndef INCLUDE_REPORT_H_ +#define INCLUDE_REPORT_H_ #include #include @@ -24,10 +24,9 @@ class Report { private: - int iVerbose; - - struct item_struct { + int iVerbose; + struct item_struct { std::string entry; int nTimes; float timing_total; @@ -35,9 +34,8 @@ class Report { int iLevel; int iStringPosBefore; int iLastEntry; - }; - + std::vector entries; int nEntries; std::string current_entry; @@ -45,9 +43,8 @@ class Report { std::string divider; int divider_length; - - int iLevel; + int iLevel; }; -#endif // AETHER_INCLUDE_REPORT_H_ +#endif // INCLUDE_REPORT_H_ diff --git a/include/sizes.h b/include/sizes.h index 16101afe..1b1749a9 100644 --- a/include/sizes.h +++ b/include/sizes.h @@ -1,14 +1,18 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_SIZES_H_ -#define AETHER_INCLUDE_SIZES_H_ +#ifndef INCLUDE_SIZES_H_ +#define INCLUDE_SIZES_H_ + +// This is the file that defines the number of grid points in each +// direction. The entire code is based on these numbers, so you need +// to recompile if you change these numbers. // This is for the geographic grid: #define nGeoGhosts 2 -#define nGeoAlts 50 +#define nGeoAlts 60 #define nGeoAltsG nGeoGhosts + nGeoAlts + nGeoGhosts #define iGeoAltStart_ nGeoGhosts // Inclusive!!! #define iGeoAltEnd_ nGeoGhosts + nGeoAlts - 1 // Inclusive!!! @@ -47,4 +51,4 @@ #define nCharsShort 20 #define nCharsLong 400 -#endif // AETHER_INCLUDE_SIZES_H_ +#endif // INCLUDE_SIZES_H_ diff --git a/include/solvers.h b/include/solvers.h index 0728b8da..7dc22399 100644 --- a/include/solvers.h +++ b/include/solvers.h @@ -1,21 +1,24 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_SOLVERS_H_ -#define AETHER_INCLUDE_SOLVERS_H_ +#ifndef INCLUDE_SOLVERS_H_ +#define INCLUDE_SOLVERS_H_ -#include "sizes.h" +// The armadillo library is to allow the use of 3d cubes and other +// array types, with array math built in. This eliminates loops! +#include +using namespace arma; -int solver_conduction(float value[nGeoAltsG], - float lambda[nGeoAltsG], - float front[nGeoAltsG], - float dt, - float dalt_lower[nGeoAltsG], - float *conduction); +fvec solver_conduction(fvec value, + fvec lambda, + fvec front, + float dt, + fvec dx); -float solver_chemistry(float old_density, - float source, - float loss, + +fcube solver_chemistry(fcube density, + fcube source, + fcube loss, float dt); -#endif // AETHER_INCLUDE_SOLVERS_H_ +#endif // INCLUDE_SOLVERS_H_ diff --git a/include/time_conversion.h b/include/time_conversion.h index 325246f4..2952409e 100644 --- a/include/time_conversion.h +++ b/include/time_conversion.h @@ -1,8 +1,8 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_TIME_CONVERSION_H_ -#define AETHER_INCLUDE_TIME_CONVERSION_H_ +#ifndef INCLUDE_TIME_CONVERSION_H_ +#define INCLUDE_TIME_CONVERSION_H_ #include @@ -13,5 +13,5 @@ void time_real_to_int(double timereal, std::vector &itime); int test_time_routines(); void display_itime(std::vector itime); -#endif // AETHER_INCLUDE_TIME_CONVERSION_H_ +#endif // INCLUDE_TIME_CONVERSION_H_ diff --git a/include/times.h b/include/times.h index 09b5764d..5623c8ae 100644 --- a/include/times.h +++ b/include/times.h @@ -1,8 +1,8 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_TIMES_H_ -#define AETHER_INCLUDE_TIMES_H_ +#ifndef INCLUDE_TIMES_H_ +#define INCLUDE_TIMES_H_ #include #include @@ -37,7 +37,7 @@ class Times { double current; double intermediate; double simulation; - long iStep; + int64_t iStep; float dt; @@ -62,7 +62,6 @@ class Times { time_t sys_time_start; time_t sys_time_current; float walltime; - }; -#endif // AETHER_INCLUDE_TIMES_H_ +#endif // INCLUDE_TIMES_H_ diff --git a/include/transform.h b/include/transform.h index 8bceb03c..20a77508 100644 --- a/include/transform.h +++ b/include/transform.h @@ -1,32 +1,42 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#ifndef AETHER_INCLUDE_TRANSFORM_H_ -#define AETHER_INCLUDE_TRANSFORM_H_ +#ifndef INCLUDE_TRANSFORM_H_ +#define INCLUDE_TRANSFORM_H_ + +// The armadillo library is to allow the use of 3d cubes and other +// array types, with array math built in. This eliminates loops! +#include +using namespace arma; + +void copy_cube_to_array(fcube cube_in, + float *array_out); + +void transform_llr_to_xyz_3d(fcube lat3d, + fcube lon3d, + fcube r3d, + fcube &x3d, + fcube &y3d, + fcube &z3d); void transform_llr_to_xyz(float llr_in[3], float xyz_out[3]); void transform_rot_z(float xyz_in[3], float angle_in, float xyz_out[3]); void transform_rot_y(float xyz_in[3], float angle_in, float xyz_out[3]); void transform_float_vector_to_array(std::vector input, - float output[3]); + float output[3]); void transform_vector_xyz_to_env(float xyz_in[3], - float lon, - float lat, - float env_out[3]); + float lon, + float lat, + float env_out[3]); -// These are not really transformations, but are +// These are not really transformations, but are placeholders void vector_diff(float vect_in_1[3], - float vect_in_2[3], - float vect_out[3]); - -void get_vector_component(float *vector_in_v3gc, - int iComponent, - int IsGeoGrid, - float *component_out_s3gc); + float vect_in_2[3], + float vect_out[3]); -#endif // AETHER_INCLUDE_TRANSFORM_H_ +#endif // INCLUDE_TRANSFORM_H_ diff --git a/inputs/chemistry_earth.xlsx b/inputs/chemistry_earth.xlsx index 20be9491..96f0c03c 100644 Binary files a/inputs/chemistry_earth.xlsx and b/inputs/chemistry_earth.xlsx differ diff --git a/src/Makefile b/src/Makefile index 4f19a958..e826fe1b 100644 --- a/src/Makefile +++ b/src/Makefile @@ -5,8 +5,10 @@ COMPILE.CPP = g++ LINK.CPP = g++ AR = ar -rs -# FLAGS = -O3 -ffast-math -c -I/opt/local/include -FLAGS = -c -I/opt/local/include +FLAGS = -O3 -ffast-math -c -I/opt/local/include + +# Can use this version for no optimization and no fast math: +# FLAGS = -c -I/opt/local/include .SUFFICES: .SUFFICES: .cpp .o diff --git a/src/add_sources.cpp b/src/add_sources.cpp index 7fe13044..00c2d41f 100644 --- a/src/add_sources.cpp +++ b/src/add_sources.cpp @@ -1,4 +1,4 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include @@ -9,32 +9,21 @@ #include "../include/inputs.h" #include "../include/report.h" -void Neutrals::add_sources( Times time, Report &report) { +void Neutrals::add_sources(Times time, Report &report) { - std::string function="add_sources"; + std::string function = "add_sources"; static int iFunction = -1; - report.enter(function, iFunction); + report.enter(function, iFunction); - long iLon, iLat, iAlt, index; + int64_t iLon, iLat, iAlt, index; float dt = time.get_dt(); - - for (iLon = 0; iLon < nGeoLonsG; iLon++) { - for (iLat = 0; iLat < nGeoLatsG; iLat++) { - for (iAlt = 0; iAlt < nGeoAltsG; iAlt++) { - - index = ijk_geo_s3gc(iLon,iLat,iAlt); - - temperature_s3gc[index] = - temperature_s3gc[index] + - dt * ( heating_euv_s3gc[index] + - conduction_s3gc[index]); - - } - } - } - - report.exit(function); + + temperature_scgc = + temperature_scgc + + dt * (heating_euv_scgc + + conduction_scgc); + + report.exit(function); return; - } diff --git a/src/advance.cpp b/src/advance.cpp index 5b8b79c5..5868e6e3 100644 --- a/src/advance.cpp +++ b/src/advance.cpp @@ -1,4 +1,4 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include @@ -15,21 +15,20 @@ #include "../include/report.h" #include "../include/output.h" +int advance(Planets &planet, + Grid &gGrid, + Times &time, + Euv &euv, + Neutrals &neutrals, + Ions &ions, + Chemistry &chemistry, + Indices &indices, + Inputs &input, + Report &report) { -int advance( Planets &planet, - Grid &gGrid, - Times &time, - Euv &euv, - Neutrals &neutrals, - Ions &ions, - Chemistry &chemistry, - Indices &indices, - Inputs &input, - Report &report) { + int iErr = 0; - int iErr=0; - - std::string function="advance"; + std::string function = "advance"; static int iFunction = -1; report.enter(function, iFunction); @@ -39,29 +38,30 @@ int advance( Planets &planet, neutrals.calc_mass_density(report); neutrals.calc_specific_heat(report); time.calc_dt(); - + iErr = calc_euv(planet, - gGrid, - time, - euv, - neutrals, - ions, - indices, - input, - report); - + gGrid, + time, + euv, + neutrals, + ions, + indices, + input, + report); + neutrals.calc_conduction(gGrid, time, report); neutrals.add_sources(time, report); - chemistry.calc_chemistry(neutrals, ions, time, gGrid, report); - + chemistry.calc_chemistry(neutrals, ions, time, gGrid, report); + + neutrals.set_bcs(report); + time.increment_time(); iErr = output(neutrals, ions, gGrid, time, planet, input, report); report.exit(function); return iErr; - } - + diff --git a/src/bfield.cpp b/src/bfield.cpp index d1d010f5..c5a03cde 100644 --- a/src/bfield.cpp +++ b/src/bfield.cpp @@ -1,4 +1,4 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include "../include/inputs.h" @@ -8,15 +8,15 @@ #include "../include/constants.h" bfield_info_type get_bfield(float lon, - float lat, - float alt, - Planets planet, - Inputs input, - Report &report) { + float lat, + float alt, + Planets planet, + Inputs input, + Report &report) { std::string function = "get_bfield"; static int iFunction = -1; - report.enter(function, iFunction); + report.enter(function, iFunction); if (lat > pi/2) { lat = twopi - lat; @@ -30,7 +30,7 @@ bfield_info_type get_bfield(float lon, } bfield_info_type bfield_info; - + if (input.get_bfield_type() == "none") { bfield_info.b[0] = 0.0; bfield_info.b[1] = 0.0; @@ -40,9 +40,8 @@ bfield_info_type get_bfield(float lon, } else if (input.get_bfield_type() == "dipole") { bfield_info = get_dipole(lon, lat, alt, planet, input, report); } - + report.exit(function); return bfield_info; - } diff --git a/src/calc_chemical_sources.cpp b/src/calc_chemical_sources.cpp index ef4a010f..89f879e5 100644 --- a/src/calc_chemical_sources.cpp +++ b/src/calc_chemical_sources.cpp @@ -1,4 +1,4 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include @@ -6,25 +6,26 @@ #include "../include/report.h" #include "../include/chemistry.h" -void Chemistry::calc_chemical_sources(float *neutral_density, - float *ion_density, - float Tn, - float Ti, - float Te, - Report &report) { +void Chemistry::calc_chemical_sources(Neutrals &neutrals, + Ions &ions, + Report &report) { std::string function = "Chemistry::calc_chemical_sources"; static int iFunction = -1; - report.enter(function, iFunction); + report.enter(function, iFunction); - long iReaction, iLoss, iSource; + int64_t iReaction, iLoss, iSource; float change, rate; int IsNeutral, id_; - + + // This is make change the same size as the grid: + fcube change3d = neutrals.temperature_scgc; + // Then zero is out: + change3d.zeros(); + for (iReaction = 0; iReaction < nReactions; iReaction++) { - // std::cout << "iReaction : " << iReaction << " " << nReactions << "\n"; - // display_reaction(reactions[iReaction]); + if (report.test_verbose(8)) display_reaction(reactions[iReaction]); // First calculate reaction rate: @@ -32,52 +33,41 @@ void Chemistry::calc_chemical_sources(float *neutral_density, // Second calculate the amount of change: - change = rate; + change3d.fill(rate); for (iLoss = 0; iLoss < reactions[iReaction].nLosses; iLoss++) { - IsNeutral = reactions[iReaction].losses_IsNeutral[iLoss]; id_ = reactions[iReaction].losses_ids[iLoss]; - - if (IsNeutral) change = change * neutral_density[id_]; - else change = change * ion_density[id_]; - + if (IsNeutral) + change3d = change3d % neutrals.neutrals[id_].density_scgc; + else + change3d = change3d % ions.species[id_].density_scgc; } // Third add change to the different consituents: + // (a) First to losses: for (iLoss = 0; iLoss < reactions[iReaction].nLosses; iLoss++) { - IsNeutral = reactions[iReaction].losses_IsNeutral[iLoss]; id_ = reactions[iReaction].losses_ids[iLoss]; - if (IsNeutral) - sources_and_losses.neutral_losses[id_] = - sources_and_losses.neutral_losses[id_] + change; + neutrals.neutrals[id_].losses_scgc = + neutrals.neutrals[id_].losses_scgc + change3d; else - sources_and_losses.ion_losses[id_] = - sources_and_losses.ion_losses[id_] + change; - + ions.species[id_].losses_scgc = + ions.species[id_].losses_scgc + change3d; } + // (b) Second to sources: for (iSource = 0; iSource < reactions[iReaction].nSources; iSource++) { - IsNeutral = reactions[iReaction].sources_IsNeutral[iSource]; id_ = reactions[iReaction].sources_ids[iSource]; - if (IsNeutral) - sources_and_losses.neutral_sources[id_] = - sources_and_losses.neutral_sources[id_] + change; + neutrals.neutrals[id_].sources_scgc = + neutrals.neutrals[id_].sources_scgc + change3d; else - sources_and_losses.ion_sources[id_] = - sources_and_losses.ion_sources[id_] + change; - - } - - } - + ions.species[id_].sources_scgc = + ions.species[id_].sources_scgc + change3d; + } // for iSource + } // for iReaction report.exit(function); - - return; - } - diff --git a/src/calc_chemistry.cpp b/src/calc_chemistry.cpp index 20f83082..da020753 100644 --- a/src/calc_chemistry.cpp +++ b/src/calc_chemistry.cpp @@ -1,4 +1,4 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include @@ -14,110 +14,73 @@ #include "../include/solvers.h" void Chemistry::calc_chemistry(Neutrals &neutrals, - Ions &ions, - Times time, - Grid grid, - Report &report) { + Ions &ions, + Times time, + Grid grid, + Report &report) { - long nLons, nLats, nAlts, iLon, iLat, iAlt, index; int iSpecies; std::string function = "Chemistry::calc_chemistry"; static int iFunction = -1; - report.enter(function, iFunction); - - if (grid.get_IsGeoGrid()) { - nLons = nGeoLonsG; - nLats = nGeoLatsG; - nAlts = nGeoAltsG; - } else { - nLons = nMagLonsG; - nLats = nMagLatsG; - nAlts = nMagAltsG; - } - - float neutral_density[nSpecies]; - float ion_density[nIons+1]; // add one for electron density - float neutral_sources[nSpecies]; - float neutral_losses[nSpecies]; - float ion_sources[nIons]; - float ion_losses[nIons]; - - float Tn, Te, Ti; + report.enter(function, iFunction); float dt = time.get_dt(); - float old_density, source, loss, new_density; - float heat_neutrals, heat_ions, heat_electrons; - + // ------------------------------------ // Calculate electron densities // ------------------------------------ - ions.fill_electrons(grid, report); - - // Don't do chemistry in the ghostcells! - - for (iLon = 0; iLon < nLons; iLon++) { - for (iLat = 0; iLat < nLats; iLat++) { - for (iAlt = 0; iAlt < nAlts; iAlt++) { - - if (grid.get_IsGeoGrid()) { - index = ijk_geo_s3gc(iLon,iLat,iAlt); - } else { - index = ijk_mag_s3gc(iLon,iLat,iAlt); - } - - // Use the private variable sources_and_losses, so we don't - // have to pass it around - - for (iSpecies=0; iSpecies < nSpecies; iSpecies++) { - neutral_density[iSpecies] = neutrals.neutrals[iSpecies].density_s3gc[index]; - sources_and_losses.neutral_losses[iSpecies] = neutrals.neutrals[iSpecies].ionization_s3gc[index]; - sources_and_losses.neutral_sources[iSpecies] = 0.0; - } - - for (iSpecies=0; iSpecies < nIons; iSpecies++) { - ion_density[iSpecies] = ions.species[iSpecies].density_s3gc[index]; - sources_and_losses.ion_sources[iSpecies] = ions.species[iSpecies].ionization_s3gc[index]; - sources_and_losses.ion_losses[iSpecies] = 0.0; - } - ion_density[nIons] = ions.density_s3gc[index]; - - Tn = neutrals.temperature_s3gc[index]; - Ti = ions.ion_temperature_s3gc[index]; - Te = ions.electron_temperature_s3gc[index]; - - // If we wanted to do a higher-order solver, we would probably - // put it starting here: (If we do that, we may want to have a - // code that calculates the reaction rates first, then take - // this outside of the higher-order solver, since the reaction - // rates involve a lot of powers, which seem like there are - // quite slow in C.) - - calc_chemical_sources(neutral_density, - ion_density, - Tn, Ti, Te, report); - - for (iSpecies=0; iSpecies < nSpecies; iSpecies++) { - old_density = neutral_density[iSpecies]; - source = sources_and_losses.neutral_sources[iSpecies]; - loss = sources_and_losses.neutral_losses[iSpecies]; - neutrals.neutrals[iSpecies].density_s3gc[index] = - solver_chemistry(old_density, source, loss, dt); - } - - for (iSpecies=0; iSpecies < nIons; iSpecies++) { - old_density = ion_density[iSpecies]; - source = sources_and_losses.ion_sources[iSpecies]; - loss = sources_and_losses.ion_losses[iSpecies]; - ions.species[iSpecies].density_s3gc[index] = - solver_chemistry(old_density, source, loss, dt); - } - - } - } + ions.fill_electrons(report); + + // ---------------------------------------------------------- + // Initialize the sources and losses with EUV stuff: + // ---------------------------------------------------------- + + for (iSpecies=0; iSpecies < nSpecies; iSpecies++) { + neutrals.neutrals[iSpecies].losses_scgc = + neutrals.neutrals[iSpecies].ionization_scgc; + neutrals.neutrals[iSpecies].sources_scgc.zeros(); + } + + for (iSpecies=0; iSpecies < nIons; iSpecies++) { + ions.species[iSpecies].losses_scgc.zeros(); + ions.species[iSpecies].sources_scgc = + ions.species[iSpecies].ionization_scgc; + } + + // ---------------------------------------------------- + // Calculate the chemical sources and losses + // ---------------------------------------------------- + + calc_chemical_sources(neutrals, ions, report); + + // --------------------------------------------------------- + // Once sources and losses are done, solve for new densities + // --------------------------------------------------------- + + for (iSpecies=0; iSpecies < nSpecies; iSpecies++) { + neutrals.neutrals[iSpecies].density_scgc = + solver_chemistry(neutrals.neutrals[iSpecies].density_scgc, + neutrals.neutrals[iSpecies].sources_scgc, + neutrals.neutrals[iSpecies].losses_scgc, + dt); } - + + for (iSpecies=0; iSpecies < nIons; iSpecies++) { + ions.species[iSpecies].density_scgc = + solver_chemistry(ions.species[iSpecies].density_scgc, + ions.species[iSpecies].sources_scgc, + ions.species[iSpecies].losses_scgc, + dt); + } + + // --------------------------------------------------------- + // Recalculate electrons + // --------------------------------------------------------- + + ions.fill_electrons(report); + report.exit(function); return; } diff --git a/src/calc_euv.cpp b/src/calc_euv.cpp index e973e5dc..e8d728bb 100644 --- a/src/calc_euv.cpp +++ b/src/calc_euv.cpp @@ -1,4 +1,4 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include @@ -12,36 +12,34 @@ #include "../include/neutrals.h" #include "../include/ions.h" -int calc_euv( Planets planet, - Grid grid, - Times time, - Euv &euv, - Neutrals &neutrals, - Ions &ions, - Indices indices, - Inputs args, - Report &report) { - - int iErr=0; - - if (time.check_time_gate(args.get_dt_euv())) { +int calc_euv(Planets planet, + Grid grid, + Times time, + Euv &euv, + Neutrals &neutrals, + Ions &ions, + Indices indices, + Inputs args, + Report &report) { + + int iErr = 0; - std::string function="Euv::calc_euv"; + if (time.check_time_gate(args.get_dt_euv())) { + std::string function = "Euv::calc_euv"; static int iFunction = -1; - report.enter(function, iFunction); - + report.enter(function, iFunction); + // Chapman integrals for EUV energy deposition: neutrals.calc_chapman(grid, report); - + iErr = euv.euvac(time, indices, report); iErr = euv.scale_from_1au(planet, time); - + neutrals.calc_ionization_heating(euv, ions, report); report.exit(function); } return iErr; - } diff --git a/src/calc_neutral_derived.cpp b/src/calc_neutral_derived.cpp index 5480c938..adbe7ff6 100644 --- a/src/calc_neutral_derived.cpp +++ b/src/calc_neutral_derived.cpp @@ -1,4 +1,4 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include @@ -17,102 +17,70 @@ void Neutrals::calc_mass_density(Report &report) { - long iLon, iLat, iAlt, index, iSpecies; + int64_t iLon, iLat, iAlt, index, iSpecies; - std::string function="Neutrals::calc_mass_density"; + std::string function = "Neutrals::calc_mass_density"; static int iFunction = -1; - report.enter(function, iFunction); - - for (iLon = 0; iLon < nGeoLonsG; iLon++) { - for (iLat = 0; iLat < nGeoLatsG; iLat++) { - for (iAlt = 0; iAlt < nGeoAltsG; iAlt++) { - - index = ijk_geo_s3gc(iLon,iLat,iAlt); - - density_s3gc[index] = 0.0; - rho_s3gc[index] = 0.0; - - for (iSpecies=0; iSpecies < nSpecies; iSpecies++) { - density_s3gc[index] = density_s3gc[index] + - neutrals[iSpecies].density_s3gc[index]; - rho_s3gc[index] = rho_s3gc[index] + - neutrals[iSpecies].mass * - neutrals[iSpecies].density_s3gc[index]; - } - - mean_major_mass_s3gc[index] = rho_s3gc[index] / density_s3gc[index]; - pressure_s3gc[index] = - density_s3gc[index] * - boltzmanns_constant * - temperature_s3gc[index]; - } - } + report.enter(function, iFunction); + + rho_scgc.zeros(); + density_scgc.zeros(); + + for (iSpecies=0; iSpecies < nSpecies; iSpecies++) { + rho_scgc = rho_scgc + + neutrals[iSpecies].mass * neutrals[iSpecies].density_scgc; + density_scgc = density_scgc + neutrals[iSpecies].density_scgc; } - + + mean_major_mass_scgc = rho_scgc / density_scgc; + pressure_scgc = boltzmanns_constant * density_scgc % temperature_scgc; + report.exit(function); return; - } - + //---------------------------------------------------------------------- // //---------------------------------------------------------------------- void Neutrals::calc_specific_heat(Report &report) { - long iLon, iLat, iAlt, index, iSpecies; + int64_t iLon, iLat, iAlt, index, iSpecies; double t, p, r; - - std::string function="Neutrals::calc_specific_heat"; - static int iFunction = -1; - report.enter(function, iFunction); - - for (iLon = 0; iLon < nGeoLonsG; iLon++) { - for (iLat = 0; iLat < nGeoLatsG; iLat++) { - for (iAlt = 0; iAlt < nGeoAltsG; iAlt++) { - - index = ijk_geo_s3gc(iLon,iLat,iAlt); - - Cv_s3gc[index] = 0.0; - gamma_s3gc[index] = 0.0; - kappa_s3gc[index] = 0.0; - - t = temperature_s3gc[index]; - - for (iSpecies=0; iSpecies < nSpecies; iSpecies++) { - p = neutrals[iSpecies].thermal_exp; - r = pow(t,p); // temp ^ exp; <--- This line slows the code way down. - Cv_s3gc[index] = Cv_s3gc[index] + - (neutrals[iSpecies].vibe - 2) * - neutrals[iSpecies].density_s3gc[index] * - boltzmanns_constant / neutrals[iSpecies].mass; - gamma_s3gc[index] = gamma_s3gc[index] + - neutrals[iSpecies].density_s3gc[index] / (neutrals[iSpecies].vibe-2); - kappa_s3gc[index] = kappa_s3gc[index] + - neutrals[iSpecies].density_s3gc[index] * - neutrals[iSpecies].thermal_cond * - r; - } - - Cv_s3gc[index] = Cv_s3gc[index] / (2*density_s3gc[index]); - gamma_s3gc[index] = gamma_s3gc[index] * 2.0 / density_s3gc[index] + 1.0; - kappa_s3gc[index] = kappa_s3gc[index] / density_s3gc[index]; - - r = gamma_s3gc[index] * - boltzmanns_constant * - temperature_s3gc[index] / - mean_major_mass_s3gc[index]; - - sound_s3gc[index] = sqrt(r); // <---- Same with this line - } - } + std::string function = "Neutrals::calc_specific_heat"; + static int iFunction = -1; + report.enter(function, iFunction); + + Cv_scgc.zeros(); + gamma_scgc.zeros(); + kappa_scgc.zeros(); + + for (iSpecies=0; iSpecies < nSpecies; iSpecies++) { + Cv_scgc = Cv_scgc + + (neutrals[iSpecies].vibe - 2) * + neutrals[iSpecies].density_scgc * + boltzmanns_constant / neutrals[iSpecies].mass; + gamma_scgc = gamma_scgc + + neutrals[iSpecies].density_scgc / (neutrals[iSpecies].vibe-2); + kappa_scgc = kappa_scgc + + neutrals[iSpecies].thermal_cond * + neutrals[iSpecies].density_scgc % + pow(temperature_scgc, neutrals[iSpecies].thermal_exp); } + Cv_scgc = Cv_scgc / (2*density_scgc); + gamma_scgc = gamma_scgc * 2.0 / density_scgc + 1.0; + kappa_scgc = kappa_scgc / density_scgc; + + sound_scgc = sqrt(boltzmanns_constant * + gamma_scgc % + temperature_scgc / + mean_major_mass_scgc); + report.exit(function); return; - } //---------------------------------------------------------------------- @@ -124,18 +92,18 @@ void Neutrals::calc_specific_heat(Report &report) { void Neutrals::calc_chapman(Grid grid, Report &report) { - long iAlt, iLon, iLat, index, indexp; + int64_t iAlt, iLon, iLat, index, indexp; float H; // scale height // This is all from Smith and Smith, JGR 1972, vol. 77, page 3592 // "Numerical evaluation of chapman's grazing incidence integral ch(X,x)" // Xp is supposed to be R/H - // JMB Update: 05/2017. Corrected a small error in the y-calc for + // JMB Update: 05/2017. Corrected a small error in the y-calc for // erfc(y) // // Also Updated the Grazing Integral for SZA > 90.0 // We now do log-linear interpolation for smoother transitions - + double a = 1.06069630; double b = 0.55643831; double c = 1.06198960; @@ -143,152 +111,140 @@ void Neutrals::calc_chapman(Grid grid, Report &report) { double f = 0.56498823; double g = 0.06651874; - double integral[nGeoAltsG], xp[nGeoAltsG], erfcy[nGeoAltsG]; - double log_int[nGeoAltsG]; double y, dy; float Hp_up, Hp_dn, grad_hs, grad_xp, grad_in, Hg, Xg, in, int_g, int_p; - long index_bottom, iindex, iindexp, iiAlt; - - std::string function="Neutrals::calc_chapman"; + int64_t index_bottom, iindex, iindexp, iiAlt; + + std::string function = "Neutrals::calc_chapman"; static int iFunction = -1; - report.enter(function, iFunction); + report.enter(function, iFunction); + + int64_t nLons = grid.get_nLons(); + int64_t nLats = grid.get_nLats(); + int64_t nAlts = grid.get_nAlts(); + int64_t nGCs = grid.get_nGCs(); + + // New way of doing it with 3D arrays: + + fcube integral3d(nLons, nLats, nAlts); + fcube log_int3d(nLons, nLats, nAlts); + fcube xp3d(nLons, nLats, nAlts); + fcube y3d(nLons, nLats, nAlts); + fcube erfcy3d(nLons, nLats, nAlts); + + fvec integral1d(nAlts); + fvec log_int1d(nAlts); + fvec xp1d(nAlts); + fvec y1d(nAlts); + fvec erfcy1d(nAlts); + fvec dAlt1d(nAlts); + fvec sza1d(nAlts); + fvec radius1d(nAlts); + fvec H1d(nAlts); for (int iSpecies=0; iSpecies < nSpecies; iSpecies++) { - for (iLon = 0; iLon < nGeoLonsG; iLon++) { - for (iLat = 0; iLat < nGeoLatsG; iLat++) { - - iAlt = nGeoAltsG-1; - - index = ijk_geo_s3gc(iLon,iLat,iAlt); - - // The integral from the top to infinity is the density * H (scale height) - // So calculate the scale height: - - H = calc_scale_height(iSpecies, index, grid); - - integral[iAlt] = neutrals[iSpecies].density_s3gc[index] * H; - log_int[iAlt] = log(integral[iAlt]); - - // if we wanted to do this properly, then the integral should be - // with cell edge values and the distances from cell centers to - // cell centers. But, we are approximating here: - - for (int iAlt=nGeoAltsG-1; iAlt>=0; iAlt--) { - - index = ijk_geo_s3gc(iLon,iLat,iAlt); - indexp = ijk_geo_s3gc(iLon,iLat,iAlt+1); - - if (iAlt < nGeoAltsG-1) { - integral[iAlt] = integral[iAlt+1] + - neutrals[iSpecies].density_s3gc[index] * grid.dalt_lower_s3gc[indexp]; - log_int[iAlt] = log(integral[iAlt]); - } - - H = calc_scale_height(iSpecies, index, grid); - - xp[iAlt] = grid.radius_s3gc[index] / H; - - // Eqn (10) Smith & Smith - y = sqrt(0.5 * xp[iAlt]) * fabs(grid.cos_sza_s3gc[index]); - - // Eqn (12) Smith and Smith - if (y < 8) erfcy[iAlt] = (a + b*y) / (c + d*y + y*y); - else erfcy[iAlt] = f / (g + y); - - } - - // Don't need chapman integrals in the lower ghostcells: - - for (int iAlt=0; iAlt < nGeoGhosts; iAlt++) { - index = ijk_geo_s3gc(iLon,iLat,iAlt); - neutrals[iSpecies].chapman_s3gc[index] = max_chapman; - } - - // Rest of domain: - - for (int iAlt=nGeoGhosts; iAlt < nGeoAltsG; iAlt++) { - - index = ijk_geo_s3gc(iLon,iLat,iAlt); - - // This is on the dayside: - if (grid.sza_s3gc[index] < pi/2 || grid.sza_s3gc[index] > 3*pi/2) { - - neutrals[iSpecies].chapman_s3gc[index] = - integral[iAlt] * sqrt(0.5 * pi * xp[iAlt]) * erfcy[iAlt]; - - } else { - - // This is on the nghtside of the terminator: - - y = grid.radius_s3gc[index] * abs(cos(grid.sza_s3gc[index]-pi/2)); - - // This sort of assumes that nGeoGhosts >= 2: - index_bottom = ijk_geo_s3gc(iLon,iLat,nGeoGhosts); - if (y > grid.radius_s3gc[index_bottom]) { - - iiAlt = iAlt; - iindex = ijk_geo_s3gc(iLon,iLat,iiAlt-1); - while (grid.radius_s3gc[iindex]>y) { - iiAlt--; - iindex = ijk_geo_s3gc(iLon,iLat,iiAlt-1); - } - iiAlt--; - - iindexp = ijk_geo_s3gc(iLon,iLat,iiAlt+1); - iindex = ijk_geo_s3gc(iLon,iLat,iiAlt); - - Hp_up = calc_scale_height(iSpecies, iindexp, grid); - Hp_dn = calc_scale_height(iSpecies, iindex, grid); - - // make sure to use the proper cell spacing (iiAlt+1 & lower): - grad_hs = (Hp_up - Hp_dn) / grid.dalt_lower_s3gc[iindexp]; - grad_xp = (xp[iiAlt+1]-xp[iiAlt]) / grid.dalt_lower_s3gc[iindexp]; - grad_in = (log_int[iiAlt+1] - log_int[iiAlt]) / grid.dalt_lower_s3gc[iindexp]; - - // Linearly interpolate H and X: - dy = y - grid.radius_s3gc[iindex]; - Hg = Hp_dn + grad_hs * dy; - Xg = xp[iiAlt] + grad_xp * dy; - in = log_int[iiAlt] + grad_in * dy; - - int_g = exp(in); - int_p = integral[iAlt]; - // Equation (19) Smith & Smith - neutrals[iSpecies].chapman_s3gc[index] = - sqrt(0.5 * pi * Xg) * (2.0 * int_g - int_p * erfcy[iAlt]); + neutrals[iSpecies].scale_height_scgc = + boltzmanns_constant * temperature_scgc / + (neutrals[iSpecies].mass * grid.gravity_scgc); - if (neutrals[iSpecies].chapman_s3gc[index] > max_chapman) - neutrals[iSpecies].chapman_s3gc[index] = max_chapman; - - } else { + xp3d = grid.radius_scgc / neutrals[iSpecies].scale_height_scgc; + y3d = sqrt(0.5 * xp3d) % abs(grid.cos_sza_scgc); + iAlt = nAlts-1; - // This says that we are in the shadow of the planet: + integral3d.slice(iAlt) = + neutrals[iSpecies].density_scgc.slice(iAlt) % + neutrals[iSpecies].scale_height_scgc.slice(iAlt); - neutrals[iSpecies].chapman_s3gc[index] = max_chapman; - - } - - } - - if (report.test_verbose(10)) - std::cout << "iSpecies, iAlt, chap : " << iSpecies << " " << iAlt << " " << - grid.sza_s3gc[index]*rtod << " " << - xp[iAlt] << " " << - erfcy[iAlt] << " " << - neutrals[iSpecies].chapman_s3gc[index] << " " << integral[iAlt] << "\n"; + for (iAlt = nAlts-1; iAlt >= 0; iAlt--) { + if (iAlt < nAlts-1) { + integral3d.slice(iAlt) = integral3d.slice(iAlt+1) + + neutrals[iSpecies].density_scgc.slice(iAlt) % + grid.dalt_lower_scgc.slice(iAlt+1); + } + } - } // iAlt + erfcy3d = (a + b * y3d) / (c + d*y3d + y3d % y3d); + for (iLon = 0; iLon < nLons ; iLon++) + for (iLat = 0; iLat < nLats ; iLat++) + for (iAlt = 0; iAlt < nAlts ; iAlt++) + if (y3d(iLon, iLat, iAlt) >= 8.0) + erfcy3d(iLon, iLat, iAlt) = f / (g + y3d(iLon, iLat, iAlt)); + + log_int3d = log(integral3d); + + // Set chapman integrals to max in the lower ghostcells + + for (int iAlt=0; iAlt < nGCs; iAlt++) + neutrals[iSpecies].chapman_scgc.slice(iAlt).fill(max_chapman); + + for (iLon = 0; iLon < nLons ; iLon++) { + for (iLat = 0; iLat < nLats ; iLat++) { + + dAlt1d = grid.dalt_lower_scgc.tube(iLon, iLat); + sza1d = grid.sza_scgc.tube(iLon, iLat); + integral1d = integral3d.tube(iLon, iLat); + log_int1d = log_int3d.tube(iLon, iLat); + xp1d = xp3d.tube(iLon, iLat); + y1d = y3d.tube(iLon, iLat); + erfcy1d = erfcy3d.tube(iLon, iLat); + radius1d = grid.radius_scgc.tube(iLon, iLat); + H1d = neutrals[iSpecies].scale_height_scgc.tube(iLon, iLat); + + for (iAlt = nGCs; iAlt < nAlts; iAlt++) { + // This is on the dayside: + if (sza1d(iAlt) < pi/2 || sza1d(iAlt) > 3*pi/2) { + neutrals[iSpecies].chapman_scgc(iLon, iLat, iAlt) = + integral1d(iAlt) * sqrt(0.5 * pi * xp1d(iAlt)) * erfcy1d(iAlt); + } else { + // This is on the nghtside of the terminator: + + y = radius1d(iAlt) * abs(cos(sza1d(iAlt)-pi/2)); + + // This sort of assumes that nGeoGhosts >= 2: + if (y > radius1d(nGCs)) { + + iiAlt = iAlt; + while (radius1d(iiAlt-1) > y) iiAlt--; + iiAlt--; + + Hp_up = H1d(iiAlt+1); + Hp_dn = H1d(iiAlt); + + // make sure to use the proper cell spacing (iiAlt+1 & lower): + grad_hs = (Hp_up - Hp_dn) / dAlt1d(iiAlt+1); + grad_xp = (xp1d(iiAlt+1) - xp1d(iiAlt)) / dAlt1d(iiAlt+1); + grad_in = (log_int1d(iiAlt+1) - log_int1d(iiAlt)) / dAlt1d(iiAlt+1); + + // Linearly interpolate H and X: + dy = y - radius1d(iiAlt); + Hg = Hp_dn + grad_hs * dy; + Xg = xp1d(iiAlt) + grad_xp * dy; + in = log_int1d(iiAlt) + grad_in * dy; + + int_g = exp(in); + int_p = integral1d(iAlt); + // Equation (19) Smith & Smith + neutrals[iSpecies].chapman_scgc(iLon, iLat, iAlt) = + sqrt(0.5 * pi * Xg) * (2.0 * int_g - int_p * erfcy1d(iAlt)); + + if (neutrals[iSpecies].chapman_scgc(iLon, iLat, iAlt) > max_chapman) + neutrals[iSpecies].chapman_scgc(iLon, iLat, iAlt) = max_chapman; + + } else { + // This says that we are in the shadow of the planet: + + neutrals[iSpecies].chapman_scgc(iLon, iLat, iAlt) = max_chapman; + } + } + } // iAlt + } // iLat + } // iLon + } // iSpecies - } // iLat - } // iLon - - } // iSpecies - report.exit(function); return; - } // ----------------------------------------------------------------------------- @@ -297,102 +253,56 @@ void Neutrals::calc_chapman(Grid grid, Report &report) { void Neutrals::calc_conduction(Grid grid, Times time, Report &report) { - float prandtl[nGeoAltsG]; - float rhocv[nGeoAltsG]; - float lambda[nGeoAltsG]; - float conduction[nGeoAltsG]; - float temp[nGeoAltsG]; - float dalt_lower[nGeoAltsG]; float dt; - - long iLon, iLat, iAlt, index; - std::string function="Neutrals::calc_conduction"; + int64_t iLon, iLat, iAlt, index; + + std::string function = "Neutrals::calc_conduction"; static int iFunction = -1; - report.enter(function, iFunction); - - for (iLon = 0; iLon < nGeoLonsG; iLon++) { - for (iLat = 0; iLat < nGeoLatsG; iLat++) { - for (iAlt=0; iAlt < nGeoAltsG; iAlt++) { - index = ijk_geo_s3gc(iLon,iLat,iAlt); - conduction_s3gc[index] = 0.0; - } - } - } - - for (iLon = iGeoLonStart_; iLon < iGeoLonEnd_; iLon++) { - for (iLat = iGeoLatStart_; iLat < iGeoLatEnd_; iLat++) { - - // Treat each altitude slice individually: - - for (iAlt=0; iAlt < nGeoAltsG; iAlt++) { - index = ijk_geo_s3gc(iLon,iLat,iAlt); - - rhocv[iAlt] = rho_s3gc[index] * Cv_s3gc[index]; - // rhocv needs to be scaled by radius squared: - rhocv[iAlt] = rhocv[iAlt] * grid.radius_sq_s3gc[index]; - - // Need to make this eddy * rho * cv: - prandtl[iAlt] = 0.0; - - lambda[iAlt] = kappa_s3gc[index] + prandtl[iAlt]; - // lambda needs to be scaled by radius squared: - lambda[iAlt] = lambda[iAlt] * grid.radius_sq_s3gc[index]; - - conduction[iAlt] = 0.0; - - temp[iAlt] = temperature_s3gc[index]; - - dalt_lower[iAlt] = grid.dalt_lower_s3gc[index]; - - } + report.enter(function, iFunction); - dt = time.get_dt(); + int64_t nLons = grid.get_nLons(); + int64_t nLats = grid.get_nLats(); + int64_t nAlts = grid.get_nAlts(); - // if (iLon == nGeoLonsG/2 && iLat == nGeoLatsG/2) { - // for (iAlt=0; iAlt < nGeoAltsG; iAlt++) { - // std::cout << "before conduction : " - // << iLon << " " << iLat << " " << iAlt << " " - // << temp[iAlt] << " " - // << lambda[iAlt] << " " - // << rhocv[iAlt] << " " - // << dt << " " - // << dalt_lower[iAlt] << "\n"; - // } - // } - - solver_conduction(temp, lambda, rhocv, dt, dalt_lower, conduction); - - // if (iLon == nGeoLonsG/2 && iLat == nGeoLatsG/2) { - // for (iAlt=0; iAlt < nGeoAltsG; iAlt++) { - // std::cout << "after conduction : " - // << iLon << " " << iLat << " " << iAlt << " " - // << conduction[iAlt] << "\n"; - // } - // } - // We want the sources to be in terms of dT/dt, while the - // conduction actually solves for Tnew-Told, so divide by dt + fcube rhocvr23d(nLons, nLats, nAlts); + fcube lambda3d(nLons, nLats, nAlts); + fcube prandtl3d(nLons, nLats, nAlts); - for (iAlt=0; iAlt < nGeoAltsG; iAlt++) { - conduction_s3gc[index] = conduction[iAlt]/dt; + rhocvr23d = rho_scgc % Cv_scgc % grid.radius2_scgc; + // Need to make this eddy * rho * cv: + prandtl3d.zeros(); + lambda3d = (kappa_scgc + prandtl3d) % grid.radius2_scgc; - if (report.test_verbose(10)) - std::cout << "conduction : " << index - << " " << conduction_s3gc[index]*seconds_per_day << " deg/day\n"; - - } + fvec temp1d(nAlts); + fvec lambda1d(nAlts); + fvec rhocvr21d(nAlts); + fvec dalt1d(nAlts); + fvec conduction1d(nAlts); - // End of calculation for each altitude + for (iLon = 0; iLon < nLons; iLon++) { + for (iLat = 0; iLat < nLats; iLat++) { - } // lat + temp1d = temperature_scgc.tube(iLon, iLat); + lambda1d = lambda3d.tube(iLon, iLat); + rhocvr21d = rhocvr23d.tube(iLon, iLat); + dalt1d = grid.dalt_lower_scgc.tube(iLon, iLat); + conduction1d.zeros(); - } // lon + dt = time.get_dt(); - report.exit(function); + conduction1d = solver_conduction(temp1d, lambda1d, rhocvr21d, dt, dalt1d); + + // We want the sources to be in terms of dT/dt, while the + // conduction actually solves for Tnew-Told, so divide by dt + conduction_scgc.tube(iLon, iLat) = conduction1d / dt; + } // lat + } // lon + report.exit(function); } - + // ----------------------------------------------------------------------------- // Calculate EUV driven ionization and heating rates @@ -400,146 +310,82 @@ void Neutrals::calc_conduction(Grid grid, Times time, Report &report) { void Neutrals::calc_ionization_heating(Euv euv, Ions &ions, Report &report) { - long iAlt, iLon, iLat, iWave, iSpecies, index, indexp; + int64_t iAlt, iLon, iLat, iWave, iSpecies, index, indexp; int i_, idion_, ideuv_, nIonizations, iIon, iIonization; float tau, intensity, photoion; float ionization; - std::string function="calc_ionization_heating"; + std::string function = "calc_ionization_heating"; static int iFunction = -1; - report.enter(function, iFunction); + report.enter(function, iFunction); // Zero out all source terms: - - for (iLon = 0; iLon < nGeoLonsG; iLon++) { - for (iLat = 0; iLat < nGeoLatsG; iLat++) { - for (iAlt = 0; iAlt < nGeoAltsG; iAlt++) { - index = ijk_geo_s3gc(iLon,iLat,iAlt); + heating_euv_scgc.zeros(); + for (iSpecies=0; iSpecies < nSpecies; iSpecies++) + neutrals[iSpecies].ionization_scgc.zeros(); - heating_euv_s3gc[index] = 0.0; - for (iSpecies=0; iSpecies < nSpecies; iSpecies++) - neutrals[iSpecies].ionization_s3gc[index] = 0.0; + int64_t nLons = heating_euv_scgc.n_rows; + int64_t nLats = heating_euv_scgc.n_cols; + int64_t nAlts = heating_euv_scgc.n_slices; + fmat tau2d = heating_euv_scgc.slice(0); + fmat intensity2d = heating_euv_scgc.slice(0); + fmat ionization2d = heating_euv_scgc.slice(0); + + for (iAlt = 2; iAlt < nAlts-2; iAlt++) { + for (iWave=0; iWave < euv.nWavelengths; iWave++) { + + tau2d.zeros(); + + for (iSpecies=0; iSpecies < nSpecies; iSpecies++) { + if (neutrals[iSpecies].iEuvAbsId_ > -1) { + i_ = neutrals[iSpecies].iEuvAbsId_; + tau2d = tau2d + + euv.waveinfo[i_].values[iWave] * + neutrals[iSpecies].chapman_scgc.slice(iAlt); + } } - } - } - - for (iLon = iGeoLonStart_; iLon < iGeoLonEnd_; iLon++) { - for (iLat = iGeoLatStart_; iLat < iGeoLatEnd_; iLat++) { - for (iAlt = iGeoAltStart_; iAlt < iGeoAltEnd_; iAlt++) { - - index = ijk_geo_s3gc(iLon,iLat,iAlt); - - heating_euv_s3gc[index] = 0.0; - for (iSpecies=0; iSpecies < nSpecies; iSpecies++) - neutrals[iSpecies].ionization_s3gc[index] = 0.0; - - for (iWave=0; iWave < euv.nWavelengths; iWave++) { - - // Need to calculate the intensity at particular wavelength - // (iWave), for each species (iSpecies), which is dependent - // on the chapman integrals at that particular location (index) - // and the cross section (i_): - - tau = 0.0; - for (iSpecies=0; iSpecies < nSpecies; iSpecies++) { - if (neutrals[iSpecies].iEuvAbsId_ > -1) { - i_ = neutrals[iSpecies].iEuvAbsId_; - tau = tau + - euv.waveinfo[i_].values[iWave] * - neutrals[iSpecies].chapman_s3gc[index]; - } - } - - intensity = euv.wavelengths_intensity_top[iWave] * exp(-1.0*tau); - - for (iSpecies=0; iSpecies < nSpecies; iSpecies++) { - - // Calculate Photo-Absorbtion for each species and add them up: - i_ = neutrals[iSpecies].iEuvAbsId_; // index of photo abs cross section - if (i_ > -1) { - heating_euv_s3gc[index] = heating_euv_s3gc[index] + - intensity * - euv.wavelengths_energy[iWave] * - euv.waveinfo[i_].values[iWave] * // cross section - neutrals[iSpecies].density_s3gc[index]; - } - - for (iIonization = 0; iIonization < neutrals[iSpecies].nEuvIonSpecies; iIonization++) { - - i_ = neutrals[iSpecies].iEuvIonId_[iIonization]; - // std::cout << iSpecies << " " << iIonization << " " << i_ << "\n"; - - ionization = - intensity * - euv.waveinfo[i_].values[iWave] * // cross section - neutrals[iSpecies].density_s3gc[index]; - - neutrals[iSpecies].ionization_s3gc[index] = - neutrals[iSpecies].ionization_s3gc[index] + ionization; - - iIon = neutrals[iSpecies].iEuvIonSpecies_[iIonization]; - ions.species[iIon].ionization_s3gc[index] = - ions.species[iIon].ionization_s3gc[index] + ionization; - - } - - } // Each species - - - } // Each wavelength - - // Scale heating with efficiency, and - // convert energy deposition to change in temperature: - - heating_euv_s3gc[index] = - heating_efficiency * - heating_euv_s3gc[index] / rho_s3gc[index] / Cv_s3gc[index]; - - if (report.test_verbose(10)) - std::cout << "heating : " << index - << " " << heating_euv_s3gc[index]*seconds_per_day << " deg/day\n"; - - } // Alts - } // Lats - } // Lons - -// // We need to do things here: -// // - Identify where the ionization cross section is stored -// // - Identify which ion gets the ionization -// // - Keep track of neutral losses -// -// // How many ionizations do we have: -// nIonizations = neutrals[iSpecies].iEuvIonId_.size(); -// for (int iIon=0; iIon < nIonizations; iIon++) { -// -// // Find cross section: -// ideuv_ = neutrals[iSpecies].iEuvIonId_[iIon]; -// // Find ion to put sources: -// idion_ = neutrals[iSpecies].iEuvIonSpecies_[iIon]; -// -// photoion = euv.waveinfo[ideuv_].values[iWave]; -// -// // This is the ionization rate: -// ionization = intensity * photoion * -// neutrals[iSpecies].density[iAlt]; -// -// // Keep track of losses for neutrals: -// neutrals[iSpecies].ionization[iAlt] = -// neutrals[iSpecies].ionization[iAlt] + -// ionization; -// -// // Keep track of sources for ions: -// ions.species[idion_].ionization[iAlt] = -// ions.species[idion_].ionization[iAlt] + -// ionization; -// -// } -// + + intensity2d = euv.wavelengths_intensity_top[iWave] * exp(-1.0*tau2d); + + for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) { + // Calculate Photo-Absorbtion for each species and add them up: + i_ = neutrals[iSpecies].iEuvAbsId_; // index of photo abs cross section + if (i_ > -1) { + heating_euv_scgc.slice(iAlt) = heating_euv_scgc.slice(iAlt) + + euv.wavelengths_energy[iWave] * + euv.waveinfo[i_].values[iWave] * // cross section + (intensity2d % + neutrals[iSpecies].density_scgc.slice(iAlt) ); + } + + for (iIonization = 0; + iIonization < neutrals[iSpecies].nEuvIonSpecies; + iIonization++) { + + i_ = neutrals[iSpecies].iEuvIonId_[iIonization]; + + ionization2d = + euv.waveinfo[i_].values[iWave] * // cross section + intensity2d % + neutrals[iSpecies].density_scgc.slice(iAlt); + + neutrals[iSpecies].ionization_scgc.slice(iAlt) = + neutrals[iSpecies].ionization_scgc(iAlt) + ionization2d; + + iIon = neutrals[iSpecies].iEuvIonSpecies_[iIonization]; + ions.species[iIon].ionization_scgc.slice(iAlt) = + ions.species[iIon].ionization_scgc.slice(iAlt) + ionization2d; + } // iIonization + } // iSpecies + } // iWave + } // iAlt + + heating_euv_scgc = + heating_efficiency * heating_euv_scgc / rho_scgc / Cv_scgc; report.exit(function); return; - } diff --git a/src/chemistry.cpp b/src/chemistry.cpp index ddfbc471..de25e644 100644 --- a/src/chemistry.cpp +++ b/src/chemistry.cpp @@ -1,4 +1,4 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include @@ -16,21 +16,20 @@ // ----------------------------------------------------------------------------- Chemistry::Chemistry(Neutrals neutrals, - Ions ions, - Inputs args, - Report &report) { + Ions ions, + Inputs args, + Report &report) { std::string function = "Chemistry::Chemistry"; static int iFunction = -1; - report.enter(function, iFunction); + report.enter(function, iFunction); int iErr = 0; read_chemistry_file(neutrals, ions, args, report); - + report.exit(function); return; - } // ----------------------------------------------------------------------------- @@ -38,13 +37,13 @@ Chemistry::Chemistry(Neutrals neutrals, // ----------------------------------------------------------------------------- int Chemistry::read_chemistry_file(Neutrals neutrals, - Ions ions, - Inputs args, - Report &report) { + Ions ions, + Inputs args, + Report &report) { std::string function = "Chemistry::read_chemistry_file"; static int iFunction = -1; - report.enter(function, iFunction); + report.enter(function, iFunction); std::ifstream infile_ptr; int iErr = 0; @@ -66,35 +65,28 @@ int Chemistry::read_chemistry_file(Neutrals neutrals, int nLines = csv.size(); if (nLines <= 2) { - iErr = 1; + iErr = 1; } else { - - // Skip 2 lines of headers! - - nReactions = 0; - - for (int iLine = 2; iLine < nLines; iLine++) { - - // Some final rows can have comments in them, so we want to - // skip anything where the length of the string in column 2 - // is == 0: - - if (csv[iLine][1].length() > 0) { - reaction = interpret_reaction_line(neutrals, ions, csv[iLine], report); - if (reaction.nLosses > 0 && reaction.nSources > 0) { - reactions.push_back(reaction); - nReactions++; - } - } - - } - + // Skip 2 lines of headers! + + nReactions = 0; + + for (int iLine = 2; iLine < nLines; iLine++) { + // Some final rows can have comments in them, so we want to + // skip anything where the length of the string in column 2 + // is == 0: + if (csv[iLine][1].length() > 0) { + reaction = interpret_reaction_line(neutrals, ions, + csv[iLine], report); + if (reaction.nLosses > 0 && reaction.nSources > 0) { + reactions.push_back(reaction); + nReactions++; + } + } + } } - } - } - report.exit(function); return iErr; } @@ -104,22 +96,22 @@ int Chemistry::read_chemistry_file(Neutrals neutrals, // ----------------------------------------------------------------------------- Chemistry::reaction_type Chemistry::interpret_reaction_line(Neutrals neutrals, - Ions ions, - std::vector line, - Report &report) { + Ions ions, + std::vector line, + Report &report) { std::string function = "Chemistry::interpret_reaction_line"; static int iFunction = -1; - report.enter(function, iFunction); + report.enter(function, iFunction); reaction_type reaction; int i; int id_, IsNeutral; - + // Losses (left side) first: reaction.nLosses = 0; - for (i=0; i<3; i++) { + for (i = 0; i < 3; i++) { find_species_id(line[i], neutrals, ions, id_, IsNeutral, report); if (id_ >= 0) { reaction.losses_names.push_back(line[i]); @@ -131,7 +123,7 @@ Chemistry::reaction_type Chemistry::interpret_reaction_line(Neutrals neutrals, // Sources (right side) second: reaction.nSources = 0; - for (i=4; i<7; i++) { + for (i = 4; i < 7; i++) { find_species_id(line[i], neutrals, ions, id_, IsNeutral, report); if (id_ >= 0) { reaction.sources_names.push_back(line[i]); @@ -149,7 +141,7 @@ Chemistry::reaction_type Chemistry::interpret_reaction_line(Neutrals neutrals, // energy released as exo-thermic reaction: reaction.energy = stof(line[9]); - + report.exit(function); return reaction; } @@ -159,40 +151,39 @@ Chemistry::reaction_type Chemistry::interpret_reaction_line(Neutrals neutrals, // ----------------------------------------------------------------------------- void Chemistry::find_species_id(std::string name, - Neutrals neutrals, - Ions ions, - int &id_, - int &IsNeutral, - Report &report) { + Neutrals neutrals, + Ions ions, + int &id_, + int &IsNeutral, + Report &report) { std::string function = "Chemistry::find_species_id"; static int iFunction = -1; - report.enter(function, iFunction); + report.enter(function, iFunction); int iSpecies; - + id_ = -1; IsNeutral = -1; if (name.length() > 0) { // Check Neutrals: - for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) + for (iSpecies = 0; iSpecies < nSpecies; iSpecies++) if (name == neutrals.neutrals[iSpecies].cName) { - id_ = iSpecies; - IsNeutral = 1; - break; + id_ = iSpecies; + IsNeutral = 1; + break; } - + if (id_ == -1) { // Check Ions: - for (iSpecies = 0; iSpecies <= nIons; iSpecies++) - if (name == ions.species[iSpecies].cName) { - id_ = iSpecies; - IsNeutral = 0; - break; - } + for (iSpecies = 0; iSpecies <= nIons; iSpecies++) + if (name == ions.species[iSpecies].cName) { + id_ = iSpecies; + IsNeutral = 0; + break; + } } - } report.exit(function); @@ -216,12 +207,14 @@ void Chemistry::display_reaction(Chemistry::reaction_type reaction) { for (i = 0; i < reaction.nSources; i++) std::cout << reaction.sources_names[i] << " + "; std::cout << " ( RR : " << reaction.rate << ")\n"; - + for (i = 0; i < reaction.nLosses; i++) - std::cout << reaction.losses_ids[i] << "(" << reaction.losses_IsNeutral[i] << ")" << " + "; + std::cout << reaction.losses_ids[i] + << "(" << reaction.losses_IsNeutral[i] << ")" << " + "; std::cout << " -> "; for (i = 0; i < reaction.nSources; i++) - std::cout << reaction.sources_ids[i] << "(" << reaction.sources_IsNeutral[i] << ")" << " + "; + std::cout << reaction.sources_ids[i] + << "(" << reaction.sources_IsNeutral[i] + << ")" << " + "; std::cout << " ( RR : " << reaction.rate << ")\n"; - } diff --git a/src/dipole.cpp b/src/dipole.cpp index d9de113d..752fbcef 100644 --- a/src/dipole.cpp +++ b/src/dipole.cpp @@ -1,4 +1,4 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include @@ -12,22 +12,22 @@ #include "../include/constants.h" bfield_info_type get_dipole(float lon, - float lat, - float alt, - Planets planet, - Inputs input, - Report &report) { + float lat, + float alt, + Planets planet, + Inputs input, + Report &report) { std::string function = "dipole"; static int iFunction = -1; - report.enter(function, iFunction); + report.enter(function, iFunction); bfield_info_type bfield_info; - + float llr[3]; float xyz[3]; float radius = planet.get_radius(lat); - + llr[0] = lon; llr[1] = lat; llr[2] = alt + radius; @@ -53,14 +53,14 @@ bfield_info_type get_dipole(float lon, float yzpp2 = (pos_rot_zy[1]*pos_rot_zy[1] + pos_rot_zy[2]*pos_rot_zy[2]); float xyzpp2 = (pos_rot_zy[0]*pos_rot_zy[0] + - pos_rot_zy[1]*pos_rot_zy[1] + - pos_rot_zy[2]*pos_rot_zy[2]); + pos_rot_zy[1]*pos_rot_zy[1] + + pos_rot_zy[2]*pos_rot_zy[2]); float xypp = sqrt(xypp2); float xzpp = sqrt(xzpp2); float yzpp = sqrt(yzpp2); float xyzpp = sqrt(xyzpp2); - + float normal_r = (radius / xyzpp); float r3 = normal_r * normal_r * normal_r; @@ -70,7 +70,7 @@ bfield_info_type get_dipole(float lon, // surface. But, to simplify things to begin with (and make it // planet agnostic), we use the classic definition of L-Shell, which // is with respect to the planetary radius. - + float cos_lat = xypp/xyzpp; float lShell = 1.0 / normal_r / (cos_lat * cos_lat); @@ -85,17 +85,18 @@ bfield_info_type get_dipole(float lon, float b[3]; b[0] = dipole_strength * r3 * 3 * pos_rot_zy[0] * pos_rot_zy[2] / xyzpp2; b[1] = dipole_strength * r3 * 3 * pos_rot_zy[2] * pos_rot_zy[1] / xyzpp2; - b[2] = dipole_strength * r3 / xyzpp2 * (2 * pos_rot_zy[2]*pos_rot_zy[2] - xypp2); - + b[2] = dipole_strength * r3 / xyzpp2 * + (2 * pos_rot_zy[2]*pos_rot_zy[2] - xypp2); + float b_rot_y[3]; float b_rot_yz[3]; - float b_env[3]; // env = East, North, Vertical + float b_env[3]; // env = East, North, Vertical transform_rot_y(b, magnetic_pole_tilt, b_rot_y); transform_rot_z(b_rot_y, magnetic_pole_rotation, b_rot_yz); transform_vector_xyz_to_env(b_rot_yz, lon, lat, b_env); - + bfield_info.b[0] = b_env[0]; bfield_info.b[1] = b_env[1]; bfield_info.b[2] = b_env[2]; diff --git a/src/euv.cpp b/src/euv.cpp index 02f087db..a3a1d188 100644 --- a/src/euv.cpp +++ b/src/euv.cpp @@ -1,10 +1,10 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include #include #include -#include +#include #include #include "../include/constants.h" @@ -35,14 +35,14 @@ Euv::Euv(Inputs args, Report report) { // This means we found both long and short wavelengths: if (!iErr) { for (int iWave = 0; iWave < nWavelengths; iWave++) { - ave = (wavelengths_short[iWave] + wavelengths_long[iWave])/2.0; - wavelengths_energy.push_back(planck_constant * - speed_light / - (ave * atom)); - // We simply want to initialize these vectors to make them the - // correct lenght: - wavelengths_intensity_1au.push_back(0.0); - wavelengths_intensity_top.push_back(0.0); + ave = (wavelengths_short[iWave] + wavelengths_long[iWave])/2.0; + wavelengths_energy.push_back(planck_constant * + speed_light / + (ave * atom)); + // We simply want to initialize these vectors to make them the + // correct lenght: + wavelengths_intensity_1au.push_back(0.0); + wavelengths_intensity_top.push_back(0.0); } } @@ -75,63 +75,56 @@ int Euv::read_file(Inputs args, Report report) { std::cout << "Could not open euv file!\n"; iErr = 1; } else { - nLines = 0; if (infile_ptr.good()) { - int IsFirstTime = 1; - while (getline(infile_ptr,line)) { - - report.print(5, line); - std::stringstream ss(line); - - // This is just to count the number of wavelengths. - // We assume that all of the lines have the same number of wavelengths. - if (IsFirstTime) { - std::stringstream ssdummy(line); - nWavelengths = 0; - while (getline(ssdummy, col, ',')) { - nWavelengths++; - } - // There are 6 extra items in each line: - nWavelengths -= 6; - } - - getline(ss, tmp.name, ','); - report.print(5, tmp.name); - getline(ss, tmp.to, ','); - getline(ss, tmp.type, ','); - getline(ss, col, ','); - mulfac = stof(col); - getline(ss, tmp.units, ','); - - for (int iWavelength=0; iWavelength < nWavelengths; iWavelength++) { - getline(ss, col, ','); - if (IsFirstTime) tmp.values.push_back(stof(col) * mulfac); - else tmp.values[iWavelength] = stof(col) * mulfac; - } - getline(ss, tmp.note, ','); - - waveinfo.push_back(tmp); - nLines++; - IsFirstTime = 0; - + while (getline(infile_ptr, line)) { + report.print(5, line); + std::stringstream ss(line); + + // This is just to count the number of wavelengths. + // We assume that all of the lines have the same number of wavelengths. + if (IsFirstTime) { + std::stringstream ssdummy(line); + nWavelengths = 0; + while (getline(ssdummy, col, ',')) { + nWavelengths++; + } + // There are 6 extra items in each line: + nWavelengths -= 6; + } + + getline(ss, tmp.name, ','); + report.print(5, tmp.name); + getline(ss, tmp.to, ','); + getline(ss, tmp.type, ','); + getline(ss, col, ','); + mulfac = stof(col); + getline(ss, tmp.units, ','); + + for (int iWavelength=0; iWavelength < nWavelengths; iWavelength++) { + getline(ss, col, ','); + if (IsFirstTime) tmp.values.push_back(stof(col) * mulfac); + else + tmp.values[iWavelength] = stof(col) * mulfac; + } + getline(ss, tmp.note, ','); + + waveinfo.push_back(tmp); + nLines++; + IsFirstTime = 0; } } else { - iErr = 1; - } infile_ptr.close(); - } return iErr; - } // --------------------------------------------------------------------------- @@ -139,9 +132,9 @@ int Euv::read_file(Inputs args, Report report) { // --------------------------------------------------------------------------- int Euv::slot_euv(std::string item, - std::string item2, - std::vector &values, - Report report) { + std::string item2, + std::vector &values, + Report report) { int iErr = 0; int iLine; @@ -162,7 +155,6 @@ int Euv::slot_euv(std::string item, if (iLine >= nLines) { iErr = 1; } else { - if (report.test_verbose(2)) { std::cout << "Found : " << waveinfo[iLine].name; if (!IgnoreItem2) std::cout << " with " << waveinfo[iLine].to; @@ -173,11 +165,8 @@ int Euv::slot_euv(std::string item, for (int iWavelength=0; iWavelength < nWavelengths; iWavelength++) { values.push_back(waveinfo[iLine].values[iWavelength]); } - } - return iErr; - } @@ -187,17 +176,13 @@ int Euv::slot_euv(std::string item, // -------------------------------------------------------------------------- int Euv::scale_from_1au(Planets planet, - Times time) { - + Times time) { int iErr = 0; float d = planet.get_star_to_planet_dist(time); float scale = 1.0 / (d*d); - for (int iWave = 0; iWave < nWavelengths; iWave++) wavelengths_intensity_top[iWave] = scale * wavelengths_intensity_1au[iWave]; - return iErr; - } // -------------------------------------------------------------------------- @@ -205,51 +190,46 @@ int Euv::scale_from_1au(Planets planet, // -------------------------------------------------------------------------- int Euv::euvac(Times time, - Indices indices, - Report &report) { + Indices indices, + Report &report) { int iErr = 0; float slope; - std::string function="Euv::euvac"; + std::string function = "Euv::euvac"; static int iFunction = -1; - report.enter(function, iFunction); - + report.enter(function, iFunction); + float f107 = indices.get_f107(time.get_current()); float f107a = indices.get_f107a(time.get_current()); f107 = 100.0; f107a = 100.0; - + float mean_f107 = (f107 + f107a)/2.0; if (report.test_verbose(7)) std::cout << "F107 & F107a : " << f107 << " " << f107a << "\n"; for (int iWave = 0; iWave < nWavelengths; iWave++) { - slope = 1.0 + euvac_afac[iWave] * (mean_f107 - 80.0); if (slope < 0.8) slope = 0.8; wavelengths_intensity_1au[iWave] = euvac_f74113[iWave] * slope * pcm2topm2; - } if (report.test_verbose(8)) { - std::cout << "EUVAC output : " - << f107 << " " << f107a - << " -> " << mean_f107 << "\n"; + << f107 << " " << f107a + << " -> " << mean_f107 << "\n"; for (int iWave = 0; iWave < nWavelengths; iWave++) { std::cout << " " << iWave << " " - << wavelengths_short[iWave] << " " - << wavelengths_long[iWave] << " " - << wavelengths_intensity_1au[iWave] << "\n"; + << wavelengths_short[iWave] << " " + << wavelengths_long[iWave] << " " + << wavelengths_intensity_1au[iWave] << "\n"; } - } - report.exit(function); + report.exit(function); return iErr; - } diff --git a/src/file_input.cpp b/src/file_input.cpp index 82159d9a..891d2205 100644 --- a/src/file_input.cpp +++ b/src/file_input.cpp @@ -1,10 +1,10 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include #include #include -#include +#include #include // ------------------------------------------------------------------- @@ -13,13 +13,11 @@ // ------------------------------------------------------------------- std::string make_lower(std::string instring) { - - std::string outstring=instring; - for (int i = 0; i < outstring.length(); i++) + std::string outstring = instring; + int len = instring.length(); + for (int i = 0; i < len; i++) outstring[i] = tolower(outstring[i]); - return outstring; - } // ------------------------------------------------------------------- @@ -28,23 +26,22 @@ std::string make_lower(std::string instring) { // ------------------------------------------------------------------- std::string strip_string_end(std::string instring) { - int i; int iStart = -1; std::string outstring; - for (i = 0; i < instring.length(); i++) { + int len = instring.length(); + for (i = 0; i < len; i++) { if (instring[i] == ' ') { if (i > 0 && iStart == i-1) break; - else iStart = i; + else + iStart = i; } } - - if (iStart > -1 && i < instring.length()) - outstring = instring.substr(0,iStart); - else outstring = instring; - + if (iStart > -1 && i < len) + outstring = instring.substr(0, iStart); + else + outstring = instring; return outstring; - } // ------------------------------------------------------------------- @@ -52,23 +49,20 @@ std::string strip_string_end(std::string instring) { // ------------------------------------------------------------------- std::string strip_spaces(std::string instring) { - int i; - int iStart = -1; std::string outstring; - for (i = 0; i < instring.length(); i++) { + int len = instring.length(); + for (i = 0; i < len; i++) { if (instring[i] != ' ') { outstring = outstring+instring[i]; } } - return outstring; - } // ------------------------------------------------------------------- -// -// +// +// // ------------------------------------------------------------------- // ------------------------------------------------------------------- @@ -77,24 +71,22 @@ std::string strip_spaces(std::string instring) { std::string read_string(std::ifstream &file_ptr, std::string hash) { - std::string line=""; + std::string line = ""; if (!file_ptr.is_open()) { std::cout << "File is not open (read_string)!\n"; std::cout << "hash : " << hash << "\n"; } else { - - getline(file_ptr,line); + getline(file_ptr, line); // Clearly these could all be combined, but I like separating them: line = strip_string_end(line); line = strip_spaces(line); - if (line.length() < 1) { std::cout << "Issue in read_inputs!\n"; std::cout << "Should be:\n"; std::cout << hash << "\n"; std::cout << "string variable with length > 0\n"; - } + } } return line; } @@ -105,18 +97,15 @@ std::string read_string(std::ifstream &file_ptr, std::string hash) { int read_int(std::ifstream &file_ptr, std::string hash) { - std::string line=""; + std::string line = ""; int output = -1; - + if (!file_ptr.is_open()) { std::cout << "File is not open (read_int)!\n"; std::cout << "hash : " << hash << "\n"; } else { - - getline(file_ptr,line); - // Clearly these could all be combined, but I like separating them: + getline(file_ptr, line); line = strip_string_end(line); - try { output = stoi(line); } @@ -126,7 +115,7 @@ int read_int(std::ifstream &file_ptr, std::string hash) { std::cout << hash << "\n"; std::cout << "Trying to read an integer, but got this: "; std::cout << line << "\n"; - } + } } return output; } @@ -142,17 +131,14 @@ std::string find_next_hash(std::ifstream &file_ptr) { if (!file_ptr.is_open()) { std::cout << "File is not open (find_next_hash)!\n"; } else { - - while (getline(file_ptr,line)) { + while (getline(file_ptr, line)) { // Clearly these could all be combined, but I like separating them: line = strip_string_end(line); line = strip_spaces(line); line = make_lower(line); if (line[0] == '#') return line; } - } - line = ""; return line; } @@ -169,40 +155,35 @@ std::vector> read_csv(std::ifstream &file_ptr) { std::vector> data; std::vector row; - std::string line, col; - line = " "; if (!file_ptr.is_open()) { std::cout << "File is not open (read_csv)!\n"; } else { - // This assumes that the CSV file's layout is perfect - that the // number of columns is the same in each row. If that is not the // case, then bad stuff happens. I need to add more debugging // stuff in here. - int IsFirstTime = 1; - while (getline(file_ptr,line) && line.length() > 1) { + while (getline(file_ptr, line) && line.length() > 1) { line = strip_string_end(line); line = strip_spaces(line); std::stringstream ss(line); int j = 0; while (getline(ss, col, ',')) { - if (IsFirstTime) { - row.push_back(col); - } else row[j]=col; - j++; + if (IsFirstTime) { + row.push_back(col); + } else { + row[j] = col; + } + j++; } data.push_back(row); IsFirstTime = 0; } - } - return data; - } // ------------------------------------------------------------------- @@ -211,17 +192,21 @@ std::vector> read_csv(std::ifstream &file_ptr) { std::vector read_itime(std::ifstream &file_ptr, std::string hash) { int iErr = 0; - std::vector itime(7,0); + std::vector itime(7, 0); // Read in the series of numbers: std::vector> stimes = read_csv(file_ptr); // Interpret whether it is a column (first) or row (second): - if (stimes.size() == 6) // column - for (int i=0; i<6; i++) itime[i] = stoi(stimes[i][0]); - else if (stimes.size() == 1 && stimes[0].size() >= 6) // row - for (int i=0; i<6; i++) itime[i] = stoi(stimes[0][i]); - else iErr = 1; + if (stimes.size() == 6) { // column + for (int i = 0; i < 6; i++) itime[i] = stoi(stimes[i][0]); + } else { + if (stimes.size() == 1 && stimes[0].size() >= 6) { // row + for (int i = 0; i < 6; i++) itime[i] = stoi(stimes[0][i]); + } else { + iErr = 1; + } + } if (iErr == 0) { itime[6] = 0; @@ -242,8 +227,5 @@ std::vector read_itime(std::ifstream &file_ptr, std::string hash) { std::cout << hash << "\n"; std::cout << "year, month, day, hour, minute, second\n"; } - return itime; - } - diff --git a/src/fill_grid.cpp b/src/fill_grid.cpp index d05268e5..6278fb98 100644 --- a/src/fill_grid.cpp +++ b/src/fill_grid.cpp @@ -1,4 +1,4 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include @@ -19,38 +19,21 @@ void Grid::calc_sza(Planets planet, Times time, Report &report) { - long iLon, iLat, iAlt, index; - float local_time; - std::string function = "Grid::calc_sza"; static int iFunction = -1; - report.enter(function, iFunction); + report.enter(function, iFunction); float lon_offset = planet.get_longitude_offset(time); float sin_dec = planet.get_sin_dec(time); float cos_dec = planet.get_cos_dec(time); - for (iLon = 0; iLon < nGeoLonsG; iLon++) { - for (iLat = 0; iLat < nGeoLatsG; iLat++) { - for (iAlt = 0; iAlt < nGeoAltsG; iAlt++) { - - index = ijk_geo_s3gc(iLon,iLat,iAlt); - - // This is in radians: - local_time = fmod(lon_offset + geoLon_s3gc[index] + twopi, twopi); - - cos_sza_s3gc[index] = - sin_dec * sin(geoLat_s3gc[index]) + - cos_dec * cos(geoLat_s3gc[index]) * cos(local_time-pi); - - sza_s3gc[index] = acos(cos_sza_s3gc[index]); - - } - } - } + fcube local_time3d = geoLon_scgc + lon_offset; + cos_sza_scgc = + sin_dec * sin(geoLat_scgc) + + cos_dec * cos(geoLat_scgc) % cos(local_time3d-pi); + sza_scgc = acos(cos_sza_scgc); report.exit(function); - } // ----------------------------------------------------------------------------- @@ -61,54 +44,30 @@ void Grid::fill_grid_bfield(Planets planet, Inputs input, Report &report) { std::string function = "Grid::fill_grid_bfield"; static int iFunction = -1; - report.enter(function, iFunction); + report.enter(function, iFunction); - long nLons, nLats, nAlts, iLon, iLat, iAlt, iDim, index, indexv; + int64_t iLon, iLat, iAlt, iDim; float lon, lat, alt; bfield_info_type bfield_info; - - if (IsGeoGrid) { - nLons = nGeoLonsG; - nLats = nGeoLatsG; - nAlts = nGeoAltsG; - } else { - nLons = nMagLonsG; - nLats = nMagLatsG; - nAlts = nMagAltsG; - } for (iLon = 0; iLon < nLons; iLon++) { for (iLat = 0; iLat < nLats; iLat++) { for (iAlt = 0; iAlt < nAlts; iAlt++) { - if (IsGeoGrid) { - index = ijk_geo_s3gc(iLon,iLat,iAlt); - } else { - index = ijk_mag_s3gc(iLon,iLat,iAlt); - } - - lon = geoLon_s3gc[index]; - lat = geoLat_s3gc[index]; - alt = geoAlt_s3gc[index]; - - bfield_info = get_bfield(lon, lat, alt, planet, input, report); - - magLat_s3gc[index] = bfield_info.lat; - magLon_s3gc[index] = bfield_info.lon; - - for (iDim = 0; iDim < 3; iDim++) { - if (IsGeoGrid) { - indexv = ijkl_geo_v3gc(iLon,iLat,iAlt,iDim); - } else { - indexv = ijkl_mag_v3gc(iLon,iLat,iAlt,iDim); - } - bfield_v3gc[indexv] = bfield_info.b[iDim]; - } - + lon = geoLon_scgc(iLon, iLat, iAlt); + lat = geoLat_scgc(iLon, iLat, iAlt); + alt = geoAlt_scgc(iLon, iLat, iAlt); + + bfield_info = get_bfield(lon, lat, alt, planet, input, report); + + magLat_scgc(iLon, iLat, iAlt) = bfield_info.lat; + magLon_scgc(iLon, iLat, iAlt) = bfield_info.lon; + + for (iDim = 0; iDim < 3; iDim++) + bfield_vcgc[iDim](iLon, iLat, iAlt) = bfield_info.b[iDim]; } } } - report.exit(function); return; } @@ -119,42 +78,29 @@ void Grid::fill_grid_bfield(Planets planet, Inputs input, Report &report) { void Grid::fill_grid_radius(Planets planet, Report &report) { - long iLon, iLat, iAlt, index; + int64_t iLon, iLat, iAlt; float radius0; float mu = planet.get_mu(); report.print(3, "starting fill_grid_radius"); - - for (iLon = 0; iLon < nGeoLonsG; iLon++) { - for (iLat = 0; iLat < nGeoLatsG; iLat++) { - for (iAlt = 0; iAlt < nGeoAltsG; iAlt++) { - - if (IsGeoGrid) { - index = ijk_geo_s3gc(iLon,iLat,iAlt); - } else { - index = ijk_mag_s3gc(iLon,iLat,iAlt); - } - // Set radius and radius^2: + // Just in case we have a latitude-dependent planetary radius + fvec radius0_1d(nLats); + for (iLat = 0; iLat < nLats; iLat++) + radius0_1d(iLat) = planet.get_radius(geoLat_scgc(0, iLat, 0)); - radius0 = planet.get_radius(geoLat_s3gc[index]); + for (iLon=0; iLon < nLons; iLon++) + for (iAlt=0; iAlt < nAlts; iAlt++) + radius_scgc.subcube(iLon, 0, iAlt, iLon, nLats-1, iAlt) = radius0_1d; + radius_scgc = radius_scgc + geoAlt_scgc; - radius_s3gc[index] = - radius0 + geoAlt_s3gc[index]; - radius_sq_s3gc[index] = - radius_s3gc[index] * radius_s3gc[index]; - radius_inv_sq_s3gc[index] = - 1.0/radius_sq_s3gc[index]; + radius2_scgc = radius_scgc % radius_scgc; + radius2i_scgc = 1.0 / radius2_scgc; - gravity_s3gc[index] = mu * radius_inv_sq_s3gc[index]; + gravity_scgc = mu * radius2i_scgc; - } - } - } - report.print(3, "ending fill_grid_radius"); - } // ----------------------------------------------------------------------------- @@ -163,114 +109,26 @@ void Grid::fill_grid_radius(Planets planet, Report &report) { void Grid::fill_grid(Planets planet, Report &report) { - long iLon, iLat, iAlt, index; + int64_t iLon, iLat, iAlt; - long indexp, indexm; - long nLons, nLats, nAlts; - report.print(3, "starting fill_grid"); - - if (IsGeoGrid) { - nLons = nGeoLonsG; - nLats = nGeoLatsG; - nAlts = nGeoAltsG; - } else { - nLons = nMagLonsG; - nLats = nMagLatsG; - nAlts = nMagAltsG; - } float xyz[3], llr[3]; - - for (iLon = 0; iLon < nLons; iLon++) { - for (iLat = 0; iLat < nLats; iLat++) { - for (iAlt = 0; iAlt < nAlts; iAlt++) { - - if (IsGeoGrid) { - index = ijk_geo_s3gc(iLon,iLat,iAlt); - } else { - index = ijk_mag_s3gc(iLon,iLat,iAlt); - } - - // Find XYZ coordinates (Geo): - - llr[0] = geoLon_s3gc[index]; - llr[1] = geoLat_s3gc[index]; - llr[2] = geoAlt_s3gc[index]; - - transform_llr_to_xyz(llr, xyz); - - geoX_s3gc[index] = xyz[0]; - geoY_s3gc[index] = xyz[1]; - geoZ_s3gc[index] = xyz[2]; - - //// Find XYZ coordinates (Mag): - // - //llr[0] = magLon_s3gc[index]; - //llr[1] = magLat_s3gc[index]; - //llr[2] = magAlt_s3gc[index]; - // - //transform_llr_to_xyz(llr, xyz); - // - //magX_s3gc[index] = xyz[0]; - //magY_s3gc[index] = xyz[1]; - //magZ_s3gc[index] = xyz[2]; - - if (iAlt > 0 && iAlt < nAlts-1) { - if (IsGeoGrid) { - indexp = ijk_geo_s3gc(iLon,iLat,iAlt+1); - indexm = ijk_geo_s3gc(iLon,iLat,iAlt-1); - } else { - indexp = ijk_mag_s3gc(iLon,iLat,iAlt+1); - indexm = ijk_mag_s3gc(iLon,iLat,iAlt-1); - } - - dalt_center_s3gc[index] = - (geoAlt_s3gc[indexp] - geoAlt_s3gc[indexm])/2; - dalt_lower_s3gc[index] = - geoAlt_s3gc[index] - geoAlt_s3gc[indexm]; - - } - - } - } - } - - report.print(4, "filling edges"); - - for (iLon = 0; iLon < nLons; iLon++) { - for (iLat = 0; iLat < nLats; iLat++) { - - iAlt = 0; - - if (IsGeoGrid) { - index = ijk_geo_s3gc(iLon,iLat,iAlt); - indexp = ijk_geo_s3gc(iLon,iLat,iAlt+1); - } else { - index = ijk_mag_s3gc(iLon,iLat,iAlt); - indexp = ijk_mag_s3gc(iLon,iLat,iAlt+1); - } - - dalt_center_s3gc[index] = dalt_center_s3gc[indexp]; - dalt_lower_s3gc[index] = dalt_lower_s3gc[indexp]; - iAlt = nAlts-1; - - if (IsGeoGrid) { - index = ijk_geo_s3gc(iLon,iLat,iAlt); - indexm = ijk_geo_s3gc(iLon,iLat,iAlt-1); - } else { - index = ijk_mag_s3gc(iLon,iLat,iAlt); - indexm = ijk_mag_s3gc(iLon,iLat,iAlt-1); - } - - dalt_center_s3gc[index] = dalt_center_s3gc[indexm]; - dalt_lower_s3gc[index] = dalt_lower_s3gc[indexm]; - - } + for (iAlt = 1; iAlt < nAlts-1; iAlt++) { + dalt_center_scgc.slice(iAlt) = + geoAlt_scgc.slice(iAlt+1) - geoAlt_scgc.slice(iAlt-1); + dalt_lower_scgc.slice(iAlt) = + geoAlt_scgc.slice(iAlt) - geoAlt_scgc.slice(iAlt-1); } + dalt_center_scgc.slice(0) = dalt_center_scgc.slice(1); + dalt_center_scgc.slice(nAlts-1) = dalt_center_scgc.slice(nAlts-2); - report.print(3, "ending fill_grid"); + dalt_lower_scgc.slice(0) = dalt_lower_scgc.slice(1); + iAlt = nAlts-1; + dalt_lower_scgc.slice(iAlt) = + geoAlt_scgc.slice(iAlt) - geoAlt_scgc.slice(iAlt-1); - + transform_llr_to_xyz_3d(geoLon_scgc, geoLat_scgc, radius_scgc, + geoX_scgc, geoY_scgc, geoZ_scgc); } diff --git a/src/grid.cpp b/src/grid.cpp index 7c675d57..a15c9fb9 100644 --- a/src/grid.cpp +++ b/src/grid.cpp @@ -1,47 +1,59 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include +#include #include "../include/inputs.h" #include "../include/grid.h" #include "../include/sizes.h" -Grid::Grid(int nX, int nY, int nZ) { +using namespace arma; - long nTotalPoints = long(nX) * long(nY) * long(nZ); +Grid::Grid(int nX_in, int nY_in, int nZ_in, int nGCs_in) { - geoLon_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); - geoLat_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); - geoAlt_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); - geoX_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); - geoY_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); - geoZ_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); + nX = nX_in; nLons = nX; + nY = nY_in; nLats = nY; + nZ = nZ_in; nAlts = nZ; + nGCs = nGCs_in; - magLon_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); - magLat_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); - magAlt_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); - magX_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); - magY_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); - magZ_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); + int64_t nTotalPoints = int64_t(nX) * int64_t(nY) * int64_t(nZ); - magLocalTime_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); + geoLon_scgc.set_size(nX, nY, nZ); + geoLat_scgc.set_size(nX, nY, nZ); + geoAlt_scgc.set_size(nX, nY, nZ); - radius_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); - radius_sq_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); - radius_inv_sq_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); + geoX_scgc.set_size(nX, nY, nZ); + geoY_scgc.set_size(nX, nY, nZ); + geoZ_scgc.set_size(nX, nY, nZ); - gravity_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); + magLon_scgc.set_size(nX, nY, nZ); + magLat_scgc.set_size(nX, nY, nZ); + magAlt_scgc.set_size(nX, nY, nZ); - sza_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); - cos_sza_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); + magX_scgc.set_size(nX, nY, nZ); + magY_scgc.set_size(nX, nY, nZ); + magZ_scgc.set_size(nX, nY, nZ); - bfield_v3gc = (float*) malloc( 3 * nTotalPoints * sizeof(float) ); - bfield_mag_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); + magLocalTime_scgc.set_size(nX, nY, nZ); - dalt_center_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); - dalt_lower_s3gc = (float*) malloc( nTotalPoints * sizeof(float) ); - + radius_scgc.set_size(nX, nY, nZ); + radius2_scgc.set_size(nX, nY, nZ); + radius2i_scgc.set_size(nX, nY, nZ); + + dalt_center_scgc.set_size(nX, nY, nZ); + dalt_lower_scgc.set_size(nX, nY, nZ); + + sza_scgc.set_size(nX, nY, nZ); + cos_sza_scgc.set_size(nX, nY, nZ); + + fcube tmp(nX, nY, nZ); + tmp.zeros(); + bfield_vcgc.push_back(tmp); // x-component + bfield_vcgc.push_back(tmp); // y-component + bfield_vcgc.push_back(tmp); // z-component + bfield_mag_scgc.set_size(nX, nY, nZ); + bfield_mag_scgc.zeros(); } int Grid::get_IsGeoGrid() { @@ -52,11 +64,18 @@ void Grid::set_IsGeoGrid(int value) { IsGeoGrid = value; } -long Grid::get_nPointsInGrid() { - long nPoints; - if (IsGeoGrid) - nPoints = long(nGeoLons) * long(nGeoLats) * long(nGeoAlts); - else - nPoints = long(nMagLons) * long(nMagLats) * long(nMagAlts); +int64_t Grid::get_nPointsInGrid() { + int64_t nPoints; + nPoints = int64_t(nX) * int64_t(nY) * int64_t(nZ); return nPoints; } + +int64_t Grid::get_nX() { return nX; } +int64_t Grid::get_nY() { return nY; } +int64_t Grid::get_nZ() { return nZ; } + +int64_t Grid::get_nLons() { return nLons; } +int64_t Grid::get_nLats() { return nLats; } +int64_t Grid::get_nAlts() { return nAlts; } + +int64_t Grid::get_nGCs() { return nGCs; } diff --git a/src/indices.cpp b/src/indices.cpp index ce2c659d..0f8e0f6c 100644 --- a/src/indices.cpp +++ b/src/indices.cpp @@ -1,4 +1,4 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include @@ -25,13 +25,13 @@ Indices::Indices(Inputs args) { iErr = read_f107_file(file, time, f107array); if (iErr == 0) iErr = set_f107(time, f107array); - else std::cout << "ERROR in reading f107 file!!!\n"; + else + std::cout << "ERROR in reading f107 file!!!\n"; } - } int Indices::set_f107(std::vector time, - std::vector f107array) { + std::vector f107array) { int iErr = set_index(time, f107array, f107); @@ -44,14 +44,14 @@ int Indices::set_f107(std::vector time, // hours at a time double currenttime = time[0]; - long nTimes = time.size(), itime = 0; + int64_t nTimes = time.size(), itime = 0; double eightone = 81.0 * 86400.0; ind_time_pair tmp; while (currenttime < time[nTimes-1]-eightone) { - long isub = itime, nSubs = 0; + int64_t isub = itime, nSubs = 0; double sumf107 = 0, sumtime = 0; while (time[isub] < currenttime + eightone) { @@ -68,7 +68,6 @@ int Indices::set_f107(std::vector time, itime++; currenttime = time[itime]; - } // Let's hold the last 81 days constant, which means we just put one @@ -79,24 +78,19 @@ int Indices::set_f107(std::vector time, f107a.push_back(tmp); return iErr; - } float Indices:: get_f107(double time) { - return get_index(time, f107); - } float Indices:: get_f107a(double time) { - return get_index(time, f107a); - } float Indices::get_index(double time, std::vector index) { - long iLow, iMid, iHigh; + int64_t iLow, iMid, iHigh; iLow = 0; iHigh = index.size()-1; @@ -104,8 +98,6 @@ float Indices::get_index(double time, std::vector index) { while (iHigh-iLow > 1) { - // cout << "top: " << iLow << " " << iMid << " " << iHigh << " " << index[iMid].time-time <<"\n"; - // Break if iMid <= time < iMid+1 if (index[iMid].time == time) break; if (index[iMid].time <= time && index[iMid+1].time > time) break; @@ -117,9 +109,6 @@ float Indices::get_index(double time, std::vector index) { iHigh = iMid; iMid = (iHigh+iLow)/2; } - - // cout << "bot: " << iLow << " " << iMid << " " << iHigh << " " << index[iMid].time-time <<"\n"; - } // At this point, time should be between iMid and iMid+1: @@ -128,12 +117,11 @@ float Indices::get_index(double time, std::vector index) { float value = (1.0-x) * index[iMid].index + x * index[iMid+1].index; return value; - } int Indices::set_index(std::vector time, - std::vector indexarray, - std::vector &index) { + std::vector indexarray, + std::vector &index) { int iErr; ind_time_pair tmp; @@ -142,15 +130,13 @@ int Indices::set_index(std::vector time, std::cout << "In set_index. Size of time and index arrays don't match!\n"; iErr = 1; } else { - - for (int i=0; i +#include #include "../include/inputs.h" #include "../include/report.h" @@ -10,55 +11,70 @@ #include "../include/sizes.h" #include "../include/fill_grid.h" +using namespace arma; + void Grid::init_geo_grid(Planets planet, Inputs input, Report &report) { - std::string function="Grid::init_geo_grid"; + std::string function = "Grid::init_geo_grid"; static int iFunction = -1; - report.enter(function, iFunction); + report.enter(function, iFunction); // This is just an example: Inputs::grid_input_struct grid_input = input.get_grid_inputs(); - long iLon, iLat, iAlt, index; + int64_t iLon, iLat, iAlt, index; float longitude, latitude, altitude; - float altitudes[nGeoAltsG]; + // Longitudes: + // - Make a 1d vector + // - copy it into the 3d cube + fvec lon1d(nLons); + float dlon = (grid_input.lon_max - grid_input.lon_min) / nGeoLons; + for (iLon=0; iLon < nLons; iLon++) + lon1d(iLon) = (iLon-nGCs+0.5) * dlon; + for (iLat=0; iLat < nLats; iLat++) { + for (iAlt=0; iAlt < nAlts; iAlt++) { + geoLon_scgc.subcube(0, iLat, iAlt, nLons-1, iLat, iAlt) = lon1d; + } + } + + // Latitudes: + // - Make a 1d vector + // - copy it into the 3d cube + fvec lat1d(nLats); float dlat = (grid_input.lat_max - grid_input.lat_min) / nGeoLats; - float dlon = (grid_input.lon_max - grid_input.lon_min) / nGeoLons; + for (iLat=0; iLat < nLats; iLat++) + lat1d(iLat) = grid_input.lat_min + (iLat-nGCs+0.5) * dlat; - IsGeoGrid = 1; - - // Make a uniform grid in altitude: - if (grid_input.IsUniformAlt) { - for (iAlt=0; iAlt < nGeoAltsG; iAlt++) { - altitudes[iAlt] = - grid_input.alt_min + float(iAlt-nGeoGhosts)*grid_input.dalt; + for (iLon=0; iLon < nLons; iLon++) { + for (iAlt=0; iAlt < nAlts; iAlt++) { + geoLat_scgc.subcube(iLon, 0, iAlt, iLon, nLats-1, iAlt) = lat1d; } - } else { - // calc_stretched_altitudes(input, planet, altitudes); } - for (iLon = 0; iLon < nGeoLonsG; iLon++) { - longitude = grid_input.lon_min + (float(iLon-nGeoGhosts)+0.5) * dlon; - for (iLat = 0; iLat < nGeoLatsG; iLat++) { - latitude = grid_input.lat_min + (float(iLat-nGeoGhosts)+0.5) * dlat; - for (iAlt = 0; iAlt < nGeoAltsG; iAlt++) { - index = ijk_geo_s3gc(iLon,iLat,iAlt); - geoLon_s3gc[index] = longitude; - geoLat_s3gc[index] = latitude; - geoAlt_s3gc[index] = altitudes[iAlt]; - } + fvec alt1d(nAlts); + + if (grid_input.IsUniformAlt) { + for (iAlt = 0; iAlt < nAlts; iAlt++) + alt1d(iAlt) = + grid_input.alt_min + + (iAlt-nGeoGhosts) * grid_input.dalt; + } + for (iLon = 0; iLon < nLons; iLon++) { + for (iLat = 0; iLat < nLats; iLat++) { + geoAlt_scgc.tube(iLon, iLat) = alt1d; } } + IsGeoGrid = 1; + // Calculate the radius, etc: fill_grid_radius(planet, report); fill_grid_bfield(planet, input, report); - - report.exit(function); + report.exit(function); } diff --git a/src/inputs.cpp b/src/inputs.cpp index 7aba2c8f..c50d7c2e 100644 --- a/src/inputs.cpp +++ b/src/inputs.cpp @@ -1,3 +1,5 @@ +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) +// Full license can be found in License.md #include #include @@ -14,15 +16,11 @@ Inputs::Inputs(Times &time, Report &report) { - int iErr; - - iErr = 0; - // ------------------------------------------------ // Set some defaults: iVerbose = 0; - euv_model="euvac"; + euv_model = "euvac"; planet = "Earth"; // ------------------------------------------------ @@ -30,11 +28,7 @@ Inputs::Inputs(Times &time, Report &report) { grid_input.alt_file = ""; grid_input.IsUniformAlt = 1; grid_input.alt_min = 100.0 * 1000.0; - grid_input.dalt = 2.5 * 1000.0; - - //grid_input.IsUniformAlt = 0; - //grid_input.alt_min = 100.0 * 1000.0; - //grid_input.dalt = 0.33; + grid_input.dalt = 5.0 * 1000.0; if (nGeoLons == 1) { grid_input.lon_min = 0.0; @@ -61,8 +55,7 @@ Inputs::Inputs(Times &time, Report &report) { // ------------------------------------------------ // Now read the input file: - iErr = read(time, report); - + int iErr = read(time, report); } // ----------------------------------------------------------------------- @@ -119,7 +112,8 @@ float Inputs::get_n_outputs() { float Inputs::get_dt_output(int iOutput) { float value = 0.0; - if (iOutput < dt_output.size()) value = dt_output[iOutput]; + int iSize = dt_output.size(); + if (iOutput < iSize) value = dt_output[iOutput]; return value; } @@ -129,7 +123,8 @@ float Inputs::get_dt_output(int iOutput) { std::string Inputs::get_type_output(int iOutput) { std::string value = ""; - if (iOutput < dt_output.size()) value = type_output[iOutput]; + int iSize = dt_output.size(); + if (iOutput < iSize) value = type_output[iOutput]; return value; } @@ -189,7 +184,7 @@ int Inputs::read(Times &time, Report &report) { int iErr; std::string line, hash; - std::vector itime(7,0); + std::vector itime(7, 0); iErr = 0; @@ -209,15 +204,15 @@ int Inputs::read(Times &time, Report &report) { hash = find_next_hash(infile_ptr); if (report.test_verbose(3)) - std::cout << "hash : -->" << hash << "<--\n"; + std::cout << "hash : -->" << hash << "<--\n"; // --------------------------- // #debug or #verbose // --------------------------- if (hash == "#debug" || hash == "#verbose") { - iVerbose = read_int(infile_ptr, hash); - report.set_verbose(iVerbose); + iVerbose = read_int(infile_ptr, hash); + report.set_verbose(iVerbose); } // --------------------------- @@ -225,12 +220,12 @@ int Inputs::read(Times &time, Report &report) { // --------------------------- if (hash == "#starttime") { - std::vector istart = read_itime(infile_ptr, hash); - if (istart[0] > 0) time.set_times(istart); - if (report.test_verbose(3)) { - std::cout << "Starttime : "; - display_itime(istart); - } + std::vector istart = read_itime(infile_ptr, hash); + if (istart[0] > 0) time.set_times(istart); + if (report.test_verbose(3)) { + std::cout << "Starttime : "; + display_itime(istart); + } } // --------------------------- @@ -238,12 +233,12 @@ int Inputs::read(Times &time, Report &report) { // --------------------------- if (hash == "#endtime") { - std::vector iend = read_itime(infile_ptr, hash); - if (iend[0] > 0) time.set_end_time(iend); - if (report.test_verbose(3)) { - std::cout << "Endtime : "; - display_itime(iend); - } + std::vector iend = read_itime(infile_ptr, hash); + if (iend[0] > 0) time.set_end_time(iend); + if (report.test_verbose(3)) { + std::cout << "Endtime : "; + display_itime(iend); + } } // --------------------------- @@ -251,7 +246,7 @@ int Inputs::read(Times &time, Report &report) { // --------------------------- if (hash == "#f107file") { - f107_file = read_string(infile_ptr, hash); + f107_file = read_string(infile_ptr, hash); } // --------------------------- @@ -259,7 +254,7 @@ int Inputs::read(Times &time, Report &report) { // --------------------------- if (hash == "#bfield") { - bfield = read_string(infile_ptr, hash); + bfield = read_string(infile_ptr, hash); } // --------------------------- @@ -267,7 +262,7 @@ int Inputs::read(Times &time, Report &report) { // --------------------------- if (hash == "#chemistry") { - chemistry_file = read_string(infile_ptr, hash); + chemistry_file = read_string(infile_ptr, hash); } // --------------------------- @@ -275,11 +270,11 @@ int Inputs::read(Times &time, Report &report) { // --------------------------- if (hash == "#planet") { - planet = read_string(infile_ptr, hash); - if (report.test_verbose(3)) - std::cout << "Setting planet to : " << planet << "\n"; - if (planet_species_file.length() <= 1) - planet_species_file = "UA/inputs/"+planet+".in"; + planet = read_string(infile_ptr, hash); + if (report.test_verbose(3)) + std::cout << "Setting planet to : " << planet << "\n"; + if (planet_species_file.length() <= 1) + planet_species_file = "UA/inputs/"+planet+".in"; } // --------------------------- @@ -287,36 +282,27 @@ int Inputs::read(Times &time, Report &report) { // --------------------------- if (hash == "#output") { - std::vector> csv = read_csv(infile_ptr); - // comma separated values, with type, then dt: - int nOutputs = csv.size(); - int iOutput; - std::cout << "output : " << nOutputs << "\n"; - if (nOutputs > 1) { - std::cout << "output0: " << type_output[0] << "\n"; - type_output[0] = csv[0][0]; - dt_output[0] = stof(csv[0][1]); - for (iOutput = 1; iOutput < nOutputs; iOutput++) { - std::cout << "output n : " << iOutput << " " << csv[iOutput][0] << "\n"; - type_output.push_back(csv[iOutput][0]); - dt_output.push_back(stof(csv[iOutput][1])); - } - // Allow users to enter 0 for dt, so they only get the - // output at the beginning of the run: - for (iOutput = 0; iOutput < nOutputs; iOutput++) - if (dt_output[iOutput] <= 0.0) dt_output[iOutput] = 1.0e32; - } else { - std::cout << "Something wrong with #output. Need to report...\n"; - } + std::vector> csv = read_csv(infile_ptr); + // comma separated values, with type, then dt: + int nOutputs = csv.size(); + int iOutput; + if (nOutputs > 1) { + type_output[0] = csv[0][0]; + dt_output[0] = stof(csv[0][1]); + for (iOutput = 1; iOutput < nOutputs; iOutput++) { + type_output.push_back(csv[iOutput][0]); + dt_output.push_back(stof(csv[iOutput][1])); + } + // Allow users to enter 0 for dt, so they only get the + // output at the beginning of the run: + for (iOutput = 0; iOutput < nOutputs; iOutput++) + if (dt_output[iOutput] <= 0.0) dt_output[iOutput] = 1.0e32; + } else { + std::cout << "Something wrong with #output. Need to report...\n"; + } } - } - infile_ptr.close(); - } - return iErr; - - } diff --git a/src/ions.cpp b/src/ions.cpp index a865db5b..2c56fe2f 100644 --- a/src/ions.cpp +++ b/src/ions.cpp @@ -1,4 +1,4 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include @@ -15,80 +15,80 @@ // ----------------------------------------------------------------------------- -// +// // ----------------------------------------------------------------------------- -Ions::species_chars Ions::create_species() { +Ions::species_chars Ions::create_species(Grid grid) { - long iDir, iLon, iLat, iAlt, index; + int64_t iDir, iLon, iLat, iAlt, index; species_chars tmp; - long iTotal = long(nGeoLonsG) * long(nGeoLatsG) * long(nGeoAltsG); + int64_t nLons = grid.get_nLons(); + int64_t nLats = grid.get_nLats(); + int64_t nAlts = grid.get_nAlts(); + + int64_t iTotal = int64_t(nLons) * int64_t(nLats) * int64_t(nAlts); // Constants: tmp.DoAdvect = 0; - - tmp.density_s3gc = (float*) malloc( iTotal * sizeof(float) ); - tmp.par_velocity_v3gc = (float*) malloc( long(3)*iTotal * sizeof(float) ); - tmp.perp_velocity_v3gc = (float*) malloc( long(3)*iTotal * sizeof(float) ); - tmp.temperature_s3gc = (float*) malloc( iTotal * sizeof(float) ); - tmp.ionization_s3gc = (float*) malloc( iTotal * sizeof(float) ); - - for (iLon = 0; iLon < nGeoLonsG; iLon++) { - for (iLat = 0; iLat < nGeoLatsG; iLat++) { - for (iAlt = 0; iAlt < nGeoAltsG; iAlt++) { - - index = ijk_geo_s3gc(iLon,iLat,iAlt); - - tmp.density_s3gc[index] = 1.0; - tmp.temperature_s3gc[index] = 1.0e-32; - tmp.ionization_s3gc[index] = 1.0e-32; - - for (iDir = 0; iDir < 3; iDir++) { - index = ijkl_geo_v3gc(iLon,iLat,iAlt,iDir); - tmp.par_velocity_v3gc[index] = 0.0; - tmp.perp_velocity_v3gc[index] = 0.0; - } - - } - } - } - + + tmp.density_scgc.set_size(nLons, nLats, nAlts); + tmp.density_scgc.ones(); + tmp.temperature_scgc.set_size(nLons, nLats, nAlts); + tmp.temperature_scgc.ones(); + tmp.ionization_scgc.set_size(nLons, nLats, nAlts); + tmp.ionization_scgc.zeros(); + + tmp.sources_scgc.set_size(nLons, nLats, nAlts); + tmp.sources_scgc.zeros(); + tmp.losses_scgc.set_size(nLons, nLats, nAlts); + tmp.losses_scgc.zeros(); + return tmp; - } // ----------------------------------------------------------------------------- -// +// Initialize ions // ----------------------------------------------------------------------------- -Ions::Ions(Inputs input, Report report) { +Ions::Ions(Grid grid, Inputs input, Report report) { + + int64_t nLons = grid.get_nLons(); + int64_t nLats = grid.get_nLats(); + int64_t nAlts = grid.get_nAlts(); + + int64_t iTotal = int64_t(nLons) * int64_t(nLats) * int64_t(nAlts); - long iTotal = long(nGeoLonsG) * long(nGeoLatsG) * long(nGeoAltsG); species_chars tmp; int iErr; - report.print(2,"Initializing Ions"); - + report.print(2, "Initializing Ions"); + for (int iSpecies=0; iSpecies < nIons; iSpecies++) { - tmp = create_species(); + tmp = create_species(grid); species.push_back(tmp); } // Create one extra species for electrons - tmp = create_species(); + tmp = create_species(grid); species.push_back(tmp); // State variables: - density_s3gc = (float*) malloc( iTotal * sizeof(float) ); - velocity_v3gc = (float*) malloc( long(3)*iTotal * sizeof(float) ); - exb_v3gc = (float*) malloc( long(3)*iTotal * sizeof(float) ); - ion_temperature_s3gc = (float*) malloc( iTotal * sizeof(float) ); - electron_temperature_s3gc = (float*) malloc( iTotal * sizeof(float) ); + + density_scgc.set_size(nLons, nLats, nAlts); + density_scgc.ones(); + ion_temperature_scgc.set_size(nLons, nLats, nAlts); + ion_temperature_scgc.ones(); + electron_temperature_scgc.set_size(nLons, nLats, nAlts); + electron_temperature_scgc.ones(); + + tmp.sources_scgc.set_size(nLons, nLats, nAlts); + tmp.sources_scgc.zeros(); + tmp.losses_scgc.set_size(nLons, nLats, nAlts); + tmp.losses_scgc.zeros(); // This gets a bunch of the species-dependent characteristics: iErr = read_planet_file(input, report); - } // ----------------------------------------------------------------------------- @@ -101,13 +101,13 @@ int Ions::read_planet_file(Inputs input, Report report) { std::string hash; std::ifstream infile_ptr; - report.print(3,"In read_planet_file for Ions"); + report.print(3, "In read_planet_file for Ions"); infile_ptr.open(input.get_planet_species_file()); if (!infile_ptr.is_open()) { std::cout << "Could not open input file: " - << input.get_planet_species_file() << "!!!\n"; + << input.get_planet_species_file() << "!!!\n"; iErr = 1; } else { @@ -118,104 +118,61 @@ int Ions::read_planet_file(Inputs input, Report report) { hash = find_next_hash(infile_ptr); if (report.test_verbose(4)) - std::cout << "hash : -->" << hash << "<--\n"; + std::cout << "hash : -->" << hash << "<--\n"; if (hash == "#ions") { - // Read in the characteristics as CSVs: - report.print(4,"Found #ions!"); - - std::vector> lines = read_csv(infile_ptr); - - // I should totally redo the initialization of the species, - // since we could just do it here, but that is for the future. - - if (lines.size()-1 != nIons) { - std::cout << "number of ion species (nIons) defined in planet.h file : " - << nIons << "\n"; - std::cout << "number of ions defined in planet.in file : " - << lines.size() << "\n"; - std::cout << "These don't match!\n"; - iErr = 1; - } else { - - // assume order of rows right now: - // name, mass, charge, advect - - for (int iSpecies=0; iSpecies < nIons; iSpecies++) { - report.print(5, "setting ion species " + lines[iSpecies+1][0]); - species[iSpecies].cName = lines[iSpecies+1][0]; - species[iSpecies].mass = stof(lines[iSpecies+1][1])*amu; - species[iSpecies].charge = stoi(lines[iSpecies+1][2]); - species[iSpecies].DoAdvect = stoi(lines[iSpecies+1][3]); - } - - species[nIons].cName = "e-"; - species[nIons].mass = mass_electron; - species[nIons].charge = -1; - species[nIons].DoAdvect = 0; - } - + // Read in the characteristics as CSVs: + report.print(4, "Found #ions!"); + + std::vector> lines = read_csv(infile_ptr); + if (lines.size()-1 != nIons) { + std::cout << "num of ion species (nIons) defined in planet.h file : " + << nIons << "\n"; + std::cout << "number of ions defined in planet.in file : " + << lines.size() << "\n"; + std::cout << "These don't match!\n"; + iErr = 1; + } else { + // assume order of rows right now: + // name, mass, charge, advect + for (int iSpecies=0; iSpecies < nIons; iSpecies++) { + report.print(5, "setting ion species " + lines[iSpecies+1][0]); + species[iSpecies].cName = lines[iSpecies+1][0]; + species[iSpecies].mass = stof(lines[iSpecies+1][1])*amu; + species[iSpecies].charge = stoi(lines[iSpecies+1][2]); + species[iSpecies].DoAdvect = stoi(lines[iSpecies+1][3]); + } + species[nIons].cName = "e-"; + species[nIons].mass = mass_electron; + species[nIons].charge = -1; + species[nIons].DoAdvect = 0; + } } - if (infile_ptr.eof()) IsDone = 1; - } - infile_ptr.close(); - } - return iErr; - } // ----------------------------------------------------------------------------- -// +// // ----------------------------------------------------------------------------- -void Ions::fill_electrons(Grid grid, - Report &report) { +void Ions::fill_electrons(Report &report) { - long nLons, nLats, nAlts, iLon, iLat, iAlt, index; int iSpecies; - float electron_density; - + std::string function = "Ions::fill_electrons"; static int iFunction = -1; - report.enter(function, iFunction); - - if (grid.get_IsGeoGrid()) { - nLons = nGeoLonsG; - nLats = nGeoLatsG; - nAlts = nGeoAltsG; - } else { - nLons = nMagLonsG; - nLats = nMagLatsG; - nAlts = nMagAltsG; - } - - for (iLon = 0; iLon < nLons; iLon++) { - for (iLat = 0; iLat < nLats; iLat++) { - for (iAlt = 0; iAlt < nAlts; iAlt++) { + report.enter(function, iFunction); - if (grid.get_IsGeoGrid()) { - index = ijk_geo_s3gc(iLon,iLat,iAlt); - } else { - index = ijk_mag_s3gc(iLon,iLat,iAlt); - } - - electron_density = 0.0; - - for (iSpecies=0; iSpecies < nIons; iSpecies++) - electron_density = electron_density + species[iSpecies].density_s3gc[index]; - - species[nIons].density_s3gc[index] = electron_density; - density_s3gc[index] = electron_density; - - } - } - } + species[nIons].density_scgc.zeros(); + for (iSpecies=0; iSpecies < nIons; iSpecies++) + species[nIons].density_scgc = + species[nIons].density_scgc + species[iSpecies].density_scgc; + density_scgc = species[nIons].density_scgc; report.exit(function); return; diff --git a/src/main.cpp b/src/main.cpp index c745408e..0c366213 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,14 +1,12 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include - #include "../include/times.h" #include "../include/inputs.h" #include "../include/report.h" - #include "../include/neutrals.h" #include "../include/euv.h" #include "../include/grid.h" @@ -31,19 +29,19 @@ int main() { Indices indices(input); // Geo grid stuff: - Grid gGrid(nGeoLonsG, nGeoLatsG, nGeoAltsG); + Grid gGrid(nGeoLonsG, nGeoLatsG, nGeoAltsG, nGeoGhosts); gGrid.init_geo_grid(planet, input, report); gGrid.fill_grid(planet, report); // Magnetic grid stuff: - Grid mGrid(nMagLonsG, nMagLatsG, nMagAltsG); - + Grid mGrid(nMagLonsG, nMagLatsG, nMagAltsG, nMagGhosts); + Neutrals neutrals(gGrid, input, report); - Ions ions(input, report); - neutrals.pair_euv(euv, ions, report); + Ions ions(gGrid, input, report); + neutrals.pair_euv(euv, ions, report); Chemistry chemistry(neutrals, ions, input, report); - + // This is for the initial output. If it is not a restart, this will go: if (time.check_time_gate(input.get_dt_output(0))) { iErr = output(neutrals, ions, gGrid, time, planet, input, report); @@ -52,36 +50,32 @@ int main() { // This is advancing now... double dt_couple = 1800.0; - + // The way most codes are set up in the SWMF is that there are two // times, an end time which ends the simulation, and an intermediate // time, which allows coupling or something to happen. So, typically // the advance functions should only go to this intermediate time, // then a loop around that goes to the end time. Then, the code can // be made into a library and run externally. - + while (time.get_current() < time.get_end()) { time.increment_intermediate(dt_couple); while (time.get_current() < time.get_intermediate()) iErr = advance(planet, - gGrid, - time, - euv, - neutrals, - ions, - chemistry, - indices, - input, - report); + gGrid, + time, + euv, + neutrals, + ions, + chemistry, + indices, + input, + report); // Do some coupling here. But we have no coupling to do. Sad. - } - report.times(); - return iErr; - } diff --git a/src/neutrals.cpp b/src/neutrals.cpp index b0c22161..29e0c03c 100644 --- a/src/neutrals.cpp +++ b/src/neutrals.cpp @@ -1,9 +1,10 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include #include #include +#include #include "../include/constants.h" #include "../include/inputs.h" @@ -14,99 +15,104 @@ #include "../include/report.h" #include "../include/earth.h" +using namespace arma; + // ----------------------------------------------------------------------------- -// +// Create a single species by filling the species structure // ----------------------------------------------------------------------------- -Neutrals::species_chars Neutrals::create_species() { +Neutrals::species_chars Neutrals::create_species(Grid grid) { - long iDir, iLon, iLat, iAlt, index; + int64_t iDir, iLon, iLat, iAlt, index; species_chars tmp; - long iTotal = long(nGeoLonsG) * long(nGeoLatsG) * long(nGeoAltsG); + int64_t nLons = grid.get_nLons(); + int64_t nLats = grid.get_nLats(); + int64_t nAlts = grid.get_nAlts(); // Constants: tmp.DoAdvect = 0; tmp.thermal_cond = 0.0; tmp.thermal_exp = 0.0; tmp.lower_bc_density = -1.0; - - tmp.density_s3gc = (float*) malloc( iTotal * sizeof(float) ); - tmp.velocity_v3gc = (float*) malloc( long(3)*iTotal * sizeof(float) ); - tmp.chapman_s3gc = (float*) malloc( iTotal * sizeof(float) ); - tmp.ionization_s3gc = (float*) malloc( iTotal * sizeof(float) ); - - for (iLon = 0; iLon < nGeoLonsG; iLon++) { - for (iLat = 0; iLat < nGeoLatsG; iLat++) { - for (iAlt = 0; iAlt < nGeoAltsG; iAlt++) { - - index = ijk_geo_s3gc(iLon,iLat,iAlt); - - tmp.density_s3gc[index] = 1.0e-32; - tmp.chapman_s3gc[index] = 1.0e-32; - tmp.ionization_s3gc[index] = 1.0e-32; - - for (iDir = 0; iDir < 3; iDir++) { - index = ijkl_geo_v3gc(iLon,iLat,iAlt,iDir); - tmp.velocity_v3gc[index] = 0.0; - } - - } - } - } - + + tmp.density_scgc.set_size(nLons, nLats, nAlts); + tmp.chapman_scgc.set_size(nLons, nLats, nAlts); + tmp.scale_height_scgc.set_size(nLons, nLats, nAlts); + tmp.ionization_scgc.set_size(nLons, nLats, nAlts); + + tmp.density_scgc.ones(); + tmp.chapman_scgc.ones(); + tmp.scale_height_scgc.ones(); + tmp.ionization_scgc.zeros(); + + tmp.sources_scgc.set_size(nLons, nLats, nAlts); + tmp.sources_scgc.zeros(); + tmp.losses_scgc.set_size(nLons, nLats, nAlts); + tmp.losses_scgc.zeros(); + return tmp; - } // ----------------------------------------------------------------------------- -// +// Initialize neutrals // ----------------------------------------------------------------------------- Neutrals::Neutrals(Grid grid, Inputs input, Report report) { int iErr; species_chars tmp; - long iTotal = long(nGeoLonsG) * long(nGeoLatsG) * long(nGeoAltsG); + + int64_t nLons = grid.get_nLons(); + int64_t nLats = grid.get_nLats(); + int64_t nAlts = grid.get_nAlts(); report.print(2, "Initializing Neutrals"); for (int iSpecies=0; iSpecies < nSpecies; iSpecies++) { - tmp = create_species(); + tmp = create_species(grid); neutrals.push_back(tmp); } // State variables: - density_s3gc = (float*) malloc( iTotal * sizeof(float) ); - velocity_v3gc = (float*) malloc( long(3)*iTotal * sizeof(float) ); - temperature_s3gc = (float*) malloc( iTotal * sizeof(float) ); + + density_scgc.set_size(nLons, nLats, nAlts); + density_scgc.ones(); + temperature_scgc.set_size(nLons, nLats, nAlts); + temperature_scgc.ones(); // Derived quantities: - rho_s3gc = (float*) malloc( iTotal * sizeof(float) ); - mean_major_mass_s3gc = (float*) malloc( iTotal * sizeof(float) ); - pressure_s3gc = (float*) malloc( iTotal * sizeof(float) ); - sound_s3gc = (float*) malloc( iTotal * sizeof(float) ); + + rho_scgc.set_size(nLons, nLats, nAlts); + rho_scgc.ones(); + mean_major_mass_scgc.set_size(nLons, nLats, nAlts); + mean_major_mass_scgc.ones(); + pressure_scgc.set_size(nLons, nLats, nAlts); + pressure_scgc.ones(); + sound_scgc.set_size(nLons, nLats, nAlts); + sound_scgc.ones(); // Heating and cooling parameters: - Cv_s3gc = (float*) malloc( iTotal * sizeof(float) ); - gamma_s3gc = (float*) malloc( iTotal * sizeof(float) ); - kappa_s3gc = (float*) malloc( iTotal * sizeof(float) ); + Cv_scgc.set_size(nLons, nLats, nAlts); + Cv_scgc.ones(); + gamma_scgc.set_size(nLons, nLats, nAlts); + gamma_scgc.zeros(); + kappa_scgc.set_size(nLons, nLats, nAlts); + kappa_scgc.zeros(); - // Source Terms: - heating_euv_s3gc = (float*) malloc( iTotal * sizeof(float) ); - conduction_s3gc = (float*) malloc( iTotal * sizeof(float) ); + conduction_scgc.set_size(nLons, nLats, nAlts); + heating_euv_scgc.set_size(nLons, nLats, nAlts); heating_efficiency = input.get_euv_heating_eff_neutrals(); - + initial_temperatures = NULL; initial_altitudes = NULL; - + // This gets a bunch of the species-dependent characteristics: iErr = read_planet_file(input, report); // This specifies the initial conditions for the neutrals: iErr = initial_conditions(grid, input, report); - } // ----------------------------------------------------------------------------- @@ -119,13 +125,13 @@ int Neutrals::read_planet_file(Inputs input, Report report) { std::string hash; std::ifstream infile_ptr; - report.print(3,"In read_planet_file for Neutrals"); + report.print(3, "In read_planet_file for Neutrals"); infile_ptr.open(input.get_planet_species_file()); if (!infile_ptr.is_open()) { std::cout << "Could not open input file: " - << input.get_planet_species_file() << "!!!\n"; + << input.get_planet_species_file() << "!!!\n"; iErr = 1; } else { @@ -134,176 +140,179 @@ int Neutrals::read_planet_file(Inputs input, Report report) { while (!IsDone) { hash = find_next_hash(infile_ptr); - if (report.test_verbose(4)) - std::cout << "hash : -->" << hash << "<--\n"; + std::cout << "hash : -->" << hash << "<--\n"; if (hash == "#neutrals") { - - // Read in the characteristics as CSVs: - report.print(4,"Found #neutrals!"); - - std::vector> lines = read_csv(infile_ptr); - - // I should totally redo the initialization of the species, - // since we could just do it here, but that is for the future. - - if (lines.size()-1 != nSpecies) { - std::cout << "number of neutrals species defined in planet.h file : " - << nSpecies << "\n"; - std::cout << "number of species defined in planet.in file : " - << lines.size() << "\n"; - std::cout << "These don't match!\n"; - iErr = 1; - } else { - - // assume order of rows right now: - // name, mass, vibration, thermal_cond, thermal_exp, advect, lower BC - - for (int iSpecies=0; iSpecies < nSpecies; iSpecies++) { - report.print(5, "setting neutral species " + lines[iSpecies+1][0]); - neutrals[iSpecies].cName = lines[iSpecies+1][0]; - neutrals[iSpecies].mass = stof(lines[iSpecies+1][1])*amu; - neutrals[iSpecies].vibe = stof(lines[iSpecies+1][2]); - neutrals[iSpecies].thermal_cond = stof(lines[iSpecies+1][3]); - neutrals[iSpecies].thermal_exp = stof(lines[iSpecies+1][4]); - neutrals[iSpecies].DoAdvect = stoi(lines[iSpecies+1][5]); - neutrals[iSpecies].lower_bc_density = stof(lines[iSpecies+1][6]); - } - - } - - } + // Read in the characteristics as CSVs: + report.print(4, "Found #neutrals!"); + + std::vector> lines = read_csv(infile_ptr); + + // I should totally redo the initialization of the species, + // since we could just do it here, but that is for the future. + + if (lines.size()-1 != nSpecies) { + std::cout << "number of neutrals species defined in planet.h file : " + << nSpecies << "\n"; + std::cout << "number of species defined in planet.in file : " + << lines.size() << "\n"; + std::cout << "These don't match!\n"; + iErr = 1; + } else { + // assume order of rows right now: + // name, mass, vibration, thermal_cond, thermal_exp, advect, lower BC + + for (int iSpecies=0; iSpecies < nSpecies; iSpecies++) { + report.print(5, "setting neutral species " + lines[iSpecies+1][0]); + neutrals[iSpecies].cName = lines[iSpecies+1][0]; + neutrals[iSpecies].mass = stof(lines[iSpecies+1][1])*amu; + neutrals[iSpecies].vibe = stof(lines[iSpecies+1][2]); + neutrals[iSpecies].thermal_cond = stof(lines[iSpecies+1][3]); + neutrals[iSpecies].thermal_exp = stof(lines[iSpecies+1][4]); + neutrals[iSpecies].DoAdvect = stoi(lines[iSpecies+1][5]); + neutrals[iSpecies].lower_bc_density = stof(lines[iSpecies+1][6]); + } // iSpecies + } // else size + } // #neutrals if (hash == "#temperature") { - report.print(4,"Found #temperatures!"); - - std::vector> temps = read_csv(infile_ptr); - - int nTemps = temps.size()-1; - initial_temperatures = (float*) malloc(nTemps * sizeof(float) ); - initial_altitudes = (float*) malloc(nTemps * sizeof(float) ); - for (int iTemp=0; iTemp < nTemps; iTemp++) { - report.print(5, "reading initial temp alt " + temps[iTemp+1][0]); - // convert altitudes from km to m - initial_altitudes[iTemp] = stof(temps[iTemp+1][0]) * 1000; - initial_temperatures[iTemp] = stof(temps[iTemp+1][1]); - } - nInitial_temps = nTemps; - } - - if (infile_ptr.eof()) IsDone = 1; - - } - - infile_ptr.close(); - - } + report.print(4, "Found #temperatures!"); - return iErr; - -} + std::vector> temps = read_csv(infile_ptr); -// ----------------------------------------------------------------------------- -// -// ----------------------------------------------------------------------------- + int nTemps = temps.size()-1; + initial_temperatures = (float*) malloc(nTemps * sizeof(float) ); + initial_altitudes = (float*) malloc(nTemps * sizeof(float) ); + for (int iTemp=0; iTemp < nTemps; iTemp++) { + report.print(5, "reading initial temp alt " + temps[iTemp+1][0]); + // convert altitudes from km to m + initial_altitudes[iTemp] = stof(temps[iTemp+1][0]) * 1000; + initial_temperatures[iTemp] = stof(temps[iTemp+1][1]); + } // for iTemp -float Neutrals::calc_scale_height(int iSpecies, - long index, - Grid grid) { + nInitial_temps = nTemps; + } // #temperature - float g = grid.gravity_s3gc[index]; - float t = temperature_s3gc[index]; - float m = neutrals[iSpecies].mass; - float H = boltzmanns_constant * t / m / g; + if (infile_ptr.eof()) IsDone = 1; + } // while !IsDone - return H; + infile_ptr.close(); + } // else (isopen) + return iErr; } // ----------------------------------------------------------------------------- -// +// This is a place holder for creating initial conditions // ----------------------------------------------------------------------------- int Neutrals::initial_conditions(Grid grid, Inputs input, Report report) { int iErr = 0; - long iDir, iLon, iLat, iAlt, index, iA, indexm; + int64_t iDir, iLon, iLat, iAlt, index, iA, indexm; float alt, r, H; report.print(3, "Creating Neutrals initial_condition"); - + // --------------------------------------------------------------------- // This section assumes we want a hydrostatic solution given the // temperature profile in the planet.in file. // --------------------------------------------------------------------- - if (nInitial_temps > 0) { + int64_t nLons = grid.get_nLons(); + int64_t nLats = grid.get_nLats(); + int64_t nAlts = grid.get_nAlts(); + int64_t nGCs = grid.get_nGCs(); + + // Let's assume that the altitudes are not dependent on lat/lon: + + fvec alt1d(nAlts); + fvec temp1d(nAlts); + + fmat H2d(nLons, nLats); - for (iLon = 0; iLon < nGeoLonsG; iLon++) { - for (iLat = 0; iLat < nGeoLatsG; iLat++) { - for (iAlt = 0; iAlt < nGeoAltsG; iAlt++) { - - index = ijk_geo_s3gc(iLon,iLat,iAlt); - - alt = grid.geoAlt_s3gc[index]; - - // Need to make a generic linear interpolator!!! - - // Find temperatures: - if (alt <= initial_altitudes[0]) { - temperature_s3gc[index] = initial_temperatures[0]; - } else { - if (alt >= initial_altitudes[nInitial_temps-1]) { - temperature_s3gc[index] = initial_temperatures[nInitial_temps-1]; - } else { - // Linear interpolation! - iA = 0; - while (alt > initial_altitudes[iA]) iA++; - iA--; - // alt will be between iA and iA+1: - r = (alt - initial_altitudes[iA]) / - (initial_altitudes[iA+1] - initial_altitudes[iA]); - temperature_s3gc[index] = - (1.0-r) * initial_temperatures[iA ] + - ( r) * initial_temperatures[iA+1]; - } - } - - // Now do the neutrals. - - // For the bottom, set to the constant conditions: - - if (iAlt == 0) { - for (int iSpecies=0; iSpecies < nSpecies; iSpecies++) - neutrals[iSpecies].density_s3gc[index] = - neutrals[iSpecies].lower_bc_density; - } else { - - // Let's use the cell below to set the density. - // Calculate scale height and then use hydrostatic balance - // to derive the density: - - indexm = ijk_geo_s3gc(iLon,iLat,iAlt-1); - - for (int iSpecies=0; iSpecies < nSpecies; iSpecies++) { - H = calc_scale_height(iSpecies, index, grid); - - neutrals[iSpecies].density_s3gc[index] = - neutrals[iSpecies].density_s3gc[indexm] * - exp(-grid.dalt_lower_s3gc[index]/H); - - } - - } - } + alt1d = grid.geoAlt_scgc.tube(0, 0); + + if (nInitial_temps > 0) { + for (iAlt = 0; iAlt < nAlts; iAlt++) { + alt = alt1d(iAlt); + // Find temperatures: + if (alt <= initial_altitudes[0]) { + temp1d[iAlt] = initial_temperatures[0]; + } else { + if (alt >= initial_altitudes[nInitial_temps-1]) { + temp1d[iAlt] = initial_temperatures[nInitial_temps-1]; + } else { + // Linear interpolation! + iA = 0; + while (alt > initial_altitudes[iA]) iA++; + iA--; + // alt will be between iA and iA+1: + r = (alt - initial_altitudes[iA]) / + (initial_altitudes[iA+1] - initial_altitudes[iA]); + temp1d[iAlt] = + (1.0-r) * initial_temperatures[iA] + + (r) * initial_temperatures[iA+1]; + } } } + } else { + temp1d = 200.0; + } + + // spread the 1D temperature across the globe: + for (iLon = 0; iLon < nLons; iLon++) { + for (iLat = 0; iLat < nLats; iLat++) { + temperature_scgc.tube(iLon, iLat) = temp1d; + } + } + + // Set the lower boundary condition: + for (int iSpecies=0; iSpecies < nSpecies; iSpecies++) { + + neutrals[iSpecies]. + density_scgc.slice(0). + fill(neutrals[iSpecies].lower_bc_density); + + // Integrate with hydrostatic equilibrium up: + for (iAlt = 1; iAlt < nAlts; iAlt++) { + neutrals[iSpecies].scale_height_scgc.slice(iAlt) = + boltzmanns_constant * + temperature_scgc.slice(iAlt) / + (neutrals[iSpecies].mass * + grid.gravity_scgc.slice(iAlt)); + neutrals[iSpecies].density_scgc.slice(iAlt) = + neutrals[iSpecies].density_scgc.slice(iAlt-1) % + exp(-grid.dalt_lower_scgc.slice(iAlt) / + neutrals[iSpecies].scale_height_scgc.slice(iAlt)); + } } + calc_mass_density(report); return iErr; - +} + +//---------------------------------------------------------------------- +// set_bcs +//---------------------------------------------------------------------- + +void Neutrals::set_bcs(Report &report) { + + std::string function = "Neutrals::set_bcs"; + static int iFunction = -1; + report.enter(function, iFunction); + + int64_t nAlts = temperature_scgc.n_slices; + + // temperature_scgc.slice(nAlts-2).fill(800.0); + // temperature_scgc.slice(nAlts-1).fill(800.0); + + temperature_scgc.slice(nAlts-2) = temperature_scgc.slice(nAlts-3); + temperature_scgc.slice(nAlts-1) = temperature_scgc.slice(nAlts-2); + + report.exit(function); } //---------------------------------------------------------------------- @@ -320,54 +329,48 @@ int Neutrals::pair_euv(Euv euv, Ions ions, Report report) { int iErr = 0; - if (report.test_verbose(3)) - std::cout << "Pairing EUV abs/ion cross sections... (Neutrals::pair_euv) \n"; - - for (int iSpecies=0; iSpecies < nSpecies; iSpecies++) { + if (report.test_verbose(3)) std::cout << "Neutrals::pair_euv \n"; + + for (int iSpecies = 0; iSpecies < nSpecies; iSpecies++) { if (report.test_verbose(5)) std::cout << neutrals[iSpecies].cName << "\n"; neutrals[iSpecies].iEuvAbsId_ = -1; neutrals[iSpecies].nEuvIonSpecies = 0; - + // Check each row to see if the first column "name" matches: for (int iEuv=0; iEuv < euv.waveinfo.size(); iEuv++) { if (report.test_verbose(6)) - std::cout << " " << euv.waveinfo[iEuv].name << "\n"; + std::cout << " " << euv.waveinfo[iEuv].name << "\n"; // if this matches... if (neutrals[iSpecies].cName == euv.waveinfo[iEuv].name) { - // First see if we can find absorbtion: - if (euv.waveinfo[iEuv].type == "abs") { - if (report.test_verbose(5)) std::cout << " Found absorbtion\n"; - neutrals[iSpecies].iEuvAbsId_ = iEuv; - } - - // Next see if we can find ionizations: - if (euv.waveinfo[iEuv].type == "ion") { - - // Loop through the ions to see if names match: - for (int iIon=0; iIon " - << ions.species[iIon].cName << "\n"; - neutrals[iSpecies].iEuvIonId_.push_back(iEuv); - neutrals[iSpecies].iEuvIonSpecies_.push_back(iIon); - neutrals[iSpecies].nEuvIonSpecies++; - } - } - } - } - - } - - } - + // First see if we can find absorbtion: + if (euv.waveinfo[iEuv].type == "abs") { + if (report.test_verbose(5)) std::cout << " Found absorbtion\n"; + neutrals[iSpecies].iEuvAbsId_ = iEuv; + } + + // Next see if we can find ionizations: + if (euv.waveinfo[iEuv].type == "ion") { + + // Loop through the ions to see if names match: + for (int iIon = 0; iIon < nIons; iIon++) { + if (ions.species[iIon].cName == euv.waveinfo[iEuv].to) { + if (report.test_verbose(5)) + std::cout << " Found ionization!! --> " + << ions.species[iIon].cName << "\n"; + neutrals[iSpecies].iEuvIonId_.push_back(iEuv); + neutrals[iSpecies].iEuvIonSpecies_.push_back(iIon); + neutrals[iSpecies].nEuvIonSpecies++; + } // if to + } // iIon loop + } // if ionization + } // if species is name + } // for iEuv + } // for iSpecies return iErr; - } - diff --git a/src/output.cpp b/src/output.cpp index 415abe1a..333bd6a1 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -1,4 +1,4 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include @@ -15,22 +15,53 @@ using namespace netCDF; using namespace netCDF::exceptions; +void output_variable_3d(std::vector count_start, + std::vector count_end, + fcube value, + NcVar variable) { + + int64_t nX = value.n_rows; + int64_t nY = value.n_cols; + int64_t nZ = value.n_slices; + int64_t iX, iY, iZ, iTotal, index; + + iTotal = nX * nY * nZ; + + float *tmp_s3gc = (float*) malloc(iTotal * sizeof(float)); + + for (iX = 0; iX < nX; iX++) { + for (iY = 0; iY < nY; iY++) { + for (iZ = 0; iZ < nZ; iZ++) { + index = iX*nY*nZ + iY*nZ + iZ; + tmp_s3gc[index] = value(iX, iY, iZ); + } + } + } + + variable.putVar(count_start, count_end, tmp_s3gc); + free(tmp_s3gc); +} + int output(Neutrals neutrals, - Ions ions, - Grid grid, - Times time, - Planets planet, - Inputs args, - Report &report) { + Ions ions, + Grid grid, + Times time, + Planets planet, + Inputs args, + Report &report) { int iErr = 0; int nOutputs = args.get_n_outputs(); int IsGeoGrid = grid.get_IsGeoGrid(); - std::string function="output"; + int64_t nLons = grid.get_nLons(); + int64_t nLats = grid.get_nLats(); + int64_t nAlts = grid.get_nAlts(); + + std::string function = "output"; static int iFunction = -1; - report.enter(function, iFunction); + report.enter(function, iFunction); for (int iOutput = 0; iOutput < nOutputs; iOutput++) { @@ -50,22 +81,20 @@ int output(Neutrals neutrals, time_string = time.get_YMD_HMS(); file_name = file_pre + "_" + time_string + file_ext; - + // Create the file: NcFile ncdf_file(file_name, NcFile::replace); // Add dimensions: - NcDim lonDim = ncdf_file.addDim("Longitude", nGeoLonsG); - NcDim latDim = ncdf_file.addDim("Latitude", nGeoLatsG); - NcDim altDim = ncdf_file.addDim("Altitude", nGeoAltsG); + NcDim lonDim = ncdf_file.addDim("Longitude", nLons); + NcDim latDim = ncdf_file.addDim("Latitude", nLats); + NcDim altDim = ncdf_file.addDim("Altitude", nAlts); - // If we wanted 1D variables, we would do something like this, but - // since all of out variables will be 3d, skip this: // Define the Coordinate Variables - //NcVar altVar = ncdf_file.addVar("Altitude", ncFloat, altDim); - + // Define the netCDF variables for the 3D data. // First create a vector of dimensions: + std::vector dimVector; dimVector.push_back(lonDim); dimVector.push_back(latDim); @@ -75,47 +104,51 @@ int output(Neutrals neutrals, NcVar latVar = ncdf_file.addVar("Latitude", ncFloat, dimVector); NcVar altVar = ncdf_file.addVar("Altitude", ncFloat, dimVector); - lonVar.putAtt(UNITS,"radians"); - latVar.putAtt(UNITS,"radians"); - altVar.putAtt(UNITS,"meters"); + lonVar.putAtt(UNITS, "radians"); + latVar.putAtt(UNITS, "radians"); + altVar.putAtt(UNITS, "meters"); - std::vector startp,countp; + std::vector startp, countp; startp.push_back(0); startp.push_back(0); startp.push_back(0); - countp.push_back(nGeoLonsG); - countp.push_back(nGeoLatsG); - countp.push_back(nGeoAltsG); + countp.push_back(nLons); + countp.push_back(nLats); + countp.push_back(nAlts); // Output longitude, latitude, altitude 3D arrays: - lonVar.putVar(startp, countp, grid.geoLon_s3gc); - latVar.putVar(startp, countp, grid.geoLat_s3gc); - altVar.putVar(startp, countp, grid.geoAlt_s3gc); + + output_variable_3d(startp, countp, grid.geoLon_scgc, lonVar); + output_variable_3d(startp, countp, grid.geoLat_scgc, latVar); + output_variable_3d(startp, countp, grid.geoAlt_scgc, altVar); // ---------------------------------------------- // Neutral Densities and Temperature // ---------------------------------------------- if (type_output == "neutrals" || - type_output == "states") { - - // Output all species densities: - std::vector denVar; - for (int iSpecies=0; iSpecies < nSpecies; iSpecies++) { - if (report.test_verbose(3)) - std::cout << "Outputting Var : " - << neutrals.neutrals[iSpecies].cName << "\n"; - denVar.push_back(ncdf_file.addVar(neutrals.neutrals[iSpecies].cName, ncFloat, dimVector)); - denVar[iSpecies].putAtt(UNITS,neutrals.density_unit); - denVar[iSpecies].putVar(startp, countp, neutrals.neutrals[iSpecies].density_s3gc); - } - - // Output bulk temperature: - NcVar tempVar = ncdf_file.addVar(neutrals.temperature_name, ncFloat, dimVector); - tempVar.putAtt(UNITS,neutrals.temperature_unit); - tempVar.putVar(startp, countp, neutrals.temperature_s3gc); - + type_output == "states") { + + // Output all species densities: + std::vector denVar; + for (int iSpecies=0; iSpecies < nSpecies; iSpecies++) { + if (report.test_verbose(3)) + std::cout << "Outputting Var : " + << neutrals.neutrals[iSpecies].cName << "\n"; + denVar.push_back(ncdf_file.addVar(neutrals.neutrals[iSpecies].cName, + ncFloat, dimVector)); + denVar[iSpecies].putAtt(UNITS, neutrals.density_unit); + output_variable_3d(startp, countp, + neutrals.neutrals[iSpecies].density_scgc, + denVar[iSpecies]); + } + + // Output bulk temperature: + NcVar tempVar = ncdf_file.addVar(neutrals.temperature_name, + ncFloat, dimVector); + tempVar.putAtt(UNITS, neutrals.temperature_unit); + output_variable_3d(startp, countp, neutrals.temperature_scgc, tempVar); } // ---------------------------------------------- @@ -123,28 +156,25 @@ int output(Neutrals neutrals, // ---------------------------------------------- if (type_output == "ions" || - type_output == "states") { - - // Output all species densities: - std::vector ionVar; - for (int iSpecies=0; iSpecies < nIons; iSpecies++) { - if (report.test_verbose(3)) - std::cout << "Outputting Var : " - << ions.species[iSpecies].cName << "\n"; - ionVar.push_back(ncdf_file.addVar(ions.species[iSpecies].cName, ncFloat, dimVector)); - ionVar[iSpecies].putAtt(UNITS,neutrals.density_unit); - ionVar[iSpecies].putVar(startp, countp, ions.species[iSpecies].density_s3gc); - } - - ionVar.push_back(ncdf_file.addVar("e-", ncFloat, dimVector)); - ionVar[nIons].putAtt(UNITS,neutrals.density_unit); - ionVar[nIons].putVar(startp, countp, ions.density_s3gc); - - // // Output bulk temperature: - // NcVar tempVar = ncdf_file.addVar(neutrals.temperature_name, ncFloat, dimVector); - // tempVar.putAtt(UNITS,neutrals.temperature_unit); - // tempVar.putVar(startp, countp, neutrals.temperature_s3gc); - + type_output == "states") { + + // Output all species densities: + std::vector ionVar; + for (int iSpecies=0; iSpecies < nIons; iSpecies++) { + if (report.test_verbose(3)) + std::cout << "Outputting Var : " + << ions.species[iSpecies].cName << "\n"; + ionVar.push_back(ncdf_file.addVar(ions.species[iSpecies].cName, + ncFloat, dimVector)); + ionVar[iSpecies].putAtt(UNITS, neutrals.density_unit); + output_variable_3d(startp, countp, + ions.species[iSpecies].density_scgc, + ionVar[iSpecies]); + } + + ionVar.push_back(ncdf_file.addVar("e-", ncFloat, dimVector)); + ionVar[nIons].putAtt(UNITS, neutrals.density_unit); + output_variable_3d(startp, countp, ions.density_scgc, ionVar[nIons]); } // ---------------------------------------------- @@ -152,41 +182,35 @@ int output(Neutrals neutrals, // ---------------------------------------------- if (type_output == "bfield") { - NcVar mLatVar = ncdf_file.addVar("Magnetic Latitude", ncFloat, dimVector); - mLatVar.putAtt(UNITS,"radians"); - mLatVar.putVar(startp, countp, grid.magLat_s3gc); - NcVar mLonVar = ncdf_file.addVar("Magnetic Longitude", ncFloat, dimVector); - mLonVar.putAtt(UNITS,"radians"); - mLonVar.putVar(startp, countp, grid.magLon_s3gc); - - // Output magnetic field components: - float *bfield_component_s3gc; - long nPointsTotal = grid.get_nPointsInGrid(); - bfield_component_s3gc = (float*) malloc( nPointsTotal * sizeof(float) ); - - NcVar bxVar = ncdf_file.addVar("Bx", ncFloat, dimVector); - bxVar.putAtt(UNITS,"nT"); - get_vector_component(grid.bfield_v3gc, 0, IsGeoGrid, bfield_component_s3gc); - bxVar.putVar(startp, countp, bfield_component_s3gc); - - NcVar byVar = ncdf_file.addVar("By", ncFloat, dimVector); - byVar.putAtt(UNITS,"nT"); - get_vector_component(grid.bfield_v3gc, 1, IsGeoGrid, bfield_component_s3gc); - byVar.putVar(startp, countp, bfield_component_s3gc); - - NcVar bzVar = ncdf_file.addVar("Bz", ncFloat, dimVector); - bzVar.putAtt(UNITS,"nT"); - get_vector_component(grid.bfield_v3gc, 2, IsGeoGrid, bfield_component_s3gc); - bzVar.putVar(startp, countp, bfield_component_s3gc); - - } - + NcVar mLatVar = ncdf_file.addVar("Magnetic Latitude", + ncFloat, dimVector); + mLatVar.putAtt(UNITS, "radians"); + output_variable_3d(startp, countp, grid.magLat_scgc, mLatVar); + + NcVar mLonVar = ncdf_file.addVar("Magnetic Longitude", + ncFloat, dimVector); + mLonVar.putAtt(UNITS, "radians"); + output_variable_3d(startp, countp, grid.magLat_scgc, mLonVar); + + // Output magnetic field components: + + NcVar bxVar = ncdf_file.addVar("Bx", ncFloat, dimVector); + bxVar.putAtt(UNITS, "nT"); + output_variable_3d(startp, countp, grid.bfield_vcgc[0], bxVar); + + NcVar byVar = ncdf_file.addVar("By", ncFloat, dimVector); + byVar.putAtt(UNITS, "nT"); + output_variable_3d(startp, countp, grid.bfield_vcgc[1], bxVar); + + NcVar bzVar = ncdf_file.addVar("Bz", ncFloat, dimVector); + bzVar.putAtt(UNITS, "nT"); + output_variable_3d(startp, countp, grid.bfield_vcgc[2], bxVar); + } // if befield + ncdf_file.close(); + } // if time check + } // for iOutput - } - } - report.exit(function); return iErr; - } diff --git a/src/planets.cpp b/src/planets.cpp index 54e8fafe..c4ff75c0 100644 --- a/src/planets.cpp +++ b/src/planets.cpp @@ -1,9 +1,10 @@ +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) +// Full license can be found in License.md #include #include #include #include -#include // std::stringstream #include #include "../include/constants.h" @@ -104,10 +105,8 @@ float Planets::get_mu() { // ----------------------------------------------------------------------------- float Planets::get_star_to_planet_dist(Times time) { - int iErr = update(time); return planet.star_planet_distance; - } // ----------------------------------------------------------------------------- @@ -115,10 +114,8 @@ float Planets::get_star_to_planet_dist(Times time) { // ----------------------------------------------------------------------------- float Planets::get_orbit_angle(Times time) { - int iErr = update(time); return planet.orbit_angle; - } // ----------------------------------------------------------------------------- @@ -126,10 +123,8 @@ float Planets::get_orbit_angle(Times time) { // ----------------------------------------------------------------------------- float Planets::get_declination(Times time) { - int iErr = update(time); return planet.declination;; - } // ----------------------------------------------------------------------------- @@ -139,7 +134,7 @@ float Planets::get_declination(Times time) { int Planets::update(Times time) { int iErr = 0; - + // Update planetary stuff once a minute: if (time.get_current() - planet.update_time > 60.0) { @@ -186,10 +181,10 @@ int Planets::update(Times time) { float y_heliocentric = sma*sqrt(1-ecc*ecc)*sin(ecc_anomaly*dtor); float z_heliocentric = 0.0; - float true_anomaly = atan2(y_heliocentric,x_heliocentric)*rtod; + float true_anomaly = atan2(y_heliocentric, x_heliocentric)*rtod; planet.star_planet_distance = sqrt(x_heliocentric * x_heliocentric + - y_heliocentric * y_heliocentric); + y_heliocentric * y_heliocentric); // convert to J2000 coordinates with x-axis aligned with vernal equinox so // we can get solar longitude in the coorect system. We don't need z. @@ -237,12 +232,10 @@ int Planets::update(Times time) { double left_over = (rotations - int(rotations)) * 360.0; // put into radians, so it is consistent with the rest of the code: planet.longitude_offset = - fmod(planet.longitude_jb2000 + left_over + 180.0,360.0) * dtor; - + fmod(planet.longitude_jb2000 + left_over + 180.0, 360.0) * dtor; } return iErr; - } int Planets::set_planet(Inputs args, Report report) { @@ -250,12 +243,13 @@ int Planets::set_planet(Inputs args, Report report) { int iErr = 0; int IsFound = 0; - for (int i=0; i 0.01) { - std::cout << "Hmmm.... For some reason, the calculated length of year is different\n"; - std::cout << "Read in : " << planet.length_of_year/seconds_per_day << " (days) \n"; - std::cout << "Calculated from SMA : " << loy/seconds_per_day << " (days)\n"; - std::cout << "Trusting SMA version!\n"; - iErr = 1; + std::cout << "Hmmm.... For some reason, the calculated length "; + std::cout << "of year is different\n"; + std::cout << "Read in : " + << planet.length_of_year/seconds_per_day + << " (days) \n"; + std::cout << "Calculated from SMA : " + << loy/seconds_per_day << " (days)\n"; + std::cout << "Trusting SMA version!\n"; + iErr = 1; } planet.length_of_year = loy; planet.length_of_day = planets[i].length_of_day; planet.longitude_jb2000 = planets[i].longitude_jb2000; float rotrate = 2*pi/planet.length_of_day*(1.0+1.0/(loy/seconds_per_day)); - planet.omega = rotrate; // frequency (rad/s) - planet.rotation_period = 2*pi/rotrate; // (seconds) + planet.omega = rotrate; // frequency (rad/s) + planet.rotation_period = 2*pi/rotrate; // (seconds) planet.mass = planets[i].mass; planet.mu = planets[i].mass * gravitational_constant; - planet.equator_radius = planets[i].equator_radius * 1000.0; // km -> m - planet.polar_radius = planets[i].polar_radius * 1000.0; // km -> m - // Looking at Earth and Saturn, it seems like the Volumetric mean radius - // is at roughly 47 deg (cos(47)=0.68) - // Obviously an approximation... + planet.equator_radius = planets[i].equator_radius * 1000.0; // km -> m + planet.polar_radius = planets[i].polar_radius * 1000.0; // km -> m + // Looking at Earth and Saturn, it seems like the Volumetric + // mean radius is at roughly 47 deg (cos(47)=0.68) Obviously an + // approximation... planet.radius = 0.68 * planet.equator_radius + 0.32 * planet.polar_radius; if (report.test_verbose(2)) - std::cout << "Planet Radius set to : " - << planet.radius/1000.0 << " (km)\n"; + std::cout << "Planet Radius set to : " + << planet.radius/1000.0 << " (km)\n"; planet.dipole_strength = planets[i].dipole_strength; planet.dipole_rotation = planets[i].dipole_rotation; planet.dipole_tilt = planets[i].dipole_tilt; - for (int j=0; j<3; j++) - planet.dipole_center[j] = planets[i].dipole_center[j]; - + for (int j = 0; j < 3; j++) + planet.dipole_center[j] = planets[i].dipole_center[j]; + planet.update_time = -1e32; break; @@ -315,12 +313,11 @@ int Planets::set_planet(Inputs args, Report report) { if (!IsFound) { std::cout << "Can't file planet " << args.get_planet() - << " in planet file information!\n"; + << " in planet file information!\n"; iErr = 1; } return iErr; - } int Planets::read_file(Inputs args, Report report) { @@ -336,7 +333,7 @@ int Planets::read_file(Inputs args, Report report) { if (!myFile.is_open()) { std::cout << "Could not open planetary file : " - << args.get_planetary_file() << "\n"; + << args.get_planetary_file() << "\n"; iErr = 1; } else { @@ -347,57 +344,48 @@ int Planets::read_file(Inputs args, Report report) { int nLines = csv.size(); if (nLines <= 2) { - iErr = 1; + iErr = 1; } else { - for (int iLine = 2; iLine < nLines; iLine++) { - - // Some final rows can have comments in them, so we want to - // skip anything where the length of the string in column 2 - // is == 0: - - if (csv[iLine][1].length() > 0) { - - tmp.name = make_lower(csv[iLine][0]); - tmp.semimajoraxis = stof(csv[iLine][1]); - tmp.eccentricity = stof(csv[iLine][2]); - tmp.inclination = stof(csv[iLine][3]); - tmp.meanlongitude = stof(csv[iLine][4]); - tmp.perihelionlongitude = stof(csv[iLine][5]); - tmp.nodelongitude = stof(csv[iLine][6]); - tmp.rates_semimajoraxis = stof(csv[iLine][7]); - tmp.rates_eccentricity = stof(csv[iLine][8]); - tmp.rates_inclination = stof(csv[iLine][9]); - tmp.rates_meanlongitude = stof(csv[iLine][10]); - tmp.rates_perihelionlongitude = stof(csv[iLine][11]); - tmp.rates_nodelongitude = stof(csv[iLine][12]); - tmp.length_of_day = stof(csv[iLine][13])*seconds_per_hour; - tmp.length_of_year = stof(csv[iLine][14])*seconds_per_day; - tmp.longitude_jb2000 = stof(csv[iLine][15]); - tmp.mass = stof(csv[iLine][16]); - tmp.equator_radius = stof(csv[iLine][17]); - tmp.polar_radius = stof(csv[iLine][18]); - tmp.planet_tilt = stof(csv[iLine][19])*dtor; - tmp.dipole_strength = stof(csv[iLine][20]); - tmp.dipole_rotation = stof(csv[iLine][21])*dtor; - tmp.dipole_tilt = stof(csv[iLine][22])*dtor; - for (int j=0; j<3; j++) - tmp.dipole_center[j] = stof(csv[iLine][23+j]); - - planets.push_back(tmp); - - } - - } - - } - - } - + for (int iLine = 2; iLine < nLines; iLine++) { + + // Some final rows can have comments in them, so we want to + // skip anything where the length of the string in column 2 + // is == 0: + + if (csv[iLine][1].length() > 0) { + tmp.name = make_lower(csv[iLine][0]); + tmp.semimajoraxis = stof(csv[iLine][1]); + tmp.eccentricity = stof(csv[iLine][2]); + tmp.inclination = stof(csv[iLine][3]); + tmp.meanlongitude = stof(csv[iLine][4]); + tmp.perihelionlongitude = stof(csv[iLine][5]); + tmp.nodelongitude = stof(csv[iLine][6]); + tmp.rates_semimajoraxis = stof(csv[iLine][7]); + tmp.rates_eccentricity = stof(csv[iLine][8]); + tmp.rates_inclination = stof(csv[iLine][9]); + tmp.rates_meanlongitude = stof(csv[iLine][10]); + tmp.rates_perihelionlongitude = stof(csv[iLine][11]); + tmp.rates_nodelongitude = stof(csv[iLine][12]); + tmp.length_of_day = stof(csv[iLine][13])*seconds_per_hour; + tmp.length_of_year = stof(csv[iLine][14])*seconds_per_day; + tmp.longitude_jb2000 = stof(csv[iLine][15]); + tmp.mass = stof(csv[iLine][16]); + tmp.equator_radius = stof(csv[iLine][17]); + tmp.polar_radius = stof(csv[iLine][18]); + tmp.planet_tilt = stof(csv[iLine][19])*dtor; + tmp.dipole_strength = stof(csv[iLine][20]); + tmp.dipole_rotation = stof(csv[iLine][21])*dtor; + tmp.dipole_tilt = stof(csv[iLine][22])*dtor; + for (int j = 0; j < 3; j++) + tmp.dipole_center[j] = stof(csv[iLine][23+j]); + + planets.push_back(tmp); + } // if length + } // for iLine + } // else nLines + } // if good file myFile.close(); - - } - + } // else open file return iErr; - } diff --git a/src/read_f107_file.cpp b/src/read_f107_file.cpp index 16266b1f..f3f4308b 100644 --- a/src/read_f107_file.cpp +++ b/src/read_f107_file.cpp @@ -1,4 +1,4 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include @@ -10,8 +10,8 @@ #include "../include/time_conversion.h" int read_f107_file(std::string f107_file, - std::vector &time, - std::vector &f107) { + std::vector &time, + std::vector &f107) { std::ifstream myFile; int iErr; @@ -26,67 +26,59 @@ int read_f107_file(std::string f107_file, } else { int IsFound = 0; - int IsAdjusted = 0; + int IsAdjusted = 0; // At some point need to take adjusted into account. std::string line; - while (getline(myFile,line)) { - + while (getline(myFile, line)) { if (line == "#Element: adjusted") IsAdjusted = 1; if (line == "#yyyy-MM-dd HH:mm value qualifier description") { - IsFound = 1; - break; + IsFound = 1; + break; } - } // This means we found the line right before the data starts! Read // in the data! - + if (IsFound) { std::string tmp; - std::vector itime(7,0); - - while (getline(myFile,line)) { - - std::stringstream ss(line); - // year - getline(ss, tmp, '-'); - itime[0] = stoi(tmp); - // month - getline(ss, tmp, '-'); - itime[1] = stoi(tmp); - // day - getline(ss, tmp, ' '); - itime[2] = stoi(tmp); - // hour - getline(ss, tmp, ':'); - itime[3] = stoi(tmp); - // minute - getline(ss, tmp, ' '); - itime[4] = stoi(tmp); - itime[5] = 0; - itime[6] = 0; - time.push_back(time_int_to_real(itime)); - - // f107 - getline(ss, tmp, '"'); - f107.push_back(stof(tmp)); - - } - - } else { + std::vector itime(7, 0); + + while (getline(myFile, line)) { + std::stringstream ss(line); + // year + getline(ss, tmp, '-'); + itime[0] = stoi(tmp); + // month + getline(ss, tmp, '-'); + itime[1] = stoi(tmp); + // day + getline(ss, tmp, ' '); + itime[2] = stoi(tmp); + // hour + getline(ss, tmp, ':'); + itime[3] = stoi(tmp); + // minute + getline(ss, tmp, ' '); + itime[4] = stoi(tmp); + itime[5] = 0; + itime[6] = 0; + time.push_back(time_int_to_real(itime)); + + // f107 + getline(ss, tmp, '"'); + f107.push_back(stof(tmp)); + } // while + } else { iErr = 1; std::cout << "Couldn't file line #yyyy-MM etc in file " - << f107_file << "\n"; - + << f107_file << "\n"; } myFile.close(); - } return iErr; - } diff --git a/src/report.cpp b/src/report.cpp index 2331d006..29f35a2e 100644 --- a/src/report.cpp +++ b/src/report.cpp @@ -1,3 +1,5 @@ +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) +// Full license can be found in License.md #include #include @@ -7,22 +9,22 @@ // ----------------------------------------------------------------------- -// +// Initialize class Report // ----------------------------------------------------------------------- Report::Report() { - current_entry = ""; nEntries = 0; iVerbose = 0; divider = ">"; divider_length = divider.length(); - iLevel = 0; - + iLevel = 0; } // ----------------------------------------------------------------------- -// +// When the code enters a function, report this if the verbose level +// is high enough, and record the system start time of the entry, so that +// the run-time of the function can be recorded on exit. // ----------------------------------------------------------------------- void Report::enter(std::string input, int &iFunction) { @@ -34,10 +36,10 @@ void Report::enter(std::string input, int &iFunction) { int iEntry = -1; if (iFunction > -1) - if (current_entry == entries[iFunction].entry) + if (current_entry == entries[iFunction].entry) iEntry = iFunction; if (iEntry == -1) { - for (int i=0; i < nEntries; i++) + for (int i = 0; i < nEntries; i++) if (current_entry == entries[i].entry) iEntry = i; } if (iEntry == -1) { @@ -61,51 +63,48 @@ void Report::enter(std::string input, int &iFunction) { iLevel++; entries[iEntry].iLevel = iLevel; iCurrentFunction = iEntry; - - print(iLevel,"Entering function : "+current_entry); - + + print(iLevel, "Entering function : " + current_entry); } // ----------------------------------------------------------------------- -// +// This records the exit of the exit of the function, computing the +// runtime of the function. +// This assumes that the enter and exit functions are perfectly paired, +// so that the enter function sets the iCurrentFunction variable. // ----------------------------------------------------------------------- void Report::exit(std::string input) { int iEntry = -1; iEntry = iCurrentFunction; - //for (int i=0; i < nEntries; i++) - // if (current_entry == entries[i].entry) iEntry = i; - //std::cout << "iEntry : " << iEntry << " " << iCurrentFunction << " " << nEntries << "\n"; + if (iEntry == -1) { std::cout << "Report::exit Error!!! Could not find valid entry!\n"; std::cout << "current_entry : " << current_entry << "\n"; } else { + // Get current system time: unsigned long long now = std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(); + + // Calculate the difference in times, do get the total timing: entries[iEntry].timing_total = entries[iEntry].timing_total + float(now - entries[iEntry].timing_start)/1000.0; + + // Increment the total number of times that the function has been called: entries[iEntry].nTimes++; - // std::size_t pos = current_entry.find(input); - // if (pos == std::string::npos) { - // std::cout << "Report::exit Error!!! Could not find input : " << input << "!\n"; - // std::cout << " current_entry : " << current_entry << "\n"; - // } else { - // current_entry = current_entry.substr(0,pos-divider_length); - current_entry = current_entry.substr(0,entries[iEntry].iStringPosBefore); + // Pop the current function off the stack: + current_entry = current_entry.substr(0, entries[iEntry].iStringPosBefore); iCurrentFunction = entries[iEntry].iLastEntry; iLevel--; - // } } - } // ----------------------------------------------------------------------- -// +// Loop through all reported functions and report their run times // ----------------------------------------------------------------------- void Report::times() { - std::cout << "Timing Summary :\n"; for (int i=0; i < nEntries; i++) { std::cout << entries[i].entry << "\n"; @@ -114,43 +113,32 @@ void Report::times() { for (int j=0; j < entries[i].iLevel; j++) std::cout << " "; std::cout << "timing_total (s) : " << entries[i].timing_total << "\n"; } - } // ----------------------------------------------------------------------- -// +// Print string if verbose level is set high enough // ----------------------------------------------------------------------- void Report::print(int iLevel, std::string output_string) { - - if (iLevel <= iVerbose) { - - for (int iL=0;iL " << output_string << "\n"; - - } - + if (test_verbose(iLevel)) std::cout << output_string << "\n"; } // ----------------------------------------------------------------------- -// +// Test verbose level and print line starter if it is high enough // ----------------------------------------------------------------------- int Report::test_verbose(int iLevel) { - int iPass = 0; if (iLevel <= iVerbose) { iPass = 1; - for (int iL=0;iL "; } - return iPass; - } // ----------------------------------------------------------------------- -// +// // ----------------------------------------------------------------------- void Report::set_verbose(int input) { @@ -158,7 +146,7 @@ void Report::set_verbose(int input) { } // ----------------------------------------------------------------------- -// +// // ----------------------------------------------------------------------- int Report::get_verbose() { diff --git a/src/solver_chemistry.cpp b/src/solver_chemistry.cpp index 196b610b..2904f535 100644 --- a/src/solver_chemistry.cpp +++ b/src/solver_chemistry.cpp @@ -1,26 +1,16 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md +#include +using namespace arma; -float solver_chemistry(float old_density, float source, float loss, float dt) { - - float new_density; - - // This isn't a great way of doing things probably, but, what the heck: - - // if sources > losses, take an explicit time-step: - - if (source > loss) { - new_density = old_density + dt * (source - loss); - } else { - - // take implicit time-step: - - float normalized_loss = loss / (old_density + 1e-6); - new_density = (old_density + dt * source) / (1.0 + dt * normalized_loss); - - } +// ----------------------------------------------------------------------- +// This is a super-simple chemistry solver that takes an implicit time-step. +// We should create more sophisticated ones, but this is ok for now. +// ----------------------------------------------------------------------- +fcube solver_chemistry(fcube density, fcube source, fcube loss, float dt) { + fcube normalized_loss = loss / (density + 1e-6); + fcube new_density = (density + dt * source) / (1.0 + dt * normalized_loss); return new_density; - } diff --git a/src/solver_conduction.cpp b/src/solver_conduction.cpp index 74d39b0e..21ce0211 100644 --- a/src/solver_conduction.cpp +++ b/src/solver_conduction.cpp @@ -1,114 +1,90 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md -#include "../include/sizes.h" - -// Assume that lambda and front are scaled by radius squared! - -int solver_conduction(float value[nGeoAltsG], - float lambda[nGeoAltsG], - float front[nGeoAltsG], - float dt, - float dalt_lower[nGeoAltsG], - float *conduction) { - - int iErr = 0; - - int iAlt; - - float r[nGeoAltsG]; - float di[nGeoAltsG]; - float m[nGeoAltsG]; - float du[nGeoAltsG]; - float dl[nGeoAltsG]; - float du12[nGeoAltsG]; - float du22[nGeoAltsG]; - float lou[nGeoAltsG]; - float a[nGeoAltsG]; - float b[nGeoAltsG]; - float c[nGeoAltsG]; - float d[nGeoAltsG]; - float cp[nGeoAltsG]; - float dp[nGeoAltsG]; - float result[nGeoAltsG]; - - // Only solve for conduction in the middle section. - conduction[0] = 0.0; - conduction[nGeoAltsG-1] = 0.0; - - for (iAlt=1; iAlt < nGeoAltsG-1; iAlt++) { - - di[iAlt] = lambda[iAlt]; - m[iAlt] = dt / front[iAlt]; - - // This is the dalt upper (lower @ iAlt+1) - du[iAlt] = dalt_lower[iAlt+1]; - // This is the dalt lower - dl[iAlt] = dalt_lower[iAlt]; - - r[iAlt] = du[iAlt]/dl[iAlt]; - - du12[iAlt] = du[iAlt]*du[iAlt] * (1+r[iAlt])*(1+r[iAlt]); - du22[iAlt] = 0.5 * (dl[iAlt]*du[iAlt] + du[iAlt]*du[iAlt]); - lou[iAlt] = di[iAlt]/du22[iAlt]; - - } - - for (iAlt=2; iAlt < nGeoAltsG-2; iAlt++) { - - dl[iAlt] = - di[iAlt+1] - - di[iAlt-1] * r[iAlt] * r[iAlt] - - di[iAlt] * (1.0-r[iAlt]*r[iAlt]); - - a[iAlt] = di[iAlt]/du22[iAlt]*r[iAlt] - dl[iAlt]/du12[iAlt]*r[iAlt]*r[iAlt]; - c[iAlt] = di[iAlt]/du22[iAlt] + dl[iAlt]/du12[iAlt]; - b[iAlt] = - - 1.0/m[iAlt] - - di[iAlt]/du22[iAlt]*(1.0+r[iAlt]) - - dl[iAlt]/du12[iAlt]*(1.0-r[iAlt]*r[iAlt]); - d[iAlt] = -1.0*value[iAlt]/m[iAlt]; - - } - - // Lower BCs: - a[1] = 0.0; - b[1] = -1.0; - c[1] = 0.0; - d[1] = -1.0*value[1]; +#include +#include +using namespace arma; + +// This code solves the conduction equation in 1D. +// Some assumptions: +// - Assume that lambda and front are scaled by radius squared! +// - The spacing can be non-uniform. +// - At this point, the bottom BC is fixed, while the top BC is zero gradient +// - The dx variable is assumed to be distance between the CURRENT cell center +// (i) and the cell center of the cell BELOW the current one (i-1). + +fvec solver_conduction(fvec value, + fvec lambda, + fvec front, + float dt, + fvec dx) { + + int64_t nPts = value.n_elem; + + fvec di = lambda; + fvec m = dt / front; + + // These are to allow for a stretched grid: + // du is cell spacing in upper direction (so, lower, shifted by one): + fvec du(nPts); + du(span(0, nPts-2)) = dx(span(1, nPts-1)); + du(nPts-1) = du(nPts-2); + // dl is lower cell spacing: + fvec dl = dx; + fvec r = du/dl; + fvec du12 = du % du % (1 + r) % (1 + r); + fvec du22 = 0.5 * (dl % du + du % du); + fvec lou = di/du22; + + fvec conduction(nPts); + conduction.zeros(); + + int64_t i; + + for (i = 2; i < nPts-2; i++) + dl(i) = di(i+1) - di(i-1) * r(i) * r(i) - di(i) * (1.0-r(i)*r[i]); + + fvec a = di / du22 % r - dl / du12 % r % r; + fvec c = di / du22 + dl / du12; + fvec b = -1.0 / m - di / du22 % (1.0 + r) - dl / du12 % (1.0 - r % r); + fvec d = -1.0 * value / m; + + // Lower BCs (fixed value): + a(1) = 0.0; + b(1) = -1.0; + c(1) = 0.0; + d(1) = -1.0*value(1); // Upper BCs: - // This assumes a constant-gradient BC: - // (For neutral temperature, this isn't really needed, since it is iso-thermal - // in the upper thermosphere, but in the ionosphere, the electron and - // ion temperatures are typically sloped.) - iAlt = nGeoAltsG-2; - a[iAlt] = 1.0*( r[iAlt]*(1.0+r[iAlt])*di[iAlt]*m[iAlt]/du22[iAlt]); - b[iAlt] = -1.0*( 1.0 + r[iAlt]*(1+r[iAlt])*di[iAlt]*m[iAlt]/du22[iAlt]); - c[iAlt] = 0.0; - d[iAlt] = -1.0*value[iAlt]; - - cp[1] = c[1]/b[1]; - for (iAlt=2; iAlt <= nGeoAltsG-2; iAlt++) - cp[iAlt] = c[iAlt]/(b[iAlt]-cp[iAlt-1]*a[iAlt]); - - dp[1] = d[1]/b[1]; - for (iAlt=2; iAlt <= nGeoAltsG-2; iAlt++) - dp[iAlt] = (d[iAlt]-dp[iAlt-1]*a[iAlt])/(b[iAlt]-cp[iAlt-1]*a[iAlt]); - - result[nGeoAltsG-2] = dp[nGeoAltsG-2]; - for (iAlt=nGeoAltsG-3; iAlt>0; iAlt--) - result[iAlt] = dp[iAlt] - cp[iAlt]*result[iAlt+1]; - - conduction[0] = 0.0; - for (iAlt=1; iAlt 0; i--) + result(i) = dp(i) - cp(i)*result(i+1); + + conduction = result - value; + conduction(0) = 0.0; + conduction(1) = 0.0; + conduction(nPts-2) = 0.0; + conduction(nPts-1) = 0.0; + + return conduction; +} diff --git a/src/test.cpp b/src/test.cpp index d3174b48..a9a579a0 100644 --- a/src/test.cpp +++ b/src/test.cpp @@ -1,4 +1,4 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include @@ -12,11 +12,11 @@ int main() { // ------------------------------------------------------------ // Test time routines: // ------------------------------------------------------------ - + iErr = test_time_routines(); if (iErr == 0) std::cout << "Passed test_time_routines!\n"; - else std::cout << "Failed test_time_routines!\n"; - + else + std::cout << "Failed test_time_routines!\n"; + return iErr; - } diff --git a/src/time.cpp b/src/time.cpp index af549927..9ca60f33 100644 --- a/src/time.cpp +++ b/src/time.cpp @@ -1,4 +1,4 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include @@ -12,7 +12,7 @@ #include "../include/time_conversion.h" // ----------------------------------------------------------------------------- -// Initialize the time variables +// Instantiate the time variables // ----------------------------------------------------------------------------- Times::Times() { @@ -33,18 +33,16 @@ Times::Times() { } // ----------------------------------------------------------------------------- -// +// This is to initialize the time variables with real times // ----------------------------------------------------------------------------- void Times::set_times(std::vector itime) { - start = time_int_to_real(itime); current = start; iStep = -1; dt = 0; // This will initiate more variables: increment_time(); - } // ----------------------------------------------------------------------------- @@ -73,7 +71,7 @@ void Times::calc_dt() { } // ----------------------------------------------------------------------------- -// +// // ----------------------------------------------------------------------------- double Times::get_current() { @@ -81,7 +79,7 @@ double Times::get_current() { } // ----------------------------------------------------------------------------- -// +// // ----------------------------------------------------------------------------- double Times::get_end() { @@ -89,7 +87,7 @@ double Times::get_end() { } // ----------------------------------------------------------------------------- -// +// // ----------------------------------------------------------------------------- std::string Times::get_YMD_HMS() { @@ -97,7 +95,7 @@ std::string Times::get_YMD_HMS() { } // ----------------------------------------------------------------------------- -// +// // ----------------------------------------------------------------------------- double Times::get_intermediate() { @@ -105,7 +103,7 @@ double Times::get_intermediate() { } // ----------------------------------------------------------------------------- -// +// // ----------------------------------------------------------------------------- float Times::get_dt() { @@ -113,32 +111,29 @@ float Times::get_dt() { } // ----------------------------------------------------------------------------- -// +// // ----------------------------------------------------------------------------- float Times::get_orbittime() { return orbittime; } - // ----------------------------------------------------------------------------- -// +// // ----------------------------------------------------------------------------- double Times::get_julian_day() { return julian_day; } - // ----------------------------------------------------------------------------- -// +// // ----------------------------------------------------------------------------- void Times::set_end_time(std::vector itime) { end = time_int_to_real(itime); } - // ----------------------------------------------------------------------------- // Increment the time by dt and increment the iteration number (iStep) // ----------------------------------------------------------------------------- @@ -174,7 +169,7 @@ void Times::increment_time() { sYMD = std::string(tmp); sprintf(tmp, "%02d%02d%02d", hour, minute, second); sHMS = std::string(tmp); - + // Calculate Julian Day (day of year): jDay = day_of_year(year, month, day); @@ -182,20 +177,22 @@ void Times::increment_time() { julian_day = time_int_to_jday(iCurrent); // Calculate UT (in hours): - ut = float(iCurrent[3]) // hours - + float(iCurrent[4])/60.0 // minutes - + (float(iCurrent[5]) + float(iCurrent[6])/1000)/3600.0; + ut = static_cast(iCurrent[3]) // hours + + static_cast(iCurrent[4])/60.0 // minutes + + (static_cast(iCurrent[5]) + + static_cast(iCurrent[6])/1000)/3600.0; // Calculate orbital parameters based on E Standish, // Solar System Dynamics, JPL,. // No constant orbital speed assumption - float day_number = 367.0*float(iCurrent[0]) - - 7.0*(float(iCurrent[0])+(float(iCurrent[1])+9.0)/12.0)/4.0 - + 275.0 * float(iCurrent[1])/9.0 - + float(iCurrent[2]) + float day_number = 367.0 * static_cast(iCurrent[0]) + - 7.0*(static_cast(iCurrent[0]) + + (static_cast(iCurrent[1]) + 9.0) / 12.0) / 4.0 + + 275.0 * static_cast(iCurrent[1]) / 9.0 + + static_cast(iCurrent[2]) - 730531.5 - + ut/24.0; + + ut / 24.0; orbittime = day_number/36525.0; @@ -208,19 +205,24 @@ void Times::increment_time() { void Times::display() { + std::string units = " (s)"; time(&sys_time_current); - walltime = double(sys_time_current) - double(sys_time_start); - //cout << "Elapsed walltime : " << walltime << "\n"; + walltime = + static_cast(sys_time_current) - + static_cast(sys_time_start); + if (walltime > 120) { + if (walltime > 7200) { + walltime = walltime/3600.0; + units = " (h)"; + } else { + walltime = walltime/60.0; + units = " (m)"; + } + } - //cout << "current time (double) : " << current << "\n"; std::cout << "Current Time : "; - for (int i=0;i<7;i++) std::cout << iCurrent[i] << " "; - std::cout << "\n"; - // cout << " (as julian day): " << julian_day << "\n"; - // cout << " (as julian day since j2000): " << julian_day-j2000 << "\n"; - - return; - + display_itime(iCurrent); + std::cout << "Wall Time : " << walltime << units; } // ----------------------------------------------------------------------------- diff --git a/src/time_conversion.cpp b/src/time_conversion.cpp index 4a5116ec..0c9c0edf 100644 --- a/src/time_conversion.cpp +++ b/src/time_conversion.cpp @@ -1,3 +1,5 @@ +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) +// Full license can be found in License.md #include #include @@ -9,10 +11,9 @@ // ----------------------------------------------------------------------------- void display_itime(std::vector itime) { - - for (int i=0; i itime) { // ----------------------------------------------------------------------------- int day_of_year(int year, int month, int day) { - - int doy; - - doy = 0; - - for (int i=1; i itime) { int nDays = day_of_year(itime[0], itime[1], itime[2])-1; double timereal = - double(itime[6])/1000.0d + // milliseconds - double(itime[5]) + // seconds - double(itime[4]) * seconds_per_minute + - double(itime[3]) * seconds_per_hour + - double(nDays+nLeaps) * seconds_per_day + - double(nYears) * seconds_per_year; + static_cast(itime[6])/1000.0d + // milliseconds + static_cast(itime[5]) + // seconds + static_cast(itime[4]) * seconds_per_minute + + static_cast(itime[3]) * seconds_per_hour + + static_cast(nDays+nLeaps) * seconds_per_day + + static_cast(nYears) * seconds_per_year; return timereal; - } // ----------------------------------------------------------------------------- @@ -64,20 +58,27 @@ double time_int_to_real(std::vector itime) { void time_real_to_int(double timereal, std::vector &itime) { - int nYears = int(timereal/seconds_per_year); + int nYears = static_cast(timereal/seconds_per_year); int nLeaps = nYears/4; - int nDays = int((timereal - (double(nYears)*seconds_per_year))/seconds_per_day); + int nDays = static_cast ((timereal - + (static_cast(nYears) * + seconds_per_year)) / + seconds_per_day); // This can happen in some circumstances, like the first few days of the year: if (nDays < nLeaps) { - nYears = (timereal - (double(nLeaps)*seconds_per_day))/seconds_per_year; + nYears = (timereal - (static_cast(nLeaps) * + seconds_per_day))/seconds_per_year; nLeaps = nYears/4; - nDays = (timereal - (double(nYears)*seconds_per_year))/seconds_per_day; + nDays = (timereal - (static_cast(nYears) * + seconds_per_year))/seconds_per_day; // This should be very rare: if (nDays < nLeaps) { - nYears = (timereal - (double(nLeaps)*seconds_per_day))/seconds_per_year; + nYears = (timereal - (static_cast(nLeaps) * + seconds_per_day))/seconds_per_year; nLeaps = nYears/4; - nDays = (timereal - (double(nYears)*seconds_per_year))/seconds_per_day; + nDays = (timereal - (static_cast(nYears) * + seconds_per_year))/seconds_per_day; } } @@ -86,20 +87,21 @@ void time_real_to_int(double timereal, std::vector &itime) { // Calculate how much time is left, after subtracting off years and days: double timeleft = timereal - - double(nYears) * seconds_per_year - - double(nDays + nLeaps) * seconds_per_day; + - static_cast(nYears) * seconds_per_year + - static_cast(nDays + nLeaps) * seconds_per_day; // Calculate hours and subtract them: int nHours = timeleft/seconds_per_hour; - timeleft -= double(nHours) * seconds_per_hour; + timeleft -= static_cast(nHours) * seconds_per_hour; // Calculate minutes and subtract them: int nMinutes = timeleft/seconds_per_minute; - timeleft -= double(nMinutes) * seconds_per_minute; + timeleft -= static_cast(nMinutes) * seconds_per_minute; // Calculate milliseconds: int nSeconds = timeleft; - int nMillis = int((timeleft-double(nSeconds))*1000.0); + int nMillis = static_cast((timeleft - + static_cast(nSeconds)) * 1000.0); itime[0] = nYears + reference_year; @@ -109,9 +111,9 @@ void time_real_to_int(double timereal, std::vector &itime) { int nMonths = 1; int iLeap = 0; if (itime[0] % 4 == 0) iLeap=1; - std::vector add_on_days {0,iLeap,0,0,0,0,0,0,0,0,0,0}; - while (nDays > (days_of_month[nMonths-1]+add_on_days[nMonths-1])) { - nDays = nDays - (days_of_month[nMonths-1]+add_on_days[nMonths-1]); + std::vector add_on_days {0, iLeap, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + while (nDays > (days_of_month[nMonths-1] + add_on_days[nMonths-1])) { + nDays = nDays - (days_of_month[nMonths-1] + add_on_days[nMonths-1]); nMonths++; } itime[1] = nMonths; @@ -120,7 +122,6 @@ void time_real_to_int(double timereal, std::vector &itime) { itime[4] = nMinutes; itime[5] = nSeconds; itime[6] = nMillis; - } // ----------------------------------------------------------------------------- @@ -128,16 +129,12 @@ void time_real_to_int(double timereal, std::vector &itime) { // ----------------------------------------------------------------------------- double time_int_to_jday(std::vector itime) { - // For this, we are going to use a relationship between our // time reference system and the Julian Day reference system. - double our_time = time_int_to_real(itime); double our_time_in_days = our_time/seconds_per_day; double jd = julian_day_of_reference + our_time_in_days; - return jd; - } // ----------------------------------------------------------------------------- @@ -148,41 +145,33 @@ int test_time_routines() { std::vector itime {1970, 1, 1, 0, 0, 0, 0}; double timeout, timecheck; - int i; int iErr; iErr = 0; - + timeout = time_int_to_real(itime); // This assumes a reference time of Jan. 1, 1965 0000 UT: timecheck = 1.5776640e+08; - - for (i=0; i<6; i++) std::cout << itime[i] << " "; + + display_itime(itime); std::cout << " --> " << timeout << " compares to : " << timecheck << "\n"; if (abs(timecheck-timeout) > 1.0) { iErr = 1; std::cout << "Fails!!!\n"; - } else std::cout << "Passes!!!\n"; + } else { + std::cout << "Passes!!!\n"; + } time_real_to_int(timecheck, itime); - for (i=0; i<6; i++) std::cout << itime[i] << " "; - std::cout << "\n"; + display_itime(itime); - itime[0]=2000; - itime[1]=1; - itime[2]=1; + itime[0] = 2000; + itime[1] = 1; + itime[2] = 1; double jd_test = time_int_to_jday(itime); std::cout << "Test Julian Day = " << jd_test << "\n"; std::cout << "Julian Day 2000 = " << j2000 << "\n"; - + return iErr; - } - - - -// ----------------------------------------------------------------------------- -// -// ----------------------------------------------------------------------------- - diff --git a/src/transform.cpp b/src/transform.cpp index a80c568e..9962cee1 100644 --- a/src/transform.cpp +++ b/src/transform.cpp @@ -1,26 +1,68 @@ -// (c) 2020, the Aether Development Team (see doc/dev_team.md for members) +// Copyright 2020, the Aether Development Team (see doc/dev_team.md for members) // Full license can be found in License.md #include #include +#include #include "../include/sizes.h" #include "../include/grid.h" +using namespace arma; + + +// ----------------------------------------------------------------------- +// copy from armidillo cube to 3d c-native array +// ----------------------------------------------------------------------- + +void copy_cube_to_array(fcube cube_in, + float *array_out) { + + int64_t nX = cube_in.n_rows; + int64_t nY = cube_in.n_cols; + int64_t nZ = cube_in.n_slices; + int64_t iX, iY, iZ, index; + + for (iX = 0; iX < nX; iX++) { + for (iY = 0; iY < nY; iY++) { + for (iZ = 0; iZ < nZ; iZ++) { + index = iX*nY*nZ + iY*nZ + iZ; + array_out[index] = cube_in(iX, iY, iZ); + } + } + } +} + + + // ----------------------------------------------------------------------- // Transform Longitude, Latitude, Radius to X, Y, Z +// Use armidillo cubes // ----------------------------------------------------------------------- -void transform_llr_to_xyz(float llr_in[3], float xyz_out[3]) { +void transform_llr_to_xyz_3d(fcube lat3d, + fcube lon3d, + fcube r3d, + fcube &x3d, + fcube &y3d, + fcube &z3d) { + + x3d = r3d % cos(lat3d) % cos(lon3d); + y3d = r3d % cos(lat3d) % sin(lon3d); + z3d = r3d % sin(lat3d); +} + +// ----------------------------------------------------------------------- +// Transform Longitude, Latitude, Radius to X, Y, Z +// ----------------------------------------------------------------------- +void transform_llr_to_xyz(float llr_in[3], float xyz_out[3]) { // llr_in[0] = longitude (in radians) // llr_in[1] = latitude (in radians) // llr_in[2] = radius - xyz_out[0] = llr_in[2] * cos(llr_in[1]) * cos(llr_in[0]); xyz_out[1] = llr_in[2] * cos(llr_in[1]) * sin(llr_in[0]); xyz_out[2] = llr_in[2] * sin(llr_in[1]); - } // ----------------------------------------------------------------------- @@ -29,16 +71,11 @@ void transform_llr_to_xyz(float llr_in[3], float xyz_out[3]) { // ----------------------------------------------------------------------- void transform_rot_z(float xyz_in[3], float angle_in, float xyz_out[3]) { - float ca = cos(angle_in); float sa = sin(angle_in); - xyz_out[0] = xyz_in[0] * ca + xyz_in[1] * sa; xyz_out[1] = -xyz_in[0] * sa + xyz_in[1] * ca; xyz_out[2] = xyz_in[2]; - - return; - } // ----------------------------------------------------------------------- @@ -47,16 +84,11 @@ void transform_rot_z(float xyz_in[3], float angle_in, float xyz_out[3]) { // ----------------------------------------------------------------------- void transform_rot_y(float xyz_in[3], float angle_in, float xyz_out[3]) { - float ca = cos(angle_in); float sa = sin(angle_in); - xyz_out[0] = xyz_in[0] * ca - xyz_in[2] * sa; xyz_out[1] = xyz_in[1]; xyz_out[2] = xyz_in[0] * sa + xyz_in[2] * ca; - - return; - } // ----------------------------------------------------------------------- @@ -64,11 +96,8 @@ void transform_rot_y(float xyz_in[3], float angle_in, float xyz_out[3]) { // ----------------------------------------------------------------------- void transform_float_vector_to_array(std::vector input, - float output[3]) { - - for (int i=0; i<3; i++) output[i] = input[i]; - - return; + float output[3]) { + for (int i = 0; i < 3; i++) output[i] = input[i]; } // ----------------------------------------------------------------------- @@ -76,16 +105,17 @@ void transform_float_vector_to_array(std::vector input, // ----------------------------------------------------------------------- void transform_vector_xyz_to_env(float xyz_in[3], - float lon, - float lat, - float env_out[3]) { - - env_out[2] = xyz_in[0] * cos(lat)*cos(lon) + xyz_in[1] * cos(lat)*sin(lon) + xyz_in[2]*sin(lat); - env_out[1] = -(xyz_in[0] * sin(lat)*cos(lon) + xyz_in[1] * sin(lat)*sin(lon) - xyz_in[2]*cos(lat)); - env_out[0] = - xyz_in[0] * sin(lon) + xyz_in[1] * cos(lon); - - return; - + float lon, + float lat, + float env_out[3]) { + + env_out[2] = xyz_in[0] * cos(lat)*cos(lon) + + xyz_in[1] * cos(lat)*sin(lon) + xyz_in[2]*sin(lat); + env_out[1] = -(xyz_in[0] * sin(lat)*cos(lon) + + xyz_in[1] * sin(lat)*sin(lon) - + xyz_in[2]*cos(lat)); + env_out[0] = - xyz_in[0] * sin(lon) + + xyz_in[1] * cos(lon); } // ----------------------------------------------------------------------- @@ -102,13 +132,9 @@ void transform_vector_xyz_to_env(float xyz_in[3], // ----------------------------------------------------------------------- void vector_diff(float vect_in_1[3], - float vect_in_2[3], - float vect_out[3]) { - - for (int i=0; i<3; i++) vect_out[i] = vect_in_1[i] - vect_in_2[i]; - - return; - + float vect_in_2[3], + float vect_out[3]) { + for (int i = 0; i < 3; i++) vect_out[i] = vect_in_1[i] - vect_in_2[i]; } // ----------------------------------------------------------------------- @@ -116,11 +142,11 @@ void vector_diff(float vect_in_1[3], // ----------------------------------------------------------------------- void get_vector_component(float *vector_in_v3gc, - int iComponent, - int IsGeoGrid, - float *component_out_s3gc) { + int iComponent, + int IsGeoGrid, + float *component_out_s3gc) { - long nLons, nLats, nAlts, iLon, iLat, iAlt, index, indexv; + int64_t nLons, nLats, nAlts, iLon, iLat, iAlt, index, indexv; if (IsGeoGrid) { nLons = nGeoLonsG; @@ -135,21 +161,15 @@ void get_vector_component(float *vector_in_v3gc, for (iLon = 0; iLon < nLons; iLon++) { for (iLat = 0; iLat < nLats; iLat++) { for (iAlt = 0; iAlt < nAlts; iAlt++) { - - if (IsGeoGrid) { - index = ijk_geo_s3gc(iLon,iLat,iAlt); - indexv = ijkl_geo_v3gc(iLon,iLat,iAlt,iComponent); - } else { - index = ijk_mag_s3gc(iLon,iLat,iAlt); - indexv = ijkl_mag_v3gc(iLon,iLat,iAlt,iComponent); - } - - component_out_s3gc[index] = vector_in_v3gc[indexv]; - + if (IsGeoGrid) { + index = ijk_geo_s3gc(iLon, iLat, iAlt); + indexv = ijkl_geo_v3gc(iLon, iLat, iAlt, iComponent); + } else { + index = ijk_mag_s3gc(iLon, iLat, iAlt); + indexv = ijkl_mag_v3gc(iLon, iLat, iAlt, iComponent); + } + component_out_s3gc[index] = vector_in_v3gc[indexv]; } } } - - return; - }