From 7645eb596871257a191f4e6c48d94af1f1757e39 Mon Sep 17 00:00:00 2001 From: Matt Williamson <87771120+mattw-nws@users.noreply.github.com> Date: Wed, 2 Aug 2023 17:11:08 +0000 Subject: [PATCH] Remove tshirt, tshirt_c and related tests and data --- CMakeLists.txt | 2 - .../catchment/Formulation_Constructors.hpp | 4 - .../catchment/Tshirt_C_Realization.hpp | 345 ------ .../catchment/Tshirt_Realization.hpp | 159 --- models/tshirt/include/Tshirt.h | 223 ---- models/tshirt/include/TshirtErrorCodes.h | 17 - models/tshirt/include/tshirt_c.h | 151 --- models/tshirt/include/tshirt_fluxes.h | 28 - models/tshirt/include/tshirt_params.h | 89 -- models/tshirt/include/tshirt_state.h | 26 - src/NGen.cpp | 1 - src/core/CMakeLists.txt | 3 - src/models/tshirt/CMakeLists.txt | 15 - src/models/tshirt/Tshirt.cpp | 369 ------ src/models/tshirt/tshirt_c.cpp | 616 ---------- src/realizations/catchment/CMakeLists.txt | 2 - .../catchment/Tshirt_C_Realization.cpp | 684 ----------- .../catchment/Tshirt_Realization.cpp | 337 ----- src/routing/CMakeLists.txt | 1 - test/CMakeLists.txt | 20 +- test/data/model/tshirt/expected_results.csv | 999 --------------- test/models/tshirt/include/TshirtTest.cpp | 63 - .../realizations/Formulation_Manager_Test.cpp | 1 - .../catchments/Tshirt_C_Realization_Test.cpp | 1079 ----------------- 24 files changed, 2 insertions(+), 5232 deletions(-) delete mode 100644 include/realizations/catchment/Tshirt_C_Realization.hpp delete mode 100644 include/realizations/catchment/Tshirt_Realization.hpp delete mode 100644 models/tshirt/include/Tshirt.h delete mode 100644 models/tshirt/include/TshirtErrorCodes.h delete mode 100644 models/tshirt/include/tshirt_c.h delete mode 100644 models/tshirt/include/tshirt_fluxes.h delete mode 100644 models/tshirt/include/tshirt_params.h delete mode 100644 models/tshirt/include/tshirt_state.h delete mode 100644 src/models/tshirt/CMakeLists.txt delete mode 100644 src/models/tshirt/Tshirt.cpp delete mode 100644 src/models/tshirt/tshirt_c.cpp delete mode 100644 src/realizations/catchment/Tshirt_C_Realization.cpp delete mode 100644 src/realizations/catchment/Tshirt_Realization.cpp delete mode 100644 test/data/model/tshirt/expected_results.csv delete mode 100644 test/models/tshirt/include/TshirtTest.cpp delete mode 100644 test/realizations/catchments/Tshirt_C_Realization_Test.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 7c0733e85c..7e196491bd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -225,7 +225,6 @@ if(NGEN_WITH_SQLITE3) endif() add_subdirectory("src/realizations/catchment") -add_subdirectory("src/models/tshirt") add_subdirectory("src/models/kernels/reservoir") add_subdirectory("src/models/kernels/evapotranspiration") add_subdirectory("src/forcing") @@ -236,7 +235,6 @@ target_link_libraries(ngen PUBLIC #NGen::core_catchment_giuh NGen::core_nexus NGen::geojson - NGen::models_tshirt NGen::realizations_catchment NGen::kernels_reservoir NGen::kernels_evapotranspiration diff --git a/include/realizations/catchment/Formulation_Constructors.hpp b/include/realizations/catchment/Formulation_Constructors.hpp index dfe5ada987..fe106de361 100644 --- a/include/realizations/catchment/Formulation_Constructors.hpp +++ b/include/realizations/catchment/Formulation_Constructors.hpp @@ -9,8 +9,6 @@ #include // Formulations -#include "Tshirt_Realization.hpp" -#include "Tshirt_C_Realization.hpp" #include "Bmi_Cpp_Formulation.hpp" #include "Bmi_C_Formulation.hpp" #include "Bmi_Fortran_Formulation.hpp" @@ -48,8 +46,6 @@ namespace realization { #ifdef ACTIVATE_PYTHON {"bmi_python", create_formulation_constructor()}, #endif // ACTIVATE_PYTHON - {"tshirt", create_formulation_constructor()}, - {"tshirt_c", create_formulation_constructor()}, #ifdef NGEN_LSTM_TORCH_LIB_ACTIVE , {"lstm", create_formulation_constructor()} diff --git a/include/realizations/catchment/Tshirt_C_Realization.hpp b/include/realizations/catchment/Tshirt_C_Realization.hpp deleted file mode 100644 index d613a615bc..0000000000 --- a/include/realizations/catchment/Tshirt_C_Realization.hpp +++ /dev/null @@ -1,345 +0,0 @@ -#ifndef NGEN_TSHIRT_C_REALIZATION_HPP -#define NGEN_TSHIRT_C_REALIZATION_HPP - -#include "Catchment_Formulation.hpp" -#include "CsvPerFeatureForcingProvider.hpp" -#include "core/catchment/HY_CatchmentArea.hpp" -#include "tshirt_params.h" -#include "GiuhJsonReader.h" -#include "tshirt_c.h" -#include -#include -#include -#include - -#define OUT_VAR_BASE_FLOW "base_flow" -#define OUT_VAR_GIUH_RUNOFF "giuh_runoff" -#define OUT_VAR_LATERAL_FLOW "lateral_flow" -#define OUT_VAR_RAINFALL "rainfall" -#define OUT_VAR_SURFACE_RUNOFF "surface_runoff" -#define OUT_VAR_TOTAL_DISCHARGE "total_discharge" - -using namespace tshirt; - -namespace realization { - - class Tshirt_C_Realization : public Catchment_Formulation { - - public: - - typedef long time_step_t; - - /** - * Constructor for initializing when a #tshirt::tshirt_params struct is passed for parameters and a - * #giuh::GiuhJsonReader is passed for obtaining GIUH ordinates. - * - * @param forcing_config - * @param output_stream - * @param soil_storage - * @param groundwater_storage - * @param storage_values_are_ratios Whether the storage values are given as proportional amounts of the max (or, - * when `false`, express amounts with units of meters). - * @param catchment_id - * @param giuh_json_reader - * @param params - * @param nash_storage - */ - Tshirt_C_Realization(forcing_params forcing_config, - utils::StreamHandler output_stream, - double soil_storage, - double groundwater_storage, - bool storage_values_are_ratios, - std::string catchment_id, - giuh::GiuhJsonReader &giuh_json_reader, - tshirt::tshirt_params params, - const std::vector &nash_storage); - - /** - * Constructor for initializing when a #tshirt::tshirt_params struct is passed for parameters and GIUH ordinates - * are available directly and passed within a vector. - * - * @param forcing_config - * @param output_stream - * @param soil_storage - * @param groundwater_storage - * @param storage_values_are_ratios Whether the storage values are given as proportional amounts of the max (or, - * when `false`, express amounts with units of meters). - * @param catchment_id - * @param guih_ordinates - * @param params - * @param nash_storage - */ - Tshirt_C_Realization(forcing_params forcing_config, - utils::StreamHandler output_stream, - double soil_storage, - double groundwater_storage, - bool storage_values_are_ratios, - std::string catchment_id, - std::vector giuh_ordinates, - tshirt::tshirt_params params, - const std::vector &nash_storage); - - /** - * Constructor for when model parameters are provided individually instead of within encapsulating struct and a - * #giuh::GiuhJsonReader is passed for obtaining GIUH ordinates. - * - * @param forcing_config - * @param output_stream - * @param soil_storage - * @param groundwater_storage - * @param storage_values_are_ratios Whether the storage values are given as proportional amounts of the max (or, - * when `false`, express amounts with units of meters). - * @param catchment_id - * @param giuh_json_reader - * @param maxsmc - * @param wltsmc - * @param satdk - * @param satpsi - * @param slope - * @param b - * @param multiplier - * @param alpha_fc - * @param Klf - * @param Kn - * @param nash_n - * @param Cgw - * @param expon - * @param max_gw_storage - * @param nash_storage - */ - Tshirt_C_Realization( - forcing_params forcing_config, - utils::StreamHandler output_stream, - double soil_storage, - double groundwater_storage, - bool storage_values_are_ratios, - std::string catchment_id, - giuh::GiuhJsonReader &giuh_json_reader, - double maxsmc, - double wltsmc, - double satdk, - double satpsi, - double slope, - double b, - double multiplier, - double alpha_fc, - double Klf, - double Kn, - int nash_n, - double Cgw, - double expon, - double max_gw_storage, - const std::vector &nash_storage); - - Tshirt_C_Realization(std::string id, shared_ptr forcing_provider, utils::StreamHandler output_stream); - - //[[deprecated]] - Tshirt_C_Realization(std::string id, forcing_params forcing_config, utils::StreamHandler output_stream); - - virtual ~Tshirt_C_Realization(); - - /** - * Return ``0``, as (for now) this type does not otherwise include ET within its calculations. - * - * @return ``0`` - */ - double calc_et() override; - - void create_formulation(boost::property_tree::ptree &config, geojson::PropertyMap *global = nullptr) override; - void create_formulation(geojson::PropertyMap properties) override; - - /** - * Run model formulation calculations for the next time step using the given input flux value in meters per - * second. - * - * The backing model works on a collection of time steps by receiving an associated collection of input fluxes. - * This is implemented by putting this input flux in a single-value vector and using it as the - * arg to a nested call to ``run_formulation_for_timesteps``, returning that result code. - * - * @param input_flux Input flux (typically expected to be just precipitation) in meters per second. - * @param t_delta_s The size of the time step in seconds - * @return The result code from the execution of the model time step calculations. - */ - int run_formulation_for_timestep(double input_flux, time_step_t t_delta_s); - - /** - * Run model formulation calculations for a series of time steps using the given collection of input flux values - * in meters per second. - * - * @param input_fluxes Ordered, per-time-step input flux (typically expected to be just precipitation) in meters - * per second. - * @param t_delta_s The sizes of each of the time steps in seconds - * @return The result code from the execution of the model time step calculations. - */ - int run_formulation_for_timesteps(std::vector input_fluxes, std::vector t_deltas_s); - - std::string get_formulation_type() override; - - /** - * Execute the backing model formulation for the given time step, where it is of the specified size, and - * return the total discharge. - * - * Any inputs and additional parameters must be made available as instance members. - * - * Types should clearly document the details of their particular response output. - * - * @param t_index The index of the time step for which to run model calculations. - * @param d_delta_s The duration, in seconds, of the time step for which to run model calculations. - * @return The total discharge of the model for this time step. - */ - double get_response(time_step_t t_index, time_step_t t_delta) override; - - // TODO: probably need to do better than this for granting and protecting access - double get_latest_flux_base_flow(); - - double get_latest_flux_giuh_runoff(); - - double get_latest_flux_lateral_flow(); - - double get_latest_flux_surface_runoff(); - - double get_latest_flux_total_discharge(); - - // TODO: look at moving these to the Formulation superclass - /** - * Get the number of output data variables made available from the calculations for enumerated time steps. - * - * @return The number of output data variables made available from the calculations for enumerated time steps. - */ - int get_output_item_count(); - - /** - * Get the names of the output data variables that are available from calculations for enumerated time steps. - * - * The method must deterministically return output variable names in the same order each time. - * - * @return The names of the output data variables that are available from calculations for enumerated time steps. - */ - const std::vector& get_output_var_names(); - - /** - * Get a copy of the values for the given output variable at all available time steps. - * - * @param name - * @return A vector containing copies of the output value of the variable, indexed by time step. - */ - // TODO: note that for this type, these are all doubles, but may need template function once moving to superclass - std::vector get_value(const std::string& name); - - /** - * Get a header line appropriate for a file made up of entries from this type's implementation of - * ``get_output_line_for_timestep``. - * - * Note that like the output generating function, this line does not include anything for time step. - * - * @return An appropriate header line for this type. - */ - std::string get_output_header_line(std::string delimiter=",") override; - - /** - * Get the values making up the header line from get_output_header_line(), but organized as a vector of strings. - * - * @return The values making up the header line from get_output_header_line() organized as a vector. - */ - const std::vector& get_output_header_fields(); - - /** - * Get a delimited string with all the output variable values for the given time step. - * - * This method is useful for preparing calculated data in a representation useful for output files, such as - * CSV files. - * - * The resulting string contains only the calculated output values for the time step, and not the time step - * index itself. - * - * An empty string is returned if the time step value is not in the range of valid time steps for which there - * are calculated values for all variables. - * - * The default delimiter is a comma. - * - * @param timestep The time step for which data is desired. - * @return A delimited string with all the output variable values for the given time step. - */ - std::string get_output_line_for_timestep(int timestep, std::string delimiter=",") override; - - private: - /** Id of associated catchment. */ - std::string catchment_id; - /** - * Struct for containing model parameters for Tshirt model, using "internal" framework implementation for such a - * struct. - */ - std::shared_ptr params; - /** Struct from C-style Tshirt formulation for holding soil parameter values. */ - NWM_soil_parameters c_soil_params; - - /** Vector of GIUH CDF ordinates. */ - std::vector giuh_cdf_ordinates; - /** Vector to serve as runoff queue for GIUH convolution calculations. */ - std::vector giuh_runoff_queue_per_timestep; - std::vector nash_storage; - - aorc_forcing_data c_aorc_params; - - // TODO: might want to consider having an initial time step value for reference (implied size is 1 hour) - // TODO: this probably need to be converted to use a different fluxes type that can be dealt with externally - std::vector> fluxes; - - conceptual_reservoir groundwater_conceptual_reservoir; - - std::vector OUTPUT_VARIABLE_NAMES = { - OUT_VAR_RAINFALL, - OUT_VAR_SURFACE_RUNOFF, - OUT_VAR_GIUH_RUNOFF, - OUT_VAR_LATERAL_FLOW, - OUT_VAR_BASE_FLOW, - OUT_VAR_TOTAL_DISCHARGE - }; - - std::vector OUTPUT_HEADER_FIELDS = { - "Rainfall", - "Direct Runoff", - "GIUH Runoff", - "Lateral Flow", - "Base Flow", - "Total Discharge" - }; - - conceptual_reservoir soil_conceptual_reservoir; - std::vector REQUIRED_PARAMETERS = { - "maxsmc", - "wltsmc", - "satdk", - "satpsi", - "slope", - "scaled_distribution_fn_shape_parameter", - "multiplier", - "alpha_fc", - "Klf", - "Kn", - "nash_n", - "Cgw", - "expon", - "max_groundwater_storage_meters", - "nash_storage", - "soil_storage_percentage", - "groundwater_storage_percentage", - "giuh" - }; - - std::function get_output_var_flux_extraction_func(const std::string& var_name); - - const std::vector& get_required_parameters() override; - - void init_ground_water_reservoir(double storage, bool storage_values_are_ratios); - - void init_soil_reservoir(double storage, bool storage_values_are_ratios); - - static double init_reservoir_storage(bool is_ratio, double amount, double max_amount); - - void sync_c_storage_params(); - - }; -} - -#endif diff --git a/include/realizations/catchment/Tshirt_Realization.hpp b/include/realizations/catchment/Tshirt_Realization.hpp deleted file mode 100644 index 9f8bfa4be2..0000000000 --- a/include/realizations/catchment/Tshirt_Realization.hpp +++ /dev/null @@ -1,159 +0,0 @@ -#ifndef NGEN_TSHIRT_REALIZATION_HPP -#define NGEN_TSHIRT_REALIZATION_HPP - -#include "Catchment_Formulation.hpp" -#include "CsvPerFeatureForcingProvider.hpp" -#include -#include "GIUH.hpp" -#include "GiuhJsonReader.h" -#include "tshirt/include/Tshirt.h" -#include "tshirt/include/tshirt_params.h" -#include -#include - -namespace realization { - - class Tshirt_Realization : public Catchment_Formulation { - - public: - - typedef long time_step_t; - - /* - Tshirt_Realization(forcing_params forcing_config, - utils::StreamHandler output_stream, - double soil_storage_meters, - double groundwater_storage_meters, - std::string catchment_id, - giuh::GiuhJsonReader &giuh_json_reader, - tshirt::tshirt_params params, - const std::vector& nash_storage, - time_step_t t); - - Tshirt_Realization( - forcing_params forcing_config, - utils::StreamHandler output_stream, - double soil_storage_meters, - double groundwater_storage_meters, - std::string catchment_id, - giuh::GiuhJsonReader &giuh_json_reader, - double maxsmc, - double wltsmc, - double satdk, - double satpsi, - double slope, - double b, - double multiplier, - double alpha_fc, - double Klf, - double Kn, - int nash_n, - double Cgw, - double expon, - double max_gw_storage, - const std::vector& nash_storage, - time_step_t t - ); - */ - - Tshirt_Realization( - std::string id, - shared_ptr forcing_provider, - utils::StreamHandler output_stream - ) : Catchment_Formulation(id, forcing_provider, output_stream) { - } - - void set_giuh_kernel(std::shared_ptr reader); - - virtual ~Tshirt_Realization(){}; - - double calc_et() override; - - /** - * Execute the backing model formulation for the given time step, where it is of the specified size, and - * return the total discharge. - * - * Function reads input precipitation from ``forcing`` member variable. It also makes use of the params struct - * for ET params accessible via ``get_et_params``. - * - * @param t_index The index of the time step for which to run model calculations. - * @param d_delta_s The duration, in seconds, of the time step for which to run model calculations. - * @return The total discharge for this time step. - */ - double get_response(time_step_t t_index, time_step_t t_delta_s) override; - - std::string get_formulation_type() override { - return "tshirt"; - } - - /** - * Get a formatted line of output values for the given time step as a delimited string. - * - * For this type, the output consists of only the total discharge amount per time step; i.e., the same value - * that was returned by ``get_response``. - * - * This method is useful for preparing calculated data in a representation useful for output files, such as - * CSV files. - * - * The resulting string will contain calculated values for applicable output variables for the particular - * formulation, as determined for the given time step. However, the string will not contain any - * representation of the time step itself. - * - * An empty string is returned if the time step value is not in the range of valid time steps for which there - * are calculated values for all variables. - * - * The default delimiter is a comma. - * - * @param timestep The time step for which data is desired. - * @return A delimited string with all the output variable values for the given time step. - */ - std::string get_output_line_for_timestep(int timestep, std::string delimiter=",") override; - - void create_formulation(boost::property_tree::ptree &config, geojson::PropertyMap *global = nullptr) override; - - void create_formulation(geojson::PropertyMap properties) override; - - const std::vector& get_required_parameters() override { - return REQUIRED_PARAMETERS; - } - - private: - std::string catchment_id; - std::unordered_map> state; - std::unordered_map> fluxes; - std::unordered_map > cascade_backing_storage; - tshirt::tshirt_params *params; - std::unique_ptr model; - std::shared_ptr giuh_kernel; - - //The delta time (dt) this instance is configured to use - time_step_t dt; - - std::vector REQUIRED_PARAMETERS = { - "maxsmc", - "wltsmc", - "satdk", - "satpsi", - "slope", - "scaled_distribution_fn_shape_parameter", - "multiplier", - "alpha_fc", - "Klf", - "Kn", - "nash_n", - "Cgw", - "expon", - "max_groundwater_storage_meters", - "nash_storage", - "soil_storage_percentage", - "groundwater_storage_percentage", - "timestep", - "giuh" - }; - - }; - -} - - -#endif //NGEN_TSHIRT_REALIZATION_HPP diff --git a/models/tshirt/include/Tshirt.h b/models/tshirt/include/Tshirt.h deleted file mode 100644 index d9bb9fabf1..0000000000 --- a/models/tshirt/include/Tshirt.h +++ /dev/null @@ -1,223 +0,0 @@ -#ifndef TSHIRT_H -#define TSHIRT_H - -#include "all.h" -#include "schaake_partitioning.hpp" -#include "Constants.h" -#include "reservoir/Reservoir.hpp" -#include "Pdm03.h" -#include "GIUH.hpp" -#include "reservoir/Reservoir_Exponential_Outlet.hpp" -#include "reservoir/Reservoir_Linear_Outlet.hpp" -#include "reservoir/Reservoir_Outlet.hpp" -#include "tshirt_fluxes.h" -#include "tshirt_params.h" -#include "tshirt_state.h" -#include -#include -#include -#include -#include - -namespace tshirt { - - /** - * Tshirt model class. - * - * Object oriented implementation of the Tshirt hydrological model, as documented at: - * https://github.com/NOAA-OWP/ngen/blob/master/doc/T-shirt_model_description.pdf - * - * The core logic for the model is implemented within the `run` function. The resulting state calculated by a call - * to `run` is stored by a struct referenced via the `current_state` member. A `previous_state` member also - * exists for holding the state prior to the most recent call to `run`. Calculated fluxes are held by the `fluxes` - * member, for which there is not 'previous' analog. - * - * How state members are adjusted at the start of a call to `run` should be controlled via the implementation of - * `manage_state_before_next_time_step_run`. - * - * Calls to run will return the result of a nested call to `mass_check`, which performs mass balance verification on - * the calculated state at the end of a call to `run`. The amount of acceptable difference is accessible with the - * `get_mass_check_error_bound` function. It is also possible to change the value of this (hard-coded in the - * default implementation) with the protected `set_mass_check_error_bound` function. - * - */ - class tshirt_model { - - public: - - /** - * Constructor for model object based on model parameters and initial state. - * - * @param model_params Model parameters tshirt_params struct. - * @param initial_state Shared smart pointer to tshirt_state struct hold initial state values - */ - tshirt_model(tshirt_params model_params, const shared_ptr& initial_state); - - /** - * Constructor for model object with parameters only. - * - * Constructor creates a "default" initial state, with soil_storage_meters and groundwater_storage_meters set to - * 0.0, and then otherwise proceeds as the constructor accepting the second parameter for initial state. - * - * @param model_params Model parameters tshirt_params struct. - */ - tshirt_model(tshirt_params model_params); - - /** - * Calculate losses due to evapotranspiration. - * - * @param soil_m - * @param et_params - * @return - */ - double calc_evapotranspiration(double soil_m, shared_ptr et_params); - - /** - * Calculate soil field capacity storage threshold, the level at which free drainage stops (i.e., "Sfc"). - * - * @return The calculated soil field capacity storage. - */ - double calc_soil_field_capacity_storage_threshold(); - - /** - * Return the smart pointer to the tshirt::tshirt_model struct for holding this object's current state. - * - * @return The smart pointer to the tshirt_model struct for holding this object's current state. - */ - shared_ptr get_current_state(); - - /** - * Return the shared pointer to the tshirt::tshirt_fluxes struct for holding this object's current fluxes. - * - * @return The shared pointer to the tshirt_fluxes struct for holding this object's current fluxes. - */ - shared_ptr get_fluxes(); - - double get_mass_check_error_bound(); - - /** - * Check that mass was conserved by the model's calculations of the current time step. - * - * @param input_storage_m The amount of water input to the system at the time step, in meters. - * @param timestep_s The size of the time step, in seconds. - * @return The appropriate code value indicating whether mass was conserved in the current time step's calculations. - */ - int mass_check(double input_storage_m, double timestep_s); - - /** - * Run the model to one time step, after performing initial housekeeping steps via a call to - * `manage_state_before_next_time_step_run`. - * - * @param dt the time step - * @param input_storage_m the amount water entering the system this time step, in meters - * @param et_params ET parameters struct - * @return - */ - int run(double dt, double input_storage_m, shared_ptr et_params); - - protected: - - /** - * Perform necessary steps prior to the execution of model calculations for a new time step, for managing member - * variables that contain model state. - * - * This function is intended to be run only at the start of a new execution of the tshirt_model::run method. It - * performs three housekeeping tasks needed before running the next group of time step modeling operations: - * - * * the initial maintained `current_state` is moved to `previous_state` - * * a new `current_state` is created - * * a new `fluxes` is created - */ - virtual void manage_state_before_next_time_step_run(); - - /** - * Set the mass_check_error_bound member to the absolute value of the given parameter. - * - * @param error_bound The value used to set the mass_check_error_bound member. - */ - void set_mass_check_error_bound(double error_bound); - - private: - /** Model state for the "current" time step, which may not be calculated yet. */ - shared_ptr current_state; - /** Model execution parameters. */ - tshirt_params model_params; - /** Model state from that previous time step before the current. */ - shared_ptr previous_state; - /** - * A collection of reservoirs for a Nash Cascade at the end of the lateral flow output from the subsurface soil - * reservoir. - */ - vector> soil_lf_nash_res; - //FIXME reservoir construction sorts outlets by activation_threshold - //so the fixed index assumption is invalid. However, in the current use case - //they both have the save activation_threshold (Sfc), but we do want percolation fluxes to happen first - //so make it index 0 - /** The index of the subsurface lateral flow outlet in the soil reservoir. */ - int lf_outlet_index = 1; - /** The index of the percolation flow outlet in the soil reservoir. */ - int perc_outlet_index = 0; - Reservoir::Explicit_Time::Reservoir soil_reservoir; - Reservoir::Explicit_Time::Reservoir groundwater_reservoir; - shared_ptr fluxes; - /** The size of the error bound that is acceptable when performing mass check calculations. */ - double mass_check_error_bound; - /** Soil field capacity storage threshold, or the level at which free drainage stops (i.e., "Sfc"). */ - double soil_field_capacity_storage_threshold; - - /** - * Check that the current state of this model object (which could be its provided initial state) is valid, printing - * a message and raising an error if not. - * - * The model parameter for Nash Cascade size, `nash_n`, must correspond appropriately to the size of referenced - * `current_state->nash_cascade_storeage_meters` vector. The latter holds the storage values of the individual - * reservoirs within the Nash Cascade. Note that the function will interpret any `nash_n` greater than `0` as valid - * if the vector itself is empty, and initialize such a vector to the correct size with all `0.0` values. - */ - void check_valid(); - - /** - * Initialize the subsurface groundwater reservoir for the model, in the `groundwater_reservoir` member field. - * - * Initialize the subsurface groundwater reservoir for the model as a Reservoir object, creating the - * reservoir with a single outlet. In particular, this is a Reservoir_Exponential_Outlet object, since the outlet - * requires the following be used to calculate discharge flow: - * - * Cgw * ( exp(expon * S / S_max) - 1 ); - * - * Note that this function should only be used during object construction. - * - * @see Reservoir - * @see Reservoir_Exponential_Outlet - */ - void initialize_groundwater_reservoir(); - - /** - * Initialize the subsurface soil reservoir for the model, in the `soil_reservoir` member field. - * - * Initialize the subsurface soil reservoir for the model as a Reservoir object, creating the reservoir - * with outlets for both the subsurface lateral flow and the percolation flow. This should only be used during - * object construction. - * - * Per the class type of the reservoir, outlets have an associated index value within a reservoir, and certain - * outlet-specific functionality requires having appropriate outlet index. These outlet indexes are maintained in - * for the lateral flow and percolation flow outlets are maintained this class within the lf_outlet_index and - * perc_outlet_index member variables respectively. - * - * @see Reservoir - */ - void initialize_soil_reservoir(); - - /** - * Initialize the Nash Cascade reservoirs applied to the subsurface soil reservoir's lateral flow outlet. - * - * Initialize the soil_lf_nash_res member, containing the collection of Reservoir objects used to create - * the Nash Cascade for soil_reservoir lateral flow outlet. The analogous values for Nash Cascade storage from - * previous_state are used for current storage of reservoirs at each given index. - */ - void initialize_subsurface_lateral_flow_nash_cascade(); - - }; -} - -#endif //TSHIRT_H diff --git a/models/tshirt/include/TshirtErrorCodes.h b/models/tshirt/include/TshirtErrorCodes.h deleted file mode 100644 index 6ce23ede79..0000000000 --- a/models/tshirt/include/TshirtErrorCodes.h +++ /dev/null @@ -1,17 +0,0 @@ -#ifndef NGEN_TSHIRTERRORCODES_H -#define NGEN_TSHIRTERRORCODES_H - -namespace tshirt { - - // TODO: consider combining with or differentiating from similar hymod enum - - /** - * An enumeration of codes indicating possible error states (or lack thereof) when running the Tshirt model. - */ - enum TshirtErrorCodes { - TSHIRT_NO_ERROR = 0, - TSHIRT_MASS_BALANCE_ERROR = 100 - }; -} - -#endif //NGEN_TSHIRTERRORCODES_H diff --git a/models/tshirt/include/tshirt_c.h b/models/tshirt/include/tshirt_c.h deleted file mode 100644 index 578a6be6c6..0000000000 --- a/models/tshirt/include/tshirt_c.h +++ /dev/null @@ -1,151 +0,0 @@ -#ifndef NGEN_TSHIRT_C_HPP -#define NGEN_TSHIRT_C_HPP - -#define TRUE 1 -#define FALSE 0 -#define DEBUG 1 -#define MAX_NUM_GIUH_ORDINATES 10 -#define MAX_NUM_NASH_CASCADE 3 -#define TSHIRT_C_FIXED_TIMESTEP_SIZE_S 3600 -#define MAX_NUM_RAIN_DATA 720 - -// define data structures -//-------------------------- - -struct conceptual_reservoir -{ -// this data structure describes a nonlinear reservoir having two outlets, one primary with an activation -// threshold that may be zero, and a secondary outlet with a threshold that may be zero -// this will also simulate a linear reservoir by setting the exponent parameter to 1.0 iff is_exponential==FALSE -// iff is_exponential==TRUE, then it uses the exponential discharge function from the NWM V2.0 forumulation -// as the primary discharge with a zero threshold, and does not calculate a secondary discharge. -//-------------------------------------------------------------------------------------------------- - // TODO: adjust this (and usages) to bool type once verfied and automated testing is in place. - int is_exponential; // set this true TRUE to use the exponential form of the discharge equation - double storage_max_m; // maximum storage in this reservoir - double storage_m; // state variable. - double coeff_primary; // the primary outlet - double exponent_primary; - double storage_threshold_primary_m; - double storage_threshold_secondary_m; - double coeff_secondary; - double exponent_secondary; -}; - -struct NWM_soil_parameters -{ -// using same variable names as used in NWM. - double smcmax; // effective porosity [V/V] - double wltsmc; // wilting point soil moisture content [V/V] - double satdk; // saturated hydraulic conductivity [m s-1] - double satpsi; // saturated capillary head [m] - double bb; // beta exponent on Clapp-Hornberger (1978) soil water relations [-] - double mult; // the multiplier applied to satdk to route water rapidly downslope - double slop; // this factor (0-1) modifies the gradient of the hydraulic head at the soil bottom. 0=no-flow. - double D; // soil depth [m] -}; - -//DATA STRUCTURE TO HOLD AORC FORCING DATA -struct aorc_forcing_data -{ -// struct NAME DESCRIPTION ORIGINAL AORC NAME -//____________________________________________________________________________________________________________________ - float precip_kg_per_m2; // Surface precipitation "kg/m^2" | APCP_surface - float incoming_longwave_W_per_m2 ; // Downward Long-Wave Rad. Flux at 0m height, W/m^2 | DLWRF_surface - float incoming_shortwave_W_per_m2; // Downward Short-Wave Radiation Flux at 0m height, W/m^2 | DSWRF_surface - float surface_pressure_Pa; // Surface atmospheric pressure, Pa | PRES_surface - float specific_humidity_2m_kg_per_kg; // Specific Humidity at 2m height, kg/kg | SPFH_2maboveground - float air_temperature_2m_K; // Air temparture at 2m height, K | TMP_2maboveground - float u_wind_speed_10m_m_per_s; // U-component of Wind at 10m height, m/s | UGRD_10maboveground - float v_wind_speed_10m_m_per_s; // V-component of Wind at 10m height, m/s | VGRD_10maboveground - float latitude; // degrees north of the equator. Negative south | latitude - float longitude; // degrees east of prime meridian. Negative west | longitude - long int time; //TODO: type? // seconds since 1970-01-01 00:00:00.0 0:00 | time -} ; - -struct tshirt_c_result_fluxes { - double timestep_rainfall_input_m; - double Schaake_output_runoff_m; - double giuh_runoff_m; - double nash_lateral_runoff_m; - double flux_from_deep_gw_to_chan_m; - double Qout_m; -}; - -// TODO: add converter helper functions for going between aorc_forcing_data and AORC_data from Forcing.h - -// function prototypes -// -------------------------------- -extern void Schaake_partitioning_scheme(double dt, double magic_number, double deficit, double qinsur, - double *runsrf, double *pddum); - -extern void conceptual_reservoir_flux_calc(struct conceptual_reservoir *da_reservoir, - double *primary_flux,double *secondary_flux); - -extern double convolution_integral(double runoff_m, int num_giuh_ordinates, - double *giuh_ordinates, double *runoff_queue_m_per_timestep); - -extern double nash_cascade(double flux_lat_m, int num_lateral_flow_nash_reservoirs, - double K_nash, double *nash_storage); - -extern int is_fabs_less_than_epsilon(double a,double epsilon); // returns TRUE iff fabs(a) - -namespace tshirt { - - /** - * Structure for containing and organizing the parameters of the Tshirt hydrological model. - */ - struct tshirt_params { - double maxsmc; //!< saturated soil moisture content (sometimes theta_e or smcmax) - double wltsmc; //!< wilting point soil moisture content - double satdk; //!< vertical saturated hydraulic conductivity [m s^-1] (sometimes Kperc or Ks) - double satpsi; //!< saturated capillary head [m] - // TODO: explain more what this is - double slope; //!< SLOPE parameter - double b; //!< 'b' exponent on Clapp-Hornberger soil water relations (sometime bexp) - double multiplier; //!< the multiplier applied to 'satdk' to route water rapidly downslope in subsurface (sometimes 'mult' or 'LKSATFAC') - double alpha_fc; //!< alpha constant for given soil type for relative suction head value, with respect to Hatm - double Klf; //!< lateral flow independent calibration parameter - double Kn; //!< Nash cascade linear reservoir coefficient lateral flow parameter - int nash_n; //!< number of nash cascades - double Cgw; //!< Ground water flow param - double Cschaake; //!< The Schaake adjusted magic constant by soil type - double expon; //!< Ground water flow exponent param (analogous to NWM 2.0 expon param) - double max_soil_storage_meters; //!< Subsurface soil water flow max storage param ("Ssmax"), calculated from maxsmc and depth - double max_groundwater_storage_meters; //!< Ground water flow max storage param ("Sgwmax"; analogous to NWM 2.0 zmax param) - double max_lateral_flow; //!< Max rate for subsurface lateral flow (i.e., max transmissivity) - double refkdt = 3.0; //!< Standard value taken from NOAH-MP, used for calculating Schaake magic constant - const double depth = 2.0; //!< Hard-coded, constant value for total soil column depth ('D') [m] - - /** - * Constructor for instances, initializing members that correspond one-to-one with a parameter and deriving - * remaining (non-const and hard-coded) members. - * - * @param maxsmc Initialization param for maxsmc member. - * @param wltsmc Initialization param for wltsmc member. - * @param satdk Initialization param for satdk member. - * @param satpsi Initialization param for satpsi member. - * @param slope Initialization param for slope member. - * @param b Initialization param for b member. - * @param multiplier Initialization param for multiplier member. - * @param alpha_fc Initialization param for alpha_fc member. - * @param Klf Initialization param for Klf member. - * @param Kn Initialization param for Kn member. - * @param nash_n Initialization param for nash_n member. - * @param Cgw Initialization param for Cgw member. - * @param expon Initialization param for expon member. - * @param max_gw_storage Initialization param for max_groundwater_storage_meters member. - */ - tshirt_params( - double maxsmc, - double wltsmc, - double satdk, - double satpsi, - double slope, - double b, - double multiplier, - double alpha_fc, - double Klf, - double Kn, - int nash_n, - double Cgw, - double expon, - double max_gw_storage) : - maxsmc(maxsmc), - wltsmc(wltsmc), - satdk(satdk), - satpsi(satpsi), - slope(slope), - b(b), - multiplier(multiplier), - alpha_fc(alpha_fc), - Klf(Klf), - Kn(Kn), - nash_n(nash_n), - Cgw(Cgw), - expon(expon), - max_groundwater_storage_meters(max_gw_storage) { - this->max_soil_storage_meters = this->depth * maxsmc; - this->Cschaake = refkdt * satdk / (2.0e-6); - this->max_lateral_flow = std::numeric_limits::max();//satdk * multiplier * this->max_soil_storage_meters; - } - - }; -} - -#endif //NGEN_TSHIRT_PARAMS_H diff --git a/models/tshirt/include/tshirt_state.h b/models/tshirt/include/tshirt_state.h deleted file mode 100644 index 8853445200..0000000000 --- a/models/tshirt/include/tshirt_state.h +++ /dev/null @@ -1,26 +0,0 @@ -#ifndef NGEN_TSHIRT_STATE_H -#define NGEN_TSHIRT_STATE_H - -namespace tshirt { - - /** - * Structure for containing and organizing the state values of the Tshirt hydrological model. - */ - struct tshirt_state { - // TODO: confirm this is correct - double soil_storage_meters; //!< current water storage in soil column nonlinear reservoir ("Ss") - double groundwater_storage_meters; //!< current water storage in ground water nonlinear reservoir ("Sgw") - vector nash_cascade_storeage_meters; //!< water storage in nonlinear reservoirs of Nash Cascade for lateral subsurface flow - - // I think this doesn't belong in state, and so is just in run() in the tshirt::tshirt_model class - //double column_total_soil_moisture_deficit; //!< soil column total moisture deficit - - tshirt_state(double soil_storage_meters = 0.0, double groundwater_storage_meters = 0.0, - vector nash_cascade_storeage_meters = vector()) - : soil_storage_meters(soil_storage_meters), - groundwater_storage_meters(groundwater_storage_meters), - nash_cascade_storeage_meters(std::move(nash_cascade_storeage_meters)) {} - }; -} - -#endif //NGEN_TSHIRT_STATE_H diff --git a/src/NGen.cpp b/src/NGen.cpp index 56d1bdb607..44cf5ce28d 100644 --- a/src/NGen.cpp +++ b/src/NGen.cpp @@ -14,7 +14,6 @@ #endif #include "NGenConfig.h" -#include "tshirt_params.h" #include #include diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index cd0c8fee62..67df894659 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -18,7 +18,6 @@ target_include_directories(core PUBLIC ${PROJECT_SOURCE_DIR}/models/kernels ${PROJECT_SOURCE_DIR}/models/kernels/evapotranspiration ${PROJECT_SOURCE_DIR}/models/ - ${PROJECT_SOURCE_DIR}/models/tshirt/include ${PROJECT_SOURCE_DIR}/include/geojson ${PROJECT_SOURCE_DIR}/include/forcing ${PROJECT_SOURCE_DIR}/include/utilities @@ -39,7 +38,6 @@ if (${NGEN_ACTIVATE_PYTHON}) ${PROJECT_SOURCE_DIR}/models/kernels ${PROJECT_SOURCE_DIR}/models/kernels/evapotranspiration ${PROJECT_SOURCE_DIR}/models/ - ${PROJECT_SOURCE_DIR}/models/tshirt/include ${PROJECT_SOURCE_DIR}/include/geojson ${PROJECT_SOURCE_DIR}/include/forcing ${PROJECT_SOURCE_DIR}/include/utilities @@ -65,7 +63,6 @@ else () ${PROJECT_SOURCE_DIR}/models/kernels ${PROJECT_SOURCE_DIR}/models/kernels/evapotranspiration ${PROJECT_SOURCE_DIR}/models/ - ${PROJECT_SOURCE_DIR}/models/tshirt/include ${PROJECT_SOURCE_DIR}/include/geojson ${PROJECT_SOURCE_DIR}/include/forcing ${PROJECT_SOURCE_DIR}/include/utilities diff --git a/src/models/tshirt/CMakeLists.txt b/src/models/tshirt/CMakeLists.txt deleted file mode 100644 index 9cbf2ece49..0000000000 --- a/src/models/tshirt/CMakeLists.txt +++ /dev/null @@ -1,15 +0,0 @@ -cmake_minimum_required(VERSION 3.10) -add_library(models_tshirt STATIC - Tshirt.cpp - tshirt_c.cpp) -add_library(NGen::models_tshirt ALIAS models_tshirt) -target_include_directories(models_tshirt PUBLIC - ${PROJECT_SOURCE_DIR}/include - ${PROJECT_SOURCE_DIR}/include/core/catchment/giuh - ${PROJECT_SOURCE_DIR}/include/utilities - ${PROJECT_SOURCE_DIR}/models/kernels - ${PROJECT_SOURCE_DIR}/models/tshirt/include - ) - -target_link_libraries(models_tshirt PUBLIC - NGen::kernels_reservoir) diff --git a/src/models/tshirt/Tshirt.cpp b/src/models/tshirt/Tshirt.cpp deleted file mode 100644 index 70b594b51f..0000000000 --- a/src/models/tshirt/Tshirt.cpp +++ /dev/null @@ -1,369 +0,0 @@ -#include "Tshirt.h" -#include "TshirtErrorCodes.h" -#include "tshirt_fluxes.h" -#include "tshirt_state.h" -#include "Constants.h" -#include - -namespace tshirt { - - /** - * Constructor for model object based on model parameters and initial state. - * - * @param model_params Model parameters tshirt_params struct. - * @param initial_state Shared smart pointer to tshirt_state struct hold initial state values - */ - tshirt_model::tshirt_model(tshirt_params model_params, const shared_ptr &initial_state) - : model_params(model_params), previous_state(initial_state), current_state(initial_state) - { - // ********** Start by calculating Sfc, as this will get by several other things - soil_field_capacity_storage_threshold = calc_soil_field_capacity_storage_threshold(); - - // ********** Sanity check init (in particular, size of Nash Cascade storage vector in the state parameter). - check_valid(); - - // ********** Create the vector of Nash Cascade reservoirs used at the end of the soil lateral flow outlet - initialize_subsurface_lateral_flow_nash_cascade(); - - // ********** Create the soil reservoir - initialize_soil_reservoir(); - - // ********** Create the groundwater reservoir - initialize_groundwater_reservoir(); - - // ********** Set fluxes to null for now: it is bogus until first call of run function, which initializes it - fluxes = nullptr; - - // ********** Acceptable error range for mass balance calculations; hard-coded for now to this value - mass_check_error_bound = 0.000001; - } - - /** - * Constructor for model object with parameters only. - * - * Constructor creates a "default" initial state, with soil_storage_meters and groundwater_storage_meters set to - * 0.0, and then otherwise proceeds as the constructor accepting the second parameter for initial state. - * - * @param model_params Model parameters tshirt_params struct. - */ - tshirt_model::tshirt_model(tshirt_params model_params) : - tshirt_model(model_params, make_shared(tshirt_state(0.0, 0.0))) {} - - /** - * Check that the current state of this model object (which could be its provided initial state) is valid, printing - * a message and raising an error if not. - * - * The model parameter for Nash Cascade size, `nash_n`, must correspond appropriately to the size of referenced - * `current_state->nash_cascade_storeage_meters` vector. The latter holds the storage values of the individual - * reservoirs within the Nash Cascade. Note that the function will interpret any `nash_n` greater than `0` as valid - * if the vector itself is empty, and initialize such a vector to the correct size with all `0.0` values. - */ - void tshirt_model::check_valid() - { - // We expect the Nash size model parameter 'nash_n' to be equal to the size of the - // 'nash_cascade_storeage_meters' member of the passed 'initial_state param (to which 'current_state' is set) - if (current_state->nash_cascade_storeage_meters.size() != model_params.nash_n) { - // Infer that empty vector should be initialized to a vector of size 'nash_n' with all 0.0 values - if (model_params.nash_n > 0 && current_state->nash_cascade_storeage_meters.empty()) { - current_state->nash_cascade_storeage_meters.resize(model_params.nash_n); - for (unsigned i = 0; i < model_params.nash_n; ++i) { - current_state->nash_cascade_storeage_meters[i] = 0.0; - } - } - else { - cerr << "ERROR: Nash Cascade size parameter in tshirt model init doesn't match storage vector size " - << "in state parameter" - << endl; - // TODO: return error of some kind here - } - } - } - - /** - * Initialize the subsurface groundwater reservoir for the model, in the `groundwater_reservoir` member field. - * - * Initialize the subsurface groundwater reservoir for the model as a Reservoir object, creating the - * reservoir with a single outlet. In particular, this is a Reservoir_Exponential_Outlet object, since the outlet - * requires the following be used to calculate discharge flow: - * - * Cgw * ( exp(expon * S / S_max) - 1 ); - * - * Note that this function should only be used during object construction. - * - * @see Reservoir - * @see Reservoir_Exponential_Outlet - */ - void tshirt_model::initialize_groundwater_reservoir() - { - // TODO: confirm, based on the equation, that max gw velocity doesn't need to be Cgw * (exp(expon) - 1) - // TODO: (i.e., S == S_max, thus maximizing the term passed to the exp() function, and thereby the equation) - double max_gw_velocity = std::numeric_limits::max(); - - // Build vector of pointers to outlets to pass the custom exponential outlet through - vector> gw_outlets_vector(1); - // TODO: verify activation threshold - gw_outlets_vector[0] = make_shared( - Reservoir::Explicit_Time::Reservoir_Exponential_Outlet(model_params.Cgw, model_params.expon, 0.0, max_gw_velocity)); - // Create the reservoir, passing the outlet via the vector argument - groundwater_reservoir = Reservoir::Explicit_Time::Reservoir(0.0, model_params.max_groundwater_storage_meters, - previous_state->groundwater_storage_meters, gw_outlets_vector); - } - - /** - * Initialize the subsurface soil reservoir for the model, in the `soil_reservoir` member field. - * - * Initialize the subsurface soil reservoir for the model as a Reservoir object, creating the reservoir - * with outlets for both the subsurface lateral flow and the percolation flow. This should only be used during - * object construction. - * - * Per the class type of the reservoir, outlets have an associated index value within a reservoir, and certain - * outlet-specific functionality requires having appropriate outlet index. These outlet indexes are maintained in - * for the lateral flow and percolation flow outlets are maintained this class within the lf_outlet_index and - * perc_outlet_index member variables respectively. - * - * @see Reservoir - */ - void tshirt_model::initialize_soil_reservoir() - { - // Build the vector of pointers to reservoir outlets - vector> soil_res_outlets(2); - - // init subsurface lateral flow linear outlet - soil_res_outlets[lf_outlet_index] = std::make_shared( - Reservoir::Explicit_Time::Reservoir_Linear_Outlet(model_params.Klf, soil_field_capacity_storage_threshold, model_params.max_lateral_flow)); - - // init subsurface percolation flow linear outlet - // The max perc flow should be equal to the params.satdk value - soil_res_outlets[perc_outlet_index] = std::make_shared( - Reservoir::Explicit_Time::Reservoir_Linear_Outlet(model_params.satdk * model_params.slope, soil_field_capacity_storage_threshold, - std::numeric_limits::max())); - // Create the reservoir, included the created vector of outlet pointers - soil_reservoir = Reservoir::Explicit_Time::Reservoir(0.0, model_params.max_soil_storage_meters, - previous_state->soil_storage_meters, soil_res_outlets); - } - - /** - * Initialize the Nash Cascade reservoirs applied to the subsurface soil reservoir's lateral flow outlet. - * - * Initialize the soil_lf_nash_res member, containing the collection of Reservoir objects used to create - * the Nash Cascade for soil_reservoir lateral flow outlet. The analogous values for Nash Cascade storage from - * previous_state are used for current storage of reservoirs at each given index. - */ - void tshirt_model::initialize_subsurface_lateral_flow_nash_cascade() - { - soil_lf_nash_res.resize(model_params.nash_n); - // TODO: verify correctness of activation_threshold (Sfc) and max_velocity (max_lateral_flow) arg values - for (unsigned long i = 0; i < soil_lf_nash_res.size(); ++i) { - //construct a single linear outlet reservoir - soil_lf_nash_res[i] = make_unique( - Reservoir::Explicit_Time::Reservoir(0.0, model_params.max_soil_storage_meters, - previous_state->nash_cascade_storeage_meters[i], model_params.Kn, - 0.0, model_params.max_lateral_flow)); - } - } - - /** - * Calculate losses due to evapotranspiration. - * - * @param soil_m The soil moisture measured in meters. - * @param et_params A shared pointer to the struct holding the ET parameters. - * @return The calculated loss value due to evapotranspiration. - */ - double tshirt_model::calc_evapotranspiration(double soil_m, shared_ptr et_params) { - et_params->final_height_reservoir = soil_m; - pdm03_wrapper(et_params.get()); - - return et_params->final_height_reservoir - soil_m; - } - - /** - * Calculate soil field capacity storage threshold, the level at which free drainage stops (i.e., "Sfc"). - * - * @return The calculated soil field capacity storage. - */ - double tshirt_model::calc_soil_field_capacity_storage_threshold() - { - // Calculate the suction head above water table (Hwt) - double head_above_water_table = - model_params.alpha_fc * ((double) STANDARD_ATMOSPHERIC_PRESSURE_PASCALS / WATER_SPECIFIC_WEIGHT); - // TODO: account for possibility of Hwt being less than 0.5 (though initially, it looks like this will never be the case) - - // distance from the bottom of the soil column to the center of the lowest discretization - double trigger_z_m = 0.5; - - // Previous way this was calculated, but which seems to disagree with Fred's implementation ... - //double z1 = head_above_water_table - trigger_z_m; - //double z2 = z1 + 2; - //double lower_lim = model_params.b * pow(z1, ((model_params.b - 1) / model_params.b)) / (model_params.b - 1); - //double upper_lim = model_params.b * pow(z2, ((model_params.b - 1) / model_params.b)) / (model_params.b - 1); - - // ... so changing to this to reflect Fred's code - double omega = head_above_water_table - trigger_z_m; - // Use this several times - double one_minus_one_over_b = (1.0 - (1.0 / model_params.b)); - double lower_lim = pow(omega, one_minus_one_over_b) / one_minus_one_over_b; - double upper_lim = pow(omega + model_params.depth, one_minus_one_over_b) / one_minus_one_over_b; - - double lim_diff = upper_lim - lower_lim; - - // Note that z^( 1 - (1/b) ) / (1 - (1/b)) == b * (z^( (b - 1) / b ) / (b - 1) - return model_params.maxsmc * pow((1.0 / model_params.satpsi), (-1.0 / model_params.b)) * lim_diff; - } - - /** - * Return the smart pointer to the tshirt::tshirt_model struct for holding this object's current state. - * - * @return The smart pointer to the tshirt_model struct for holding this object's current state. - */ - shared_ptr tshirt_model::get_current_state() { - return current_state; - } - - /** - * Return the shared pointer to the tshirt::tshirt_fluxes struct for holding this object's current fluxes. - * - * @return The shared pointer to the tshirt_fluxes struct for holding this object's current fluxes. - */ - shared_ptr tshirt_model::get_fluxes() { - return fluxes; - } - - /** - * Get the size of the error bound that is acceptable when performing mass check calculations. - * - * @return The size of the error bound that is acceptable when performing mass check calculations. - * @see tshirt_model::mass_check - */ - double tshirt_model::get_mass_check_error_bound() { - return mass_check_error_bound; - } - - /** - * Check that mass was conserved by the model's calculations of the current time step. - * - * @param input_storage_m The amount of water input to the system at the time step, in meters. - * @param timestep_s The size of the time step, in seconds. - * @return The appropriate code value indicating whether mass was conserved in the current time step's calculations. - */ - int tshirt_model::mass_check(double input_storage_m, double timestep_s) { - // TODO: change this to have those be part of state somehow, either of object or struct, or just make private - // Initialize both mass values from current and next states storage - double previous_storage_m = previous_state->soil_storage_meters + previous_state->groundwater_storage_meters; - double current_storage_m = current_state->soil_storage_meters + current_state->groundwater_storage_meters; - - // Add the masses of the Nash reservoirs before and after the time step - for (unsigned int i = 0; i < previous_state->nash_cascade_storeage_meters.size(); ++i) { - previous_storage_m += previous_state->nash_cascade_storeage_meters[i]; - current_storage_m += current_state->nash_cascade_storeage_meters[i]; - } - - // Increase the initial mass by input value - previous_storage_m += input_storage_m; - - // Increase final mass by calculated fluxes that leave the system (i.e., not the percolation flow) - current_storage_m += fluxes->et_loss_meters; - current_storage_m += fluxes->surface_runoff_meters_per_second * timestep_s; - current_storage_m += fluxes->soil_lateral_flow_meters_per_second * timestep_s; - current_storage_m += fluxes->groundwater_flow_meters_per_second * timestep_s; - - double abs_mass_diff_meters = abs(previous_storage_m - current_storage_m); - return abs_mass_diff_meters > get_mass_check_error_bound() ? tshirt::TSHIRT_MASS_BALANCE_ERROR - : tshirt::TSHIRT_NO_ERROR; - } - - /** - * Run the model to one time step, after performing initial housekeeping steps via a call to - * `manage_state_before_next_time_step_run`. - * - * @param dt the time step size in seconds - * @param input_storage_m the amount water entering the system this time step, in meters - * @return - */ - int tshirt_model::run(double dt, double input_storage_m, shared_ptr et_params) { - // Do resetting/housekeeping for new calculations and new state values - manage_state_before_next_time_step_run(); - - // In meters - double soil_column_moisture_deficit_m = - model_params.max_soil_storage_meters - previous_state->soil_storage_meters; - - // Perform Schaake partitioning, passing some declared references to hold the calculated values. - double surface_runoff, subsurface_infiltration_flux; - Schaake_partitioning_scheme_cpp(dt, model_params.Cschaake, soil_column_moisture_deficit_m, input_storage_m, - &surface_runoff, &subsurface_infiltration_flux); - - double subsurface_excess, nash_subsurface_excess; - double mean_timestep_infiltration_m_per_s = subsurface_infiltration_flux / dt; - soil_reservoir.response_meters_per_second(mean_timestep_infiltration_m_per_s, (int)dt, subsurface_excess); - - // lateral subsurface flow - double Qlf = soil_reservoir.velocity_meters_per_second_for_outlet(lf_outlet_index); - - // percolation flow - double Qperc = soil_reservoir.velocity_meters_per_second_for_outlet(perc_outlet_index); - - // TODO: make sure ET doesn't need to be taken out sooner - // Get new soil storage amount calculated by reservoir - double new_soil_storage_m = soil_reservoir.get_storage_height_meters(); - // Calculate and store ET - fluxes->et_loss_meters = calc_evapotranspiration(new_soil_storage_m, et_params); - // Update the current soil storage, accounting for ET - current_state->soil_storage_meters = new_soil_storage_m - fluxes->et_loss_meters; - - // Cycle through lateral flow Nash cascade of reservoirs - // loop essentially copied from Hymod logic, but with different variable names - for (unsigned long int i = 0; i < soil_lf_nash_res.size(); ++i) { - // get response water velocity of reservoir - Qlf = soil_lf_nash_res[i]->response_meters_per_second(Qlf, dt, nash_subsurface_excess); - // TODO: confirm this is correct - Qlf += nash_subsurface_excess / dt; - current_state->nash_cascade_storeage_meters[i] = soil_lf_nash_res[i]->get_storage_height_meters(); - } - - // Get response and update gw res state - double excess_gw_water; - fluxes->groundwater_flow_meters_per_second = groundwater_reservoir.response_meters_per_second(Qperc, dt, - excess_gw_water); - // update local copy of state - current_state->groundwater_storage_meters = groundwater_reservoir.get_storage_height_meters(); - - // record other fluxes in internal copy - fluxes->soil_lateral_flow_meters_per_second = Qlf; - fluxes->soil_percolation_flow_meters_per_second = Qperc; - - // Save "raw" runoff here and have realization class calculate GIUH surface runoff using that kernel - // TODO: for now add this to runoff, but later adjust calculations to limit flow into reservoir to avoid excess - fluxes->surface_runoff_meters_per_second = surface_runoff + (subsurface_excess / dt) + (excess_gw_water / dt); - //fluxes->surface_runoff_meters_per_second = surface_runoff_depth_m; - - return mass_check(input_storage_m, dt); - } - - /** - * Perform necessary steps prior to the execution of model calculations for a new time step, for managing member - * variables that contain model state. - * - * This function is intended to be run only at the start of a new execution of the tshirt_model::run method. It - * performs three housekeeping tasks needed before running the next group of time step modeling operations: - * - * * the initial maintained `current_state` is moved to `previous_state` - * * a new `current_state` is created - * * a new `fluxes` is created - */ - void tshirt_model::manage_state_before_next_time_step_run() - { - previous_state = current_state; - current_state = make_shared(tshirt_state(0.0, 0.0, vector(model_params.nash_n))); - fluxes = make_shared(tshirt_fluxes(0.0, 0.0, 0.0, 0.0, 0.0)); - } - - /** - * Set the mass_check_error_bound member to the absolute value of the given parameter. - * - * @param error_bound The value used to set the mass_check_error_bound member. - */ - void tshirt_model::set_mass_check_error_bound(double error_bound) { - mass_check_error_bound = error_bound >= 0 ? error_bound : abs(error_bound); - } - -} diff --git a/src/models/tshirt/tshirt_c.cpp b/src/models/tshirt/tshirt_c.cpp deleted file mode 100644 index bb3a3a7bf4..0000000000 --- a/src/models/tshirt/tshirt_c.cpp +++ /dev/null @@ -1,616 +0,0 @@ -#include "tshirt_c.h" -#include "Constants.h" -#include -#include -#include -#include - -/** - * Solve for the flow through a Nash Cascade to delay the arrival of some flow. - * - * @param flux_lat_m Input to the cascade in meters. - * @param num_lateral_flow_nash_reservoirs The number of internal reservoirs for the Nash Cascade. - * @param K_nash The 'K' parameter for the cascade reservoirs. - * @param nash_storage Pointer to head of storage array for cascade reservoirs. - * @return The outflow from the cascade, in meters. - */ - // TODO: update with more agnostic parameter names once fully tested and testable in the future. -extern double nash_cascade(double flux_lat_m, int num_lateral_flow_nash_reservoirs, - double K_nash, double *nash_storage) { - // TODO: update this with more modern features once fully tested - int i; - double outflow_m; - static double Q[MAX_NUM_NASH_CASCADE]; - - // Loop through reservoirs - for (i = 0; i < num_lateral_flow_nash_reservoirs; i++) { - Q[i] = K_nash * nash_storage[i]; - nash_storage[i] -= Q[i]; - - if (i == 0) nash_storage[i] += flux_lat_m; - else nash_storage[i] += Q[i - 1]; - } - - /* Get Qout */ - outflow_m = Q[num_lateral_flow_nash_reservoirs - 1]; - - // Return the flow output - return (outflow_m); - -} - -/** - * Perform GIUH convolution integral solver calculation for an array of ordinates. - * - * @param runoff_m Input runoff in meters. - * @param num_giuh_ordinates The number of GIUH ordinates. - * @param giuh_ordinates Pointer to the head of an array of GUIH ordinates. - * @param runoff_queue_m_per_timestep Pointer to array of calculated convolution queue runoff values for the ordinates. - * @return The calculated amount for current GIUH runoff in meters. - */ -extern double convolution_integral(double runoff_m, int num_giuh_ordinates, - double *giuh_ordinates, double *runoff_queue_m_per_timestep) { - //############################################################## - // This function solves the convolution integral involving N - // GIUH ordinates. - //############################################################## - double runoff_m_now; - int N, i; - - N = num_giuh_ordinates; - runoff_queue_m_per_timestep[N] = 0.0; - - for (i = 0; i < N; i++) { - runoff_queue_m_per_timestep[i] += giuh_ordinates[i] * runoff_m; - } - runoff_m_now = runoff_queue_m_per_timestep[0]; - - // shift all the entries in preparation for the next time step - for (i = 0; i < N; i++) { - runoff_queue_m_per_timestep[i] = runoff_queue_m_per_timestep[i + 1]; - } - - return (runoff_m_now); -} - -/** - * Perform calculations for the flux for a conceptual reservoir with one or two outlet - * - * Perform calculations for the flux for a conceptual reservoir (linear or nonlinear) with one or two outlets, or an - * exponential nonlinear conceptual reservoir with a single outlet. In the case of non-exponential reservoir, each - * outlet can have its own activation storage threshold. - * - * Flow for the second outlet is turned off in the conceptual reservoir by setting the discharge coefficient to 0.0. - * - * All fluxes calculated by this routine are instantaneous with units of the coefficient. - * - * Flux values are returned by reference pointers passed in as arguments. - * - * @param reservoir The #conceptual_reservoir structure. - * @param primary_flux_m Pointer to location for holding calculated flux (in meters) for primary outlet. - * @param secondary_flux_m Pointer to location for holding calculated flux (in meters) for secondary outlet. - */ -extern void conceptual_reservoir_flux_calc(struct conceptual_reservoir *reservoir, - double *primary_flux_m, double *secondary_flux_m) -{ - /* - * <<<>>> - * - struct conceptual_reservoir - { - int is_exponential; // set this true TRUE to use the exponential form of the discharge equation - double storage_max_m; - double storage_m; - double coeff_primary; - double exponent_secondary; - double storage_threshold_primary_m; - double storage_threshold_secondary_m; - double coeff_secondary; - double exponent_secondary; - }; - - */ - - // local variables - double storage_above_threshold_m; - - if (reservoir->is_exponential == - TRUE) // single outlet reservoir like the NWM V1.2 exponential conceptual gw reservoir - { - // calculate the one flux and return. - *primary_flux_m = reservoir->coeff_primary * - (exp(reservoir->exponent_primary * reservoir->storage_m / reservoir->storage_max_m) - 1.0); - *secondary_flux_m = 0.0; - return; - } - // code goes past here iff it is not a single outlet exponential deep groundwater reservoir of the NWM variety - // The vertical outlet is assumed to be primary and satisfied first. - - *primary_flux_m = 0.0; - storage_above_threshold_m = reservoir->storage_m - reservoir->storage_threshold_primary_m; - if (storage_above_threshold_m > 0.0) { - // flow is possible from the primary outlet - *primary_flux_m = reservoir->coeff_primary * - pow(storage_above_threshold_m / - (reservoir->storage_max_m - reservoir->storage_threshold_primary_m), - reservoir->exponent_primary); - if (*primary_flux_m > storage_above_threshold_m) - *primary_flux_m = storage_above_threshold_m; // limit to max. available - } - *secondary_flux_m = 0.0; - storage_above_threshold_m = reservoir->storage_m - reservoir->storage_threshold_secondary_m; - if (storage_above_threshold_m > 0.0) { - // flow is possible from the secondary outlet - *secondary_flux_m = reservoir->coeff_secondary * - pow(storage_above_threshold_m / - (reservoir->storage_max_m - reservoir->storage_threshold_secondary_m), - reservoir->exponent_secondary); - if (*secondary_flux_m > (storage_above_threshold_m - (*primary_flux_m))) - *secondary_flux_m = storage_above_threshold_m - (*primary_flux_m); // limit to max. available - } - return; -} - -/** - * Perform Schaake Runoff partitioning. - * - * This function takes water_input_depth_m and partitions it into surface_runoff_depth_m and infiltration_depth_m using - * the scheme from Schaake et al. 1996. The logic has been modified to eliminate reference to ice processes, and to - * and to de-obfuscate and use descriptive and dimensionally consistent variable names. - * - * Function "returns" the surface runoff and infiltration amounts via passed in reference pointer parameters. - * - * Note the folloing general formula for the Schaake magic constant: - * Schaake_adjusted_magic_constant_by_soil_type = C*Ks(soiltype)/Ks_ref, where C=3, and Ks_ref=2.0E-06 m/s - * - * @param timestep_h Time step size in hours - * @param Schaake_adjusted_magic_constant_by_soil_type Magic constant value for soil type. - * @param column_total_soil_moisture_deficit_m Column moisture deficit in meters. - * @param water_input_depth_m Amount of input water to soil surface this time step, in meters. - * @param surface_runoff_depth_m Pointer for partitioned surface runoff for this time step, in meters. - * @param infiltration_depth_m Pointer for partitioned infiltration depth for this time step, in meters. - */ -void Schaake_partitioning_scheme(double timestep_h, double Schaake_adjusted_magic_constant_by_soil_type, - double column_total_soil_moisture_deficit_m, double water_input_depth_m, - double *surface_runoff_depth_m, double *infiltration_depth_m) { - int k; - double timestep_d, Schaake_parenthetical_term, Ic, Px, infilt_dep_m; - - - if (0.0 < water_input_depth_m) { - if (0.0 > column_total_soil_moisture_deficit_m) { - *surface_runoff_depth_m = water_input_depth_m; - *infiltration_depth_m = 0.0; - } - else { - // partition time-step total applied water as per Schaake et al. 1996. - // change from dt in [s] to dt1 in [d] because kdt has units of [d^(-1)] - timestep_d = timestep_h / 24.0; // timestep_d is the time step in days. - - // calculate the parenthetical part of Eqn. 34 from Schaake et al. Note the magic constant has units of [d^(-1)] - - Schaake_parenthetical_term = (1.0 - exp(-Schaake_adjusted_magic_constant_by_soil_type * timestep_d)); - - // From Schaake et al. Eqn. 2., using the column total moisture deficit - // BUT the way it is used here, it is the cumulative soil moisture deficit in the entire soil profile. - // "Layer" info not used in this subroutine in noah-mp, except to sum up the total soil moisture storage. - // NOTE: when column_total_soil_moisture_deficit_m becomes zero, which occurs when the soil column is saturated, - // then Ic=0, where Ic in the Schaake paper is called the "spatially averaged infiltration capacity", - // and is defined in Eqn. 12. - - Ic = column_total_soil_moisture_deficit_m * Schaake_parenthetical_term; - - Px = water_input_depth_m; // Total water input to partitioning scheme this time step [m] - - // This is eqn 24 from Schaake et al. NOTE: this is 0 in the case of a saturated soil column, when Ic=0. - // Physically happens only if soil has no-flow lower b.c. - - *infiltration_depth_m = (Px * (Ic / (Px + Ic))); - - - if (0.0 < (water_input_depth_m - (*infiltration_depth_m))) { - *surface_runoff_depth_m = water_input_depth_m - (*infiltration_depth_m); - } else *surface_runoff_depth_m = 0.0; - - *infiltration_depth_m = water_input_depth_m - (*surface_runoff_depth_m); - } - } - else { - *surface_runoff_depth_m = 0.0; - *infiltration_depth_m = 0.0; - } - return; -} - -/** - * Test if ``fabs(a)`` is less than epsilon. - * - * @param a The fabs function argument value. - * @param epsilon The epsilon value. - * @return Whether ``fabs(a)`` is less than epsilon. - */ -// TODO: adjust to using bool type once verified and testing is in place -extern int is_fabs_less_than_epsilon(double a, double epsilon) { - if (fabs(a) < epsilon) return (TRUE); - else return (FALSE); -} - - // TODO: probably need to produce a "generic" parameters struct or object somehow - // TODO: convert from tshirt_params to NWM_soil_parameters struct via separate function in realization and test. -/** - * Run the model, recording fluxes in one or more flux structs at a specified location. - * - * @param NWM_soil_params Struct containing soil model params - * @param gw_reservoir Struct containing state for ground water reservoir - * @param soil_reservoir Struct containing state for soil reservoir - * @param num_timesteps The number of time steps, which PROBABLY should always be 1 - * @param giuh_ordinates A pointer to the head of an array of GIUH ordinate values - * @param num_giuh_ordinates The number of elements in the array pointed to by ``giuh_ordinates`` - * @param runoff_queue_m_per_timestep Pointer to head of array for GIUH convolution queue, of size ``num_giuh_ordinates`` + 1 - * @param field_capacity_atm_press_fraction Fraction of atmospheric pressure for field capacity suction pressure (also alpha_fc) - * @param water_table_slope Assumed near channel water table slope lateral flow parameter - * @param Schaake_adjusted_magic_constant_by_soil_type Schaake magic constant for soil type - * @param lateral_flow_linear_reservoir_constant Later flow reservoir constant (Klf) - * @param K_nash Nash reservoir flow constant - * @param num_lateral_flow_nash_reservoirs The number of reservoirs in the cascade for the lateral flow - * @param nash_storage Pointer to head of array of already-set Nash Cascade reservoir storage values - * @param yes_aorc Whether - * @param aorc_data Pointer to head of array of per-time-step AORC data structs, when yes_aorc is TRUE (array of size num_timesteps) - * @param rain_rate Pointer to head of array of per-time-step simple rain rate data (in m/h), when yes_aorc is FALSE (array of size num_timesteps) - * @param num_added_fluxes Reference for number of fluxes added to 'fluxes' parameter array; initialized to 0. - * @param fluxes Pointer to head of array of array per-time-step resulting flux values (array of size num_timesteps) - * @return - */ -extern int run(NWM_soil_parameters& NWM_soil_params, - conceptual_reservoir& gw_reservoir, - conceptual_reservoir& soil_reservoir, - int num_timesteps, - double* giuh_ordinates, - int num_giuh_ordinates, - double *runoff_queue_m_per_timestep, - double field_capacity_atm_press_fraction, - double water_table_slope, - double Schaake_adjusted_magic_constant_by_soil_type, - double lateral_flow_linear_reservoir_constant, - double K_nash, - int num_lateral_flow_nash_reservoirs, - double* nash_storage, - int yes_aorc, - aorc_forcing_data* aorc_data, - double* rain_rate, - int& num_added_fluxes, - tshirt_c_result_fluxes* fluxes) -{ - // Do this to be safe ... - num_added_fluxes = 0; - - /* Commenting out, since this will need to return values rather than write results to a file. - FILE *in_fptr; - FILE *out_fptr; - FILE *out_debug_fptr; - */ - - // local variables - //------------------------------------------------------------------- - int i; - int tstep; - double upper_lim; - double lower_lim; - double diff; - - double catchment_area_km2 = 5.0; // in the range of our desired size - double drainage_density_km_per_km2 = 3.5; // this is approx. the average blue line drainage density for CONUS - - double timestep_h; - - // forcing - double timestep_rainfall_input_m; - //int yes_aorc = TRUE; // change to TRUE to read in entire AORC precip/met. forcing data set. - - //Schaake partitioning function parameters - - //double Schaake_adjusted_magic_constant_by_soil_type; - double Schaake_output_runoff_m; - double infiltration_depth_m; - - // GIUH state & parameters - // The GIUH ordinates (and ordinates count) now has to be passed in, rather than hard coded - //double *giuh_ordinates = NULL; // assumed GIUH hydrograph ordinates - //double *runoff_queue_m_per_timestep = NULL; - //int num_giuh_ordinates; - double giuh_runoff_m; - double soil_reservoir_storage_deficit_m; // the available space in the soil conceptual reservoir - - // groundwater storage parameters and state - - // lateral flow function parameters - double lateral_flow_threshold_storage_m; - double field_capacity_storage_threshold_m; - //double lateral_flow_linear_reservoir_constant; - //int num_lateral_flow_nash_reservoirs; - //double K_nash; // lateral_flow_nash_cascade_reservoir_constant; - //double *nash_storage = NULL; - double assumed_near_channel_water_table_slope; // [L/L] - double nash_lateral_runoff_m; - - // calculated flux variables - double flux_overland_m; // flux of surface runoff that goes through the GIUH convolution process - double flux_perc_m = 0.0; // flux from soil to deeper groundwater reservoir - double flux_lat_m; // lateral flux in the subsurface to the Nash cascade - double flux_from_deep_gw_to_chan_m; // flux from the deep reservoir into the channels - double gw_reservoir_storage_deficit_m; // the available space in the conceptual groundwater reservoir - double primary_flux, secondary_flux; // temporary vars. - - //double field_capacity_atm_press_fraction; // [-] - double soil_water_content_at_field_capacity; // [V/V] - - double atm_press_Pa = STANDARD_ATMOSPHERIC_PRESSURE_PASCALS; - double unit_weight_water_N_per_m3 = WATER_SPECIFIC_WEIGHT; - - double H_water_table_m; // Hwt in NWM/t-shirt parameter equiv. doc discussed between Eqn's 4 and 5. - // this is the distance down to a fictitious water table from the point there the - // soil water content is assumed equal to field capacity at the middle of lowest discretization - double trigger_z_m; // this is the distance up from the bottom of the soil that triggers lateral flow when - // the soil water content equals field capacity. 0.5 for center of bottom discretization - double Omega; // The limits of integration used in Eqn. 5 of the parameter equivalence document. - - double Qout_m; // the sum of giuh, nash-cascade, and base flow outputs m per timestep - - // mass balance variables. These all store cumulative discharges per unit catchment area [m3/m2] or [m] - //----------------------- - - // TODO: probably need to eventually move these out of here to realization/formulation, since now this only handles one time step at a time - - double volstart = 0.0; - double vol_sch_runoff = 0.0; - double vol_sch_infilt = 0.0; - - double vol_out_giuh = 0.0; - double vol_end_giuh = 0.0; - - double vol_to_gw = 0.0; - double vol_in_gw_start = 0.0; - double vol_in_gw_end = 0.0; - double vol_from_gw = 0.0; - - double vol_in_nash = 0.0; - double vol_in_nash_end = 0.0; // note the nash cascade is empty at start of simulation. - double vol_out_nash = 0.0; - - double vol_soil_start = 0.0; - double vol_to_soil = 0.0; - double vol_soil_to_lat_flow = 0.0; - double vol_soil_to_gw = 0.0; // this should equal vol_to_gw - double vol_soil_end = 0.0; - // note, vol_from_soil_to_gw is same as vol_to_gw. - - double volin = 0.0; - double volout = 0.0; - double volend = 0.0; - - /* Commenting out, since this will need to return values rather than write results to a file. - if ((out_fptr = fopen("test.out", "w")) == NULL) { - printf("Can't open output file\n"); - return -1; - } - */ - - #ifdef DEBUG - /* Commenting out, since this will need to return values rather than write results to a file. - if ((out_debug_fptr = fopen("debug.out", "w")) == NULL) { - printf("Can't open output file\n"); - return -1; - } - */ - #endif - - //initialize simulation constants - //--------------------------- - timestep_h = 1.0; - - trigger_z_m = 0.5; // distance from the bottom of the soil column to the center of the lowest discretization - - // equation 3 from NWM/t-shirt parameter equivalence document - H_water_table_m = field_capacity_atm_press_fraction * atm_press_Pa / unit_weight_water_N_per_m3; - soil_water_content_at_field_capacity=NWM_soil_params.smcmax * - pow(H_water_table_m / NWM_soil_params.satpsi, (1.0 / NWM_soil_params.bb)); - - // solve the integral given by Eqn. 5 in the parameter equivalence document. - // this equation calculates the amount of water stored in the 2 m thick soil column when the water content - // at the center of the bottom discretization (trigger_z_m) is at field capacity - Omega = H_water_table_m - trigger_z_m; - lower_lim = pow(Omega, (1.0 - (1.0 / NWM_soil_params.bb))) / (1.0 - (1.0 / NWM_soil_params.bb)); - upper_lim = pow(Omega + NWM_soil_params.D, (1.0 - (1.0 / NWM_soil_params.bb))) / (1.0 - (1.0 / NWM_soil_params.bb)); - field_capacity_storage_threshold_m = NWM_soil_params.smcmax - * pow((1.0 / NWM_soil_params.satpsi), (-1.0 / NWM_soil_params.bb)) - * (upper_lim - lower_lim); - - // initialize lateral flow function parameters - //--------------------------------------------- - // TODO: think this is slope from params (and tshirt params), but verify - //assumed_near_channel_water_table_slope=0.01; // [L/L] - assumed_near_channel_water_table_slope = water_table_slope; - lateral_flow_threshold_storage_m = field_capacity_storage_threshold_m; // making them the same, but they don't have 2B - - if (num_lateral_flow_nash_reservoirs > MAX_NUM_NASH_CASCADE) { - fprintf(stdout, "Number of Nash Cascade linear reservoirs greater than MAX_NUM_NASH_CASCADE.\n"); - return -2; - } - - volstart += gw_reservoir.storage_m; // initial mass balance checks in g.w. reservoir - vol_in_gw_start = gw_reservoir.storage_m; - - volstart += soil_reservoir.storage_m; // initial mass balance checks in soil reservoir - vol_soil_start = soil_reservoir.storage_m; - - // Rather than read in forcing or rainfall data from file (as in original) function receives directly as param - if(yes_aorc == TRUE) { - for (i = 0; i < num_timesteps; i++) { - // assumed 1000 kg/m3 density of water. The AORC data result is mm/h, so also convert to m/h. - // TODO: test to confirm the conversion to meters works out correctly - rain_rate[i] = ((double) aorc_data[i].precip_kg_per_m2) / 1000; - // Since we need to deal in m/s, and just got mm/s, make adjustment - } - } - - // run the model for num_timesteps - //--------------------------------------- - - double lateral_flux; // flux from soil to lateral flow Nash cascade +to cascade - double percolation_flux; // flux from soil to gw nonlinear researvoir, +downward - double oldval; - - //################################################################################### - //############################ TIME LOOP ################################# - //################################################################################### - - // For this, just do a single time step, though leaving coding structure in place for loop (for now) to minimize changes - //num_timesteps=num_rain_dat+279; // run a little bit beyond the rain data to see what happens. - for(tstep = 0; tstep < num_timesteps; tstep++) { - if (tstep < num_timesteps) - timestep_rainfall_input_m = rain_rate[tstep]; - else - timestep_rainfall_input_m = 0.0; - - volin += timestep_rainfall_input_m; - - //################################################## - // partition rainfall using Schaake function - //################################################## - - soil_reservoir_storage_deficit_m=(NWM_soil_params.smcmax * NWM_soil_params.D - soil_reservoir.storage_m); - - Schaake_partitioning_scheme(timestep_h, Schaake_adjusted_magic_constant_by_soil_type, - soil_reservoir_storage_deficit_m, timestep_rainfall_input_m, - &Schaake_output_runoff_m, &infiltration_depth_m); - - // check to make sure that there is storage available in soil to hold the water that does not runoff - //-------------------------------------------------------------------------------------------------- - if (soil_reservoir_storage_deficit_m < infiltration_depth_m) { - Schaake_output_runoff_m += (infiltration_depth_m - - soil_reservoir_storage_deficit_m); // put won't fit back into runoff - infiltration_depth_m = soil_reservoir_storage_deficit_m; - soil_reservoir.storage_m = soil_reservoir.storage_max_m; - } - - flux_overland_m = Schaake_output_runoff_m; - - vol_sch_runoff += flux_overland_m; - vol_sch_infilt += infiltration_depth_m; - - if (flux_perc_m > soil_reservoir_storage_deficit_m) { - diff = flux_perc_m - soil_reservoir_storage_deficit_m; // the amount that there is not capacity ffor - infiltration_depth_m = soil_reservoir_storage_deficit_m; - vol_sch_runoff += diff; // send excess water back to GIUH runoff - vol_sch_infilt -= diff; // correct overprediction of infilt. - flux_overland_m += diff; // bug found by Nels. This was missing and fixes it. - } - - vol_to_soil += infiltration_depth_m; - soil_reservoir.storage_m += infiltration_depth_m; // put the infiltrated water in the soil. - - // calculate fluxes from the soil storage into the deep groundwater (percolation) and to lateral subsurface flow - //-------------------------------------------------------------------------------------------------------------- - conceptual_reservoir_flux_calc(&soil_reservoir, &percolation_flux, &lateral_flux); - - flux_perc_m = percolation_flux; // m/h <<<<<<<<<<< flux of percolation from soil to g.w. reservoir >>>>>>>>> - - flux_lat_m = lateral_flux; // m/h <<<<<<<<<<< flux into the lateral flow Nash cascade >>>>>>>> - - // calculate flux of base flow from deep groundwater reservoir to channel - //-------------------------------------------------------------------------- - //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - gw_reservoir_storage_deficit_m = gw_reservoir.storage_max_m - gw_reservoir.storage_m; - - // limit amount transferred to deficit iff there is insufficient avail. storage - if (flux_perc_m > gw_reservoir_storage_deficit_m) { - diff = flux_perc_m - gw_reservoir_storage_deficit_m; - flux_perc_m = gw_reservoir_storage_deficit_m; - vol_sch_runoff += diff; // send excess water back to GIUH runoff - vol_sch_infilt -= diff; // correct overprediction of infilt. - } - - vol_to_gw += flux_perc_m; - vol_soil_to_gw += flux_perc_m; - - gw_reservoir.storage_m += flux_perc_m; - soil_reservoir.storage_m -= flux_perc_m; - soil_reservoir.storage_m -= flux_lat_m; - vol_soil_to_lat_flow += flux_lat_m; //TODO add this to nash cascade as input - volout = volout + flux_lat_m; - - conceptual_reservoir_flux_calc(&gw_reservoir, &primary_flux, &secondary_flux); - - flux_from_deep_gw_to_chan_m = primary_flux; // m/h <<<<<<<<<< BASE FLOW FLUX >>>>>>>>> - vol_from_gw += flux_from_deep_gw_to_chan_m; - - // in the instance of calling the gw reservoir the secondary flux should be zero- verify - if (is_fabs_less_than_epsilon(secondary_flux, 1.0e-09) == FALSE) printf("problem with nonzero flux point 1\n"); - - - // adjust state of deep groundwater conceptual nonlinear reservoir - //----------------------------------------------------------------- - - gw_reservoir.storage_m -= flux_from_deep_gw_to_chan_m; - - - // Solve the convolution integral ffor this time step - - giuh_runoff_m = convolution_integral(Schaake_output_runoff_m, num_giuh_ordinates, - giuh_ordinates, runoff_queue_m_per_timestep); - vol_out_giuh += giuh_runoff_m; - - volout += giuh_runoff_m; - volout += flux_from_deep_gw_to_chan_m; - - // Route lateral flow through the Nash cascade. - nash_lateral_runoff_m = nash_cascade(flux_lat_m, num_lateral_flow_nash_reservoirs, - K_nash, nash_storage); - vol_in_nash += flux_lat_m; - vol_out_nash += nash_lateral_runoff_m; - - #ifdef DEBUG - //fprintf(out_debug_fptr,"%d %lf %lf\n",tstep,flux_lat_m,nash_lateral_runoff_m); - #endif - - Qout_m = giuh_runoff_m + nash_lateral_runoff_m + flux_from_deep_gw_to_chan_m; - - // <<<<<<<>>>>>>>> - // These fluxs are all in units of meters per time step. Multiply them by the "catchment_area_km2" variable - // to convert them into cubic meters per time step. - - fluxes[tstep].timestep_rainfall_input_m = timestep_rainfall_input_m; - fluxes[tstep].Schaake_output_runoff_m = Schaake_output_runoff_m; - fluxes[tstep].giuh_runoff_m = giuh_runoff_m; - fluxes[tstep].nash_lateral_runoff_m = nash_lateral_runoff_m; - fluxes[tstep].flux_from_deep_gw_to_chan_m = flux_from_deep_gw_to_chan_m; - fluxes[tstep].Qout_m = Qout_m; - - // Increment reference variable tracking how many fluxes got appended to passed "fluxes" reference array. - num_added_fluxes++; - } - - // - // PERFORM MASS BALANCE CHECKS AND WRITE RESULTS TO stderr. - //---------------------------------------------------------- - - volend = soil_reservoir.storage_m + gw_reservoir.storage_m; - vol_in_gw_end = gw_reservoir.storage_m; - - #ifdef DEBUG - //fclose(out_debug_fptr); - #endif - - // the GIUH queue might have water in it at the end of the simulation, so sum it up. - for(i = 0; i < num_giuh_ordinates; i++) - vol_end_giuh += runoff_queue_m_per_timestep[i]; - - for(i = 0; i < num_lateral_flow_nash_reservoirs; i++) - vol_in_nash_end += nash_storage[i]; - - vol_soil_end = soil_reservoir.storage_m; - - // TODO: need ensure mass balance check is performed correctly and verifies - - return 0; -} diff --git a/src/realizations/catchment/CMakeLists.txt b/src/realizations/catchment/CMakeLists.txt index 59f9d35a23..54d28490b1 100644 --- a/src/realizations/catchment/CMakeLists.txt +++ b/src/realizations/catchment/CMakeLists.txt @@ -16,14 +16,12 @@ target_include_directories(realizations_catchment PUBLIC ${PROJECT_SOURCE_DIR}/models ${PROJECT_SOURCE_DIR}/models/kernels ${PROJECT_SOURCE_DIR}/models/kernels/evapotranspiration - ${PROJECT_SOURCE_DIR}/models/tshirt/include ) target_link_libraries(realizations_catchment PUBLIC ${CMAKE_DL_LIBS} NGen::core_catchment NGen::core_catchment_giuh - NGen::models_tshirt NGen::geojson NGen::kernels_evapotranspiration ) diff --git a/src/realizations/catchment/Tshirt_C_Realization.cpp b/src/realizations/catchment/Tshirt_C_Realization.cpp deleted file mode 100644 index e2a1ff2eb5..0000000000 --- a/src/realizations/catchment/Tshirt_C_Realization.cpp +++ /dev/null @@ -1,684 +0,0 @@ -#include "Tshirt_C_Realization.hpp" -#include "Constants.h" -#include -#include "tshirt_c.h" -#include "GIUH.hpp" -#include -#include -#include -#include - -using namespace realization; - -Tshirt_C_Realization::Tshirt_C_Realization(forcing_params forcing_config, - utils::StreamHandler output_stream, - double soil_storage, - double groundwater_storage, - bool storage_values_are_ratios, - std::string catchment_id, - giuh::GiuhJsonReader &giuh_json_reader, - tshirt::tshirt_params params, - const std::vector &nash_storage) - : Tshirt_C_Realization::Tshirt_C_Realization(std::move(forcing_config), output_stream, soil_storage, - groundwater_storage, storage_values_are_ratios, - std::move(catchment_id), - giuh_json_reader.extract_cumulative_frequency_ordinates(catchment_id), - params, nash_storage) -{ - -} - -Tshirt_C_Realization::Tshirt_C_Realization(forcing_params forcing_config, - utils::StreamHandler output_stream, - double soil_storage, - double groundwater_storage, - bool storage_values_are_ratios, - std::string catchment_id, - std::vector giuh_ordinates, - tshirt::tshirt_params params, - const std::vector &nash_storage) - //: Catchment_Formulation(catchment_id, std::move(std::make_unique(forcing_config)), output_stream), catchment_id(std::move(catchment_id)), - : Catchment_Formulation(catchment_id, std::move(std::make_unique(forcing_config)), output_stream), catchment_id(std::move(catchment_id)), - giuh_cdf_ordinates(std::move(giuh_ordinates)), params(std::make_shared(params)), nash_storage(nash_storage), c_soil_params(NWM_soil_parameters()), - groundwater_conceptual_reservoir(conceptual_reservoir()), soil_conceptual_reservoir(conceptual_reservoir()), - c_aorc_params(aorc_forcing_data()) -{ - // initialize 0 values if necessary in Nash Cascade storage vector - if (this->nash_storage.empty()) { - for (int i = 0; i < params.nash_n; i++) { - this->nash_storage.push_back(0.0); - } - } - - // Create this with 0 values initially - giuh_runoff_queue_per_timestep = std::vector(giuh_cdf_ordinates.size() + 1); - for (int i = 0; i < giuh_cdf_ordinates.size() + 1; i++) { - giuh_runoff_queue_per_timestep.push_back(0.0); - } - - fluxes = std::vector>(); - - // Convert params to struct for C-impl - sync_c_storage_params(); - - // TODO: Convert aorc to struct for C-impl - - this->init_ground_water_reservoir(groundwater_storage, storage_values_are_ratios); - this->init_soil_reservoir(soil_storage, storage_values_are_ratios); - -} - -Tshirt_C_Realization::Tshirt_C_Realization(forcing_params forcing_config, - utils::StreamHandler output_stream, - double soil_storage, - double groundwater_storage, - bool storage_values_are_ratios, - std::string catchment_id, - giuh::GiuhJsonReader &giuh_json_reader, - double maxsmc, - double wltsmc, - double satdk, - double satpsi, - double slope, - double b, - double multiplier, - double alpha_fc, - double Klf, - double Kn, - int nash_n, - double Cgw, - double expon, - double max_gw_storage, - const std::vector &nash_storage) - : Tshirt_C_Realization::Tshirt_C_Realization(std::move(forcing_config), output_stream, soil_storage, - groundwater_storage, storage_values_are_ratios, - std::move(catchment_id), giuh_json_reader, - tshirt::tshirt_params(maxsmc, wltsmc, satdk, satpsi, slope, b, - multiplier, alpha_fc, Klf, Kn, nash_n, - Cgw, expon, max_gw_storage), - nash_storage) -{ -} - -Tshirt_C_Realization::Tshirt_C_Realization( - std::string id, - forcing_params forcing_config, - utils::StreamHandler output_stream -//) : Catchment_Formulation(std::move(id), std::move(std::make_unique(forcing_config)), output_stream) { -) : Catchment_Formulation(std::move(id), std::move(std::make_unique(forcing_config)), output_stream) { - fluxes = std::vector>(); -} - -Tshirt_C_Realization::Tshirt_C_Realization( - std::string id, - shared_ptr forcing_provider, - utils::StreamHandler output_stream -) : Catchment_Formulation(std::move(id), forcing_provider, output_stream) { - fluxes = std::vector>(); -} - -Tshirt_C_Realization::~Tshirt_C_Realization() -{ - //destructor -} - -/** - * Return ``0``, as (for now) this type does not otherwise include ET within its calculations. - * - * @return ``0`` - */ -double Tshirt_C_Realization::calc_et() { - return 0; -} - -void Tshirt_C_Realization::create_formulation(geojson::PropertyMap properties) { - // TODO: don't particularly like the idea of "creating" the formulation, parameter constructs, etc., inside this - // type but outside the constructor. - // TODO: look at creating a factory or something, rather than an instance, for doing this type of thing. - - // TODO: (if this remains and doesn't get replaced with factory) protect against this being called after calls to - // get_response have started being made, or else the reservoir values are going to be jacked up. - this->validate_parameters(properties); - - catchment_id = this->get_id(); - - //dt = options.at("timestep").as_natural_number(); - - params = std::make_shared(tshirt_params{ - properties.at("maxsmc").as_real_number(), //maxsmc FWRFH - properties.at("wltsmc").as_real_number(), //wltsmc from fred_t-shirt.c FIXME NOT USED IN TSHIRT?!?! - properties.at("satdk").as_real_number(), //satdk FWRFH - properties.at("satpsi").as_real_number(), //satpsi FIXME what is this and what should its value be? - properties.at("slope").as_real_number(), //slope - properties.at("scaled_distribution_fn_shape_parameter").as_real_number(), //b bexp? FWRFH - properties.at("multiplier").as_real_number(), //multipier FIXMME (lksatfac) - properties.at("alpha_fc").as_real_number(), //aplha_fc field_capacity_atm_press_fraction - properties.at("Klf").as_real_number(), //Klf lateral flow nash coefficient? - properties.at("Kn").as_real_number(), //Kn Kn 0.001-0.03 F Nash Cascade coeeficient - static_cast(properties.at("nash_n").as_natural_number()), //number_lateral_flow_nash_reservoirs - properties.at("Cgw").as_real_number(), //fred_t-shirt gw res coeeficient (per h) - properties.at("expon").as_real_number(), //expon FWRFH - properties.at("max_groundwater_storage_meters").as_real_number() //max_gw_storage Sgwmax FWRFH - }); - // Very important this also gets done - sync_c_storage_params(); - - init_ground_water_reservoir(properties.at("groundwater_storage_percentage").as_real_number(), true); - init_soil_reservoir(properties.at("soil_storage_percentage").as_real_number(), true); - - nash_storage = properties.at("nash_storage").as_real_vector(); - - geojson::JSONProperty giuh = properties.at("giuh"); - - // Since this implementation really just cares about the ordinates, allow them to be passed directly here, or read - // from a separate file - if (giuh.has_key("cdf_ordinates")) { - giuh_cdf_ordinates = giuh.at("cdf_ordinates").as_real_vector(); - } - else { - std::vector missing_parameters; - if (!giuh.has_key("giuh_path")) { - missing_parameters.emplace_back("giuh_path"); - } - if (!giuh.has_key("crosswalk_path")) { - missing_parameters.emplace_back("crosswalk_path"); - } - if (!missing_parameters.empty()) { - std::string message = "A giuh configuration cannot be created for '" + catchment_id + "'; the following parameters are missing: "; - - for (int missing_parameter_index = 0; missing_parameter_index < missing_parameters.size(); missing_parameter_index++) { - message += missing_parameters[missing_parameter_index]; - - if (missing_parameter_index < missing_parameters.size() - 1) { - message += ", "; - } - } - - throw std::runtime_error(message); - } - - std::unique_ptr giuh_reader = std::make_unique( - giuh.at("giuh_path").as_string(), - giuh.at("crosswalk_path").as_string() - ); - - std::shared_ptr giuh_kernel = giuh_reader->get_giuh_kernel_for_id(catchment_id); - giuh_kernel->set_interpolation_regularity_seconds(3600); - // This needs to have all but the first interpolated incremental value (which is always 0 from the kernel), so - giuh_cdf_ordinates = std::vector(giuh_kernel->get_interpolated_incremental_runoff().size() - 1); - for (int i = 0; i < giuh_cdf_ordinates.size(); ++ i) { - giuh_cdf_ordinates[i] = giuh_kernel->get_interpolated_incremental_runoff()[i + 1]; - } - } - // Create this with 0 values initially - giuh_runoff_queue_per_timestep = std::vector(giuh_cdf_ordinates.size() + 1); - for (int i = 0; i < giuh_cdf_ordinates.size() + 1; i++) { - giuh_runoff_queue_per_timestep.push_back(0.0); - } -} - -void Tshirt_C_Realization::create_formulation(boost::property_tree::ptree &config, geojson::PropertyMap *global) { - // TODO: don't particularly like the idea of "creating" the formulation, parameter constructs, etc., inside this - // type but outside the constructor. - // TODO: look at creating a factory or something, rather than an instance, for doing this type of thing. - - // TODO: (if this remains and doesn't get replaced with factory) protect against this being called after calls to - // get_response have started being made, or else the reservoir values are going to be jacked up. - - geojson::PropertyMap options = this->interpret_parameters(config, global); - - catchment_id = this->get_id(); - - //dt = options.at("timestep").as_natural_number(); - - params = std::make_shared(tshirt_params{ - options.at("maxsmc").as_real_number(), //maxsmc FWRFH - options.at("wltsmc").as_real_number(), //wltsmc from fred_t-shirt.c FIXME NOT USED IN TSHIRT?!?! - options.at("satdk").as_real_number(), //satdk FWRFH - options.at("satpsi").as_real_number(), //satpsi FIXME what is this and what should its value be? - options.at("slope").as_real_number(), //slope - options.at("scaled_distribution_fn_shape_parameter").as_real_number(), //b bexp? FWRFH - options.at("multiplier").as_real_number(), //multipier FIXMME (lksatfac) - options.at("alpha_fc").as_real_number(), //aplha_fc field_capacity_atm_press_fraction - options.at("Klf").as_real_number(), //Klf lateral flow nash coefficient? - options.at("Kn").as_real_number(), //Kn Kn 0.001-0.03 F Nash Cascade coeeficient - static_cast(options.at("nash_n").as_natural_number()), //number_lateral_flow_nash_reservoirs - options.at("Cgw").as_real_number(), //fred_t-shirt gw res coeeficient (per h) - options.at("expon").as_real_number(), //expon FWRFH - options.at("max_groundwater_storage_meters").as_real_number() //max_gw_storage Sgwmax FWRFH - }); - // Very important this also gets done - sync_c_storage_params(); - - // TODO: might need to have this handle the ET params also - - init_ground_water_reservoir(options.at("groundwater_storage_percentage").as_real_number(), true); - init_soil_reservoir(options.at("soil_storage_percentage").as_real_number(), true); - - nash_storage = options.at("nash_storage").as_real_vector(); - - geojson::JSONProperty giuh = options.at("giuh"); - - // Since this implementation really just cares about the ordinates, allow them to be passed directly here, or read - // from a separate file - if (giuh.has_key("cdf_ordinates")) { - giuh_cdf_ordinates = giuh.at("cdf_ordinates").as_real_vector(); - } - else { - std::vector missing_parameters; - if (!giuh.has_key("giuh_path")) { - missing_parameters.emplace_back("giuh_path"); - } - if (!giuh.has_key("crosswalk_path")) { - missing_parameters.emplace_back("crosswalk_path"); - } - if (!missing_parameters.empty()) { - std::string message = "A giuh configuration cannot be created for '" + catchment_id + "'; the following parameters are missing: "; - - for (int missing_parameter_index = 0; missing_parameter_index < missing_parameters.size(); missing_parameter_index++) { - message += missing_parameters[missing_parameter_index]; - - if (missing_parameter_index < missing_parameters.size() - 1) { - message += ", "; - } - } - - throw std::runtime_error(message); - } - - std::unique_ptr giuh_reader = std::make_unique( - giuh.at("giuh_path").as_string(), - giuh.at("crosswalk_path").as_string() - ); - - - std::shared_ptr giuh_kernel = giuh_reader->get_giuh_kernel_for_id(catchment_id); - giuh_kernel->set_interpolation_regularity_seconds(3600); - // This needs to have all but the first interpolated incremental value (which is always 0 from the kernel), so - giuh_cdf_ordinates = std::vector(giuh_kernel->get_interpolated_incremental_runoff().size() - 1); - for (int i = 0; i < giuh_cdf_ordinates.size(); ++ i) { - giuh_cdf_ordinates[i] = giuh_kernel->get_interpolated_incremental_runoff()[i + 1]; - } - } - // Create this with 0 values initially - giuh_runoff_queue_per_timestep = std::vector(giuh_cdf_ordinates.size() + 1); - for (int i = 0; i < giuh_cdf_ordinates.size() + 1; i++) { - giuh_runoff_queue_per_timestep.push_back(0.0); - } - -} - -std::string Tshirt_C_Realization::get_formulation_type() { - return "tshirt_c"; -} - -double Tshirt_C_Realization::get_latest_flux_base_flow() { - return fluxes.empty() ? 0.0 : fluxes.back()->flux_from_deep_gw_to_chan_m; -} - -double Tshirt_C_Realization::get_latest_flux_giuh_runoff() { - return fluxes.empty() ? 0.0 : fluxes.back()->giuh_runoff_m; -} - -double Tshirt_C_Realization::get_latest_flux_lateral_flow() { - return fluxes.empty() ? 0.0 : fluxes.back()->nash_lateral_runoff_m; -} - -double Tshirt_C_Realization::get_latest_flux_surface_runoff() { - return fluxes.empty() ? 0.0 : fluxes.back()->Schaake_output_runoff_m; -} - -double Tshirt_C_Realization::get_latest_flux_total_discharge() { - return fluxes.empty() ? 0.0 : fluxes.back()->Qout_m; -} - -/** - * Get the number of output data variables made available from the calculations for enumerated time steps. - * - * @return The number of output data variables made available from the calculations for enumerated time steps. - */ -int Tshirt_C_Realization::get_output_item_count() { - return get_output_var_names().size(); -} - -/** - * Get a header line appropriate for a file made up of entries from this type's implementation of - * ``get_output_line_for_timestep``. - * - * Note that like the output generating function, this line does not include anything for time step. - * - * @return An appropriate header line for this type. - */ -std::string Tshirt_C_Realization::get_output_header_line(std::string delimiter) { - return boost::algorithm::join(get_output_header_fields(), delimiter); -} - -/** - * Get the values making up the header line from get_output_header_line(), but organized as a vector of strings. - * - * @return The values making up the header line from get_output_header_line() organized as a vector. - */ -const std::vector& Tshirt_C_Realization::get_output_header_fields() { - return OUTPUT_HEADER_FIELDS; -} - -/** - * Get a delimited string with all the output variable values for the given time step. - * - * This method is useful for preparing calculated data in a representation useful for output files, such as - * CSV files. - * - * The resulting string contains only the calculated output values for the time step, and not the time step - * index itself. - * - * An empty string is returned if the time step value is not in the range of valid time steps for which there - * are calculated values for all variables. - * - * The default delimiter is a comma. - * - * @param timestep The time step for which data is desired. - * @return A delimited string with all the output variable values for the given time step. - */ -std::string Tshirt_C_Realization::get_output_line_for_timestep(int timestep, std::string delimiter) { - // Check if the timestep is in bounds for the fluxes vector, and handle case when it isn't - if (timestep >= fluxes.size() || fluxes[timestep] == nullptr) { - return ""; - } - std::string output_str; - tshirt_c_result_fluxes flux_for_timestep = *fluxes[timestep]; - for (const std::string& name : get_output_var_names()) { - // Get a lambda that takes a fluxes struct and returns the right (double) member value from it from the name - std::function get_val_func = get_output_var_flux_extraction_func(name); - double output_var_value = get_val_func(flux_for_timestep); - output_str += output_str.empty() ? std::to_string(output_var_value) : "," + std::to_string(output_var_value); - } - return output_str; -} - -/** - * Get the names of the output data variables that are available from calculations for enumerated time steps. - * - * @return The names of the output data variables that are available from calculations for enumerated time steps. - */ -const std::vector &Tshirt_C_Realization::get_output_var_names() { - return OUTPUT_VARIABLE_NAMES; -} - -// TODO: don't care for this, as it could have the reference locations accidentally altered (also, raw pointer => bad) -//@robertbartel is this TODO resolved with these changes? -const std::vector& Tshirt_C_Realization::get_required_parameters() { - return REQUIRED_PARAMETERS; -} - -/** - * Execute the backing model formulation for the given time step, where it is of the specified size, and - * return the total discharge. - * - * Any inputs and additional parameters must be made available as instance members. - * - * Types should clearly document the details of their particular response output. - * - * @param t_index The index of the time step for which to run model calculations. - * @param d_delta_s The duration, in seconds, of the time step for which to run model calculations. - * @return The total discharge of the model for this time step. - */ -double Tshirt_C_Realization::get_response(time_step_t t_index, time_step_t t_delta_s) { - // TODO: check that t_delta_s is of approprate size - - // TODO: add some logic for ensuring the right precip data is gathered. This will need to include considerations - // for the cases when the dt is larger, smaller and equal to a specific, discrete data point available in the - // forcing. It may actually belong within the forcing object. - - // TODO: it also needs to account for getting the right precip data point (i.e., t_index may not be "next") - //Checking the time step used is consistent with that provided in forcing data - time_t t_delta = this->forcing->record_duration(); - if (t_delta != t_delta_s) { //Checking the time step used is consistent with that provided in forcing data - throw std::invalid_argument("Getting response using insonsistent time step with provided forcing data"); - } - //Negative t_index is not allowed - if (t_index < 0) { - throw std::invalid_argument("Getting response of negative time step in Tshirt C Realization is not allowed."); - } - time_t start_time = this->forcing->get_data_start_time(); - time_t stop_time = this->forcing->get_data_stop_time(); - time_t t_current = start_time + t_index * t_delta_s; - //Ensure model run does not exceed the end time of forcing - if (t_current > stop_time) { - throw std::invalid_argument("Getting response beyond time with available forcing."); - } - - double precip; - const std::string forcing_name = CSDMS_STD_NAME_LIQUID_EQ_PRECIP_RATE; - precip = this->forcing->get_value(CatchmentAggrDataSelector(this->catchment_id, CSDMS_STD_NAME_LIQUID_EQ_PRECIP_RATE, t_current, t_delta_s, ""), data_access::SUM); - int response_result = run_formulation_for_timestep(precip, t_delta_s); - // TODO: check t_index is the next expected time step to be calculated - - return fluxes.back()->Qout_m; -} - -std::function -Tshirt_C_Realization::get_output_var_flux_extraction_func(const std::string& var_name) { - // TODO: think about making this a lazily initialized member map - if (var_name == OUT_VAR_BASE_FLOW) { - return [](tshirt_c_result_fluxes flux) { return flux.flux_from_deep_gw_to_chan_m;}; - } - else if (var_name == OUT_VAR_GIUH_RUNOFF) { - return [](tshirt_c_result_fluxes flux) { return flux.giuh_runoff_m;}; - } - else if (var_name == OUT_VAR_LATERAL_FLOW) { - return [](tshirt_c_result_fluxes flux) { return flux.nash_lateral_runoff_m;}; - } - else if (var_name == OUT_VAR_RAINFALL) { - return [](tshirt_c_result_fluxes flux) { return flux.timestep_rainfall_input_m;}; - } - else if (var_name == OUT_VAR_SURFACE_RUNOFF) { - return [](tshirt_c_result_fluxes flux) { return flux.Schaake_output_runoff_m;}; - } - else if (var_name == OUT_VAR_TOTAL_DISCHARGE) { - return [](tshirt_c_result_fluxes flux) { return flux.Qout_m;}; - } - else { - throw std::invalid_argument("Cannot get values for unrecognized variable name " + var_name); - } -} - -/** - * Get a copy of the values for the given output variable at all available time steps. - * - * @param name - * @return A vector containing copies of the output value of the variable, indexed by time step. - */ -std::vector Tshirt_C_Realization::get_value(const std::string& name) { - // Generate a lambda that takes a fluxes struct and returns the right (double) member value from it from the name - std::function get_val_func = get_output_var_flux_extraction_func(name); - // Then, assuming we don't bail, use the lambda to build the result array - std::vector outputs = std::vector(fluxes.size()); - for (int i = 0; i < fluxes.size(); ++i) { - outputs[i] = get_val_func(*fluxes[i]); - } - return outputs; -} - -/** - * Run model formulation calculations for the next time step using the given input flux value in meters per second. - * - * The backing model works on a collection of time steps by receiving an associated collection of input fluxes. This - * is implemented by putting this input flux in a single-value vector and using it as the arg to a nested call to - * ``run_formulation_for_timesteps``, returning that result code. - * - * @param input_flux Input flux (typically expected to be just precipitation) in meters per second. - * @param t_delta_s The size of the time step in seconds - * @return The result code from the execution of the model time step calculations. - */ -int Tshirt_C_Realization::run_formulation_for_timestep(double input_flux, time_step_t t_delta_s) { - std::vector input_flux_in_vector{input_flux}; - return run_formulation_for_timesteps({input_flux}, {t_delta_s}); -} - -/** - * Run model formulation calculations for a series of time steps using the given collection of input flux values - * in meters per second. - * - * @param input_fluxes Ordered, per-time-step input flux (typically expected to be just precipitation) in meters - * per second. - * @param t_delta_s The sizes of each of the time steps in seconds - * @return The result code from the execution of the model time step calculations. - */ -int Tshirt_C_Realization::run_formulation_for_timesteps(std::vector input_fluxes, - std::vector t_deltas_s) { - int num_timesteps = (int) input_fluxes.size(); - - // FIXME: verify this needs to be independent like this - // FIXME: also make sure it shouldn't be parameterized somewhere, rather than hard-coded - double assumed_near_channel_water_table_slope = 0.01; - - // TODO: need some kind of guarantee that the vectors won't be resized and current buffer arrays won't be changed or - // removed (before the below "run" call finishes) - double* giuh_ordinates = &giuh_cdf_ordinates[0]; - double* giuh_runoff_queue = &giuh_runoff_queue_per_timestep[0]; - - //aorc_forcing_data empty_forcing[num_timesteps]; - aorc_forcing_data empty_forcing[1]; - - // Since the input_fluxes param values are in meters per second, this will need to do some conversions to what gets - // passed to Tshirt_C's run(), which expects meters per time step. - std::vector input_meters_per_time_step(input_fluxes.size()); - for (int i = 0; i < input_fluxes.size(); ++i) { - input_meters_per_time_step[i] = input_fluxes[i] * t_deltas_s[i]; - } - double* input_as_array = &input_meters_per_time_step[0]; - - tshirt_c_result_fluxes output_fluxes_as_array[num_timesteps]; - - // use this to sanity check the fluxes got added as expected to array - int num_added_fluxes = 0; - - int result = run(c_soil_params, - groundwater_conceptual_reservoir, - soil_conceptual_reservoir, - num_timesteps, - giuh_ordinates, - (int)giuh_cdf_ordinates.size(), - giuh_runoff_queue, - params->alpha_fc, - assumed_near_channel_water_table_slope, - params->Cschaake, - params->Klf, - params->Kn, - params->nash_n, - &nash_storage[0], - FALSE, - empty_forcing, - input_as_array, - num_added_fluxes, - output_fluxes_as_array); - - // Move fluxes over to member data structure - for (int i = 0; i < num_added_fluxes && i < num_timesteps; i++) { - fluxes.push_back(std::make_shared(output_fluxes_as_array[i])); - } - - return result; -} - -void Tshirt_C_Realization::init_ground_water_reservoir(double storage, bool storage_values_are_ratios) { - // Populate the groundwater conceptual reservoir data structure - //----------------------------------------------------------------------- - // one outlet, 0.0 threshold, nonlinear and exponential as in NWM - groundwater_conceptual_reservoir.is_exponential=TRUE; // set this true TRUE to use the exponential form of the discharge equation - groundwater_conceptual_reservoir.storage_max_m=params->max_groundwater_storage_meters; - - groundwater_conceptual_reservoir.coeff_primary=params->Cgw; // per h - groundwater_conceptual_reservoir.exponent_primary=params->expon; // linear iff 1.0, non-linear iff > 1.0 - groundwater_conceptual_reservoir.storage_threshold_primary_m=0.0; // 0.0 means no threshold applied - - groundwater_conceptual_reservoir.storage_threshold_secondary_m=0.0; // 0.0 means no threshold applied - groundwater_conceptual_reservoir.coeff_secondary=0.0; // 0.0 means that secondary outlet is not applied - groundwater_conceptual_reservoir.exponent_secondary=1.0; // linear - - groundwater_conceptual_reservoir.storage_m = init_reservoir_storage(storage_values_are_ratios, - storage, - groundwater_conceptual_reservoir.storage_max_m); -} - -double Tshirt_C_Realization::init_reservoir_storage(bool is_ratio, double amount, double max_amount) { - // Negative amounts are always ignored and just considered emtpy - if (amount < 0.0) { - return 0.0; - } - // When not a ratio (and positive), just return the literal amount - if (!is_ratio) { - return amount; - } - // When between 0 and 1, return the simple ratio computation - if (amount <= 1.0) { - return max_amount * amount; - } - // Otherwise, just return the literal amount, and assume the is_ratio value was invalid - // TODO: is this the best way to handle this? - else { - return amount; - } -} - -void Tshirt_C_Realization::init_soil_reservoir(double storage, bool storage_values_are_ratios) { - - double trigger_z_m = 0.5; // distance from bottom of soil column to the center of the lowest discretization - - // calculate the activation storage for the secondary lateral flow outlet in the soil nonlinear reservoir. - // following the method in the NWM/t-shirt parameter equivalence document, assuming field capacity soil - // suction pressure = 1/3 atm= field_capacity_atm_press_fraction * atm_press_Pa. - - // equation 3 from NWM/t-shirt parameter equivalence document - // This may need to be changed as follows later, but for now, use the constant value - //double H_water_table_m = params->alpha_fc * forcing.get_AORC_PRES_surface_Pa() / WATER_SPECIFIC_WEIGHT; - double H_water_table_m = params->alpha_fc * STANDARD_ATMOSPHERIC_PRESSURE_PASCALS / WATER_SPECIFIC_WEIGHT; - - - // solve the integral given by Eqn. 5 in the parameter equivalence document. - // this equation calculates the amount of water stored in the 2 m thick soil column when the water content - // at the center of the bottom discretization (trigger_z_m) is at field capacity - double Omega = H_water_table_m - trigger_z_m; - double lower_lim = pow(Omega, (1.0 - 1.0 / c_soil_params.bb)) / (1.0 - 1.0 / c_soil_params.bb); - double upper_lim = pow(Omega + c_soil_params.D, (1.0 - 1.0 / c_soil_params.bb)) / (1.0 - 1.0 / c_soil_params.bb); - - // initialize lateral flow function parameters - //--------------------------------------------- - double field_capacity_storage_threshold_m = - c_soil_params.smcmax * pow(1.0 / c_soil_params.satpsi, (-1.0 / c_soil_params.bb)) * - (upper_lim - lower_lim); - double lateral_flow_threshold_storage_m = field_capacity_storage_threshold_m; // making them the same, but they don't have 2B - - // Initialize the soil conceptual reservoir data structure. Indented here to highlight different purposes - //------------------------------------------------------------------------------------------------------------- - // soil conceptual reservoir first, two outlets, two thresholds, linear (exponent=1.0). - soil_conceptual_reservoir.is_exponential = FALSE; // set this true TRUE to use the exponential form of the discharge equation - // this should NEVER be set to true in the soil reservoir. - soil_conceptual_reservoir.storage_max_m = c_soil_params.smcmax * c_soil_params.D; - // vertical percolation parameters------------------------------------------------ - // TODO: should this get parameterized somehow? - soil_conceptual_reservoir.coeff_primary = c_soil_params.satdk * c_soil_params.slop * 3600.0; // m per h - soil_conceptual_reservoir.exponent_primary = 1.0; // 1.0=linear - soil_conceptual_reservoir.storage_threshold_primary_m = field_capacity_storage_threshold_m; - // lateral flow parameters -------------------------------------------------------- - soil_conceptual_reservoir.coeff_secondary = params->Klf; // 0.0 to deactiv. else =lateral_flow_linear_reservoir_constant; // m per h - soil_conceptual_reservoir.exponent_secondary = 1.0; // 1.0=linear - soil_conceptual_reservoir.storage_threshold_secondary_m = lateral_flow_threshold_storage_m; - - soil_conceptual_reservoir.storage_m = init_reservoir_storage(storage_values_are_ratios, - storage, - soil_conceptual_reservoir.storage_max_m); - -} - -void Tshirt_C_Realization::sync_c_storage_params() { - // Convert params to struct for C-impl - c_soil_params.D = params->depth; - c_soil_params.bb = params->b; - c_soil_params.mult = params->multiplier; - c_soil_params.satdk = params->satdk; - c_soil_params.satpsi = params->satpsi; - c_soil_params.slop = params->slope; - c_soil_params.smcmax = params->maxsmc; - c_soil_params.wltsmc = params->wltsmc; -} diff --git a/src/realizations/catchment/Tshirt_Realization.cpp b/src/realizations/catchment/Tshirt_Realization.cpp deleted file mode 100644 index f12d058858..0000000000 --- a/src/realizations/catchment/Tshirt_Realization.cpp +++ /dev/null @@ -1,337 +0,0 @@ -#include "giuh_kernel.hpp" -#include "Tshirt_Realization.hpp" -#include "TshirtErrorCodes.h" -#include "Catchment_Formulation.hpp" -using namespace realization; - -/* -Tshirt_Realization::Tshirt_Realization( - forcing_params forcing_config, - utils::StreamHandler output_stream, - double soil_storage_meters, - double groundwater_storage_meters, - std::string catchment_id, - giuh::GiuhJsonReader &giuh_json_reader, - tshirt::tshirt_params params, - const vector &nash_storage, - time_step_t t) - : Catchment_Formulation(catchment_id, forcing_config, output_stream), catchment_id(catchment_id), dt(t) -{ - this->params = ¶ms; - giuh_kernel = giuh_json_reader.get_giuh_kernel_for_id(this->catchment_id); - - // If the look-up failed in the reader for some reason, and we got back a null pointer ... - if (this->giuh_kernel == nullptr) { - // ... revert to a pass-through kernel - this->giuh_kernel = std::make_shared( - giuh::giuh_kernel(this->catchment_id, giuh_json_reader.get_associated_comid(this->catchment_id))); - } - - //FIXME not really used, don't call??? - //add_time(t, params.nash_n); - state[0] = std::make_shared(tshirt::tshirt_state(soil_storage_meters, groundwater_storage_meters, nash_storage)); - //state[0]->soil_storage_meters = soil_storage_meters; - //state[0]->groundwater_storage_meters = groundwater_storage_meters; - - for (int i = 0; i < params.nash_n; ++i) { - - state[0]->nash_cascade_storeage_meters[i] = nash_storage[i]; - } - - model = make_unique(tshirt::tshirt_model(params, state[0])); -} - -Tshirt_Realization::Tshirt_Realization( - forcing_params forcing_config, - utils::StreamHandler output_stream, - double soil_storage_meters, - double groundwater_storage_meters, - std::string catchment_id, - giuh::GiuhJsonReader &giuh_json_reader, - double maxsmc, - double wltsmc, - double satdk, - double satpsi, - double slope, - double b, - double multiplier, - double alpha_fc, - double Klf, - double Kn, - int nash_n, - double Cgw, - double expon, - double max_gw_storage, - const std::vector &nash_storage, - time_step_t t -) : Tshirt_Realization::Tshirt_Realization(forcing_config, output_stream, soil_storage_meters, groundwater_storage_meters, - catchment_id, giuh_json_reader, - tshirt::tshirt_params(maxsmc, wltsmc, satdk, satpsi, slope, b, multiplier, - alpha_fc, Klf, Kn, nash_n, Cgw, expon, max_gw_storage), - nash_storage, t) { - -} -*/ - -double Tshirt_Realization::calc_et() { - return 0; -} - -/** - * Execute the backing model formulation for the given time step, where it is of the specified size, and - * return the total discharge. - * - * Function reads input precipitation from ``forcing`` member variable. It also makes use of the params struct - * for ET params accessible via ``get_et_params``. - * - * @param t_index The index of the time step for which to run model calculations. - * @param d_delta_s The duration, in seconds, of the time step for which to run model calculations. - * @return The total discharge for this time step. - */ -double Tshirt_Realization::get_response(time_step_t t_index, time_step_t t_delta_s) { - //FIXME doesn't do anything, don't call??? - //add_time(t+1, params.nash_n); - // TODO: this is problematic, because what happens if the wrong t_index is passed? - time_t t_delta = this->forcing->record_duration(); - if (t_delta != t_delta_s) { //Checking the time step used is consistent with that provided in forcing data - throw std::invalid_argument("Getting response using insonsistent time step with provided forcing data"); - } - //t_index can't be negative - if (t_index < 0) { - throw std::invalid_argument("Getting response of negative time step in Tshirt C Realization is not allowed."); - } - time_t start_time = this->forcing->get_data_start_time(); - time_t stop_time = this->forcing->get_data_stop_time(); - time_t t_current = start_time + t_index * t_delta_s; - //Ensure model run does not exceed the end time of forcing - if (t_current > stop_time) { - throw std::invalid_argument("Getting response beyond time with available forcing."); - } - - double precip; - const std::string forcing_name = CSDMS_STD_NAME_LIQUID_EQ_PRECIP_RATE; - precip = this->forcing->get_value(CatchmentAggrDataSelector(this->catchment_id, CSDMS_STD_NAME_LIQUID_EQ_PRECIP_RATE, t_current, t_delta_s, ""), data_access::SUM); - //FIXME should this run "daily" or hourly (t) which should really be dt - //Do we keep an "internal dt" i.e. this->dt and reconcile with t? - int error = model->run(t_index, precip * t_delta_s / 1000, get_et_params_ptr()); - if(error == tshirt::TSHIRT_MASS_BALANCE_ERROR){ - std::cout<<"WARNING Tshirt_Realization::model mass balance error"<get_current_state(); - fluxes[t_index] = model->get_fluxes(); - double giuh = giuh_kernel->calc_giuh_output(t_index, fluxes[t_index]->surface_runoff_meters_per_second); - return fluxes[t_index]->soil_lateral_flow_meters_per_second + fluxes[t_index]->groundwater_flow_meters_per_second + - giuh; -} - -/** - * Get a formatted line of output values for the given time step as a delimited string. - * - * For this type, the output consists of only the total discharge amount per time step; i.e., the same value that was - * returned by ``get_response``. - * - * This method is useful for preparing calculated data in a representation useful for output files, such as - * CSV files. - * - * The resulting string will contain calculated values for applicable output variables for the particular - * formulation, as determined for the given time step. However, the string will not contain any - * representation of the time step itself. - * - * An empty string is returned if the time step value is not in the range of valid time steps for which there - * are calculated values for all variables. - * - * The default delimiter is a comma. - * - * @param timestep The time step for which data is desired. - * @return A delimited string with all the output variable values for the given time step. - */ -std::string Tshirt_Realization::get_output_line_for_timestep(int timestep, std::string delimiter) { - if (timestep >= fluxes.size()) { - return ""; - } - double discharge = fluxes[timestep]->soil_lateral_flow_meters_per_second + - fluxes[timestep]->groundwater_flow_meters_per_second + - giuh_kernel->calc_giuh_output(timestep, fluxes[timestep]->surface_runoff_meters_per_second); - return std::to_string(discharge); -} - -void Tshirt_Realization::create_formulation(geojson::PropertyMap properties) { - this->validate_parameters(properties); - - this->catchment_id = this->get_id(); - this->dt = properties.at("timestep").as_natural_number(); - - tshirt::tshirt_params tshirt_params{ - properties.at("maxsmc").as_real_number(), //maxsmc FWRFH - properties.at("wltsmc").as_real_number(), //wltsmc from fred_t-shirt.c FIXME NOT USED IN TSHIRT?!?! - properties.at("satdk").as_real_number(), //satdk FWRFH - properties.at("satpsi").as_real_number(), //satpsi FIXME what is this and what should its value be? - properties.at("slope").as_real_number(), //slope - properties.at("scaled_distribution_fn_shape_parameter").as_real_number(), //b bexp? FWRFH - properties.at("multiplier").as_real_number(), //multipier FIXMME (lksatfac) - properties.at("alpha_fc").as_real_number(), //aplha_fc field_capacity_atm_press_fraction - properties.at("Klf").as_real_number(), //Klf lateral flow nash coefficient? - properties.at("Kn").as_real_number(), //Kn Kn 0.001-0.03 F Nash Cascade coeeficient - static_cast(properties.at("nash_n").as_natural_number()), //number_lateral_flow_nash_reservoirs - properties.at("Cgw").as_real_number(), //fred_t-shirt gw res coeeficient (per h) - properties.at("expon").as_real_number(), //expon FWRFH - properties.at("max_groundwater_storage_meters").as_real_number() //max_gw_storage Sgwmax FWRFH - }; - - this->params = &tshirt_params; - - double soil_storage_meters = tshirt_params.max_soil_storage_meters * properties.at("soil_storage_percentage").as_real_number(); - double ground_water_storage = tshirt_params.max_groundwater_storage_meters * properties.at("groundwater_storage_percentage").as_real_number(); - - std::vector nash_storage = properties.at("nash_storage").as_real_vector(); - - this->state[0] = std::make_shared(tshirt::tshirt_state(soil_storage_meters, ground_water_storage, nash_storage)); - - for (int i = 0; i < tshirt_params.nash_n; ++i) { - - this->state[0]->nash_cascade_storeage_meters[i] = nash_storage[i]; - } - - this->model = make_unique(tshirt::tshirt_model(tshirt_params, this->state[0])); - - geojson::JSONProperty giuh = properties.at("giuh"); - - std::vector missing_parameters; - - if (!giuh.has_key("giuh_path")) { - missing_parameters.push_back("giuh_path"); - } - - if (!giuh.has_key("crosswalk_path")) { - missing_parameters.push_back("crosswalk_path"); - } - - if (missing_parameters.size() > 0) { - std::string message = "A giuh configuration cannot be created for '" + this->get_id() + "'; the following parameters are missing: "; - - for (int missing_parameter_index = 0; missing_parameter_index < missing_parameters.size(); missing_parameter_index++) { - message += missing_parameters[missing_parameter_index]; - - if (missing_parameter_index < missing_parameters.size() - 1) { - message += ", "; - } - } - - throw std::runtime_error(message); - } - - std::unique_ptr giuh_reader = std::make_unique( - giuh.at("giuh_path").as_string(), - giuh.at("crosswalk_path").as_string() - ); - - this->giuh_kernel = giuh_reader->get_giuh_kernel_for_id(this->catchment_id); - - if (this->giuh_kernel == nullptr) { - // ... revert to a pass-through kernel - this->giuh_kernel = std::make_shared( - giuh::giuh_kernel( - this->catchment_id, - giuh_reader->get_associated_comid(this->catchment_id) - ) - ); - } -} - -void Tshirt_Realization::create_formulation(boost::property_tree::ptree &config, geojson::PropertyMap *global) { - geojson::PropertyMap options = this->interpret_parameters(config, global); - - this->catchment_id = this->get_id(); - this->dt = options.at("timestep").as_natural_number(); - - tshirt::tshirt_params tshirt_params{ - options.at("maxsmc").as_real_number(), //maxsmc FWRFH - options.at("wltsmc").as_real_number(), //wltsmc from fred_t-shirt.c FIXME NOT USED IN TSHIRT?!?! - options.at("satdk").as_real_number(), //satdk FWRFH - options.at("satpsi").as_real_number(), //satpsi FIXME what is this and what should its value be? - options.at("slope").as_real_number(), //slope - options.at("scaled_distribution_fn_shape_parameter").as_real_number(), //b bexp? FWRFH - options.at("multiplier").as_real_number(), //multipier FIXMME (lksatfac) - options.at("alpha_fc").as_real_number(), //aplha_fc field_capacity_atm_press_fraction - options.at("Klf").as_real_number(), //Klf lateral flow nash coefficient? - options.at("Kn").as_real_number(), //Kn Kn 0.001-0.03 F Nash Cascade coeeficient - static_cast(options.at("nash_n").as_natural_number()), //number_lateral_flow_nash_reservoirs - options.at("Cgw").as_real_number(), //fred_t-shirt gw res coeeficient (per h) - options.at("expon").as_real_number(), //expon FWRFH - options.at("max_groundwater_storage_meters").as_real_number() //max_gw_storage Sgwmax FWRFH - }; - - this->params = &tshirt_params; - - double soil_storage_meters = tshirt_params.max_soil_storage_meters * options.at("soil_storage_percentage").as_real_number(); - double ground_water_storage = tshirt_params.max_groundwater_storage_meters * options.at("groundwater_storage_percentage").as_real_number(); - - std::vector nash_storage = options.at("nash_storage").as_real_vector(); - - this->state[0] = std::make_shared(tshirt::tshirt_state(soil_storage_meters, ground_water_storage, nash_storage)); - - for (int i = 0; i < tshirt_params.nash_n; ++i) { - - this->state[0]->nash_cascade_storeage_meters[i] = nash_storage[i]; - } - - this->model = make_unique(tshirt::tshirt_model(tshirt_params, this->state[0])); - - geojson::JSONProperty giuh = options.at("giuh"); - - std::vector missing_parameters; - - if (!giuh.has_key("giuh_path")) { - missing_parameters.push_back("giuh_path"); - } - - if (!giuh.has_key("crosswalk_path")) { - missing_parameters.push_back("crosswalk_path"); - } - - if (missing_parameters.size() > 0) { - std::string message = "A giuh configuration cannot be created for '" + this->get_id() + "'; the following parameters are missing: "; - - for (int missing_parameter_index = 0; missing_parameter_index < missing_parameters.size(); missing_parameter_index++) { - message += missing_parameters[missing_parameter_index]; - - if (missing_parameter_index < missing_parameters.size() - 1) { - message += ", "; - } - } - - throw std::runtime_error(message); - } - - std::unique_ptr giuh_reader = std::make_unique( - giuh.at("giuh_path").as_string(), - giuh.at("crosswalk_path").as_string() - ); - - this->giuh_kernel = giuh_reader->get_giuh_kernel_for_id(this->catchment_id); - - if (this->giuh_kernel == nullptr) { - // ... revert to a pass-through kernel - this->giuh_kernel = std::make_shared( - giuh::giuh_kernel( - this->catchment_id, - giuh_reader->get_associated_comid(this->catchment_id) - ) - ); - } -} - -void Tshirt_Realization::set_giuh_kernel(std::shared_ptr reader) { - this->giuh_kernel = reader->get_giuh_kernel_for_id(this->catchment_id); - - // If the look-up failed in the reader for some reason, and we got back a null pointer ... - if (this->giuh_kernel == nullptr) { - // ... revert to a pass-through kernel - this->giuh_kernel = std::make_shared( - giuh::giuh_kernel( - this->catchment_id, - reader->get_associated_comid(this->catchment_id) - ) - ); - } -} diff --git a/src/routing/CMakeLists.txt b/src/routing/CMakeLists.txt index 843615be56..252cef610c 100644 --- a/src/routing/CMakeLists.txt +++ b/src/routing/CMakeLists.txt @@ -18,7 +18,6 @@ target_link_libraries(routing PUBLIC pybind11::embed) # ${CMAKE_DL_LIBS} # NGen::core_catchment # NGen::core_catchment_giuh -# NGen::models_tshirt # NGen::geojson # NGen::kernels_evapotranspiration # ) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 3deaee8629..3fceacd9d2 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -67,14 +67,6 @@ endfunction() ########################## Individual tests ########################## -#add_automated_test(test_hymod models/hymod/include/HymodTest.cpp) -# -#add_automated_test( -# test_tshirt -# models/tshirt/include/TshirtTest.cpp -# realizations/catchments/Tshirt_C_Realization_Test.cpp -#) - ########################## GeoJSON Unit Tests add_test(test_geojson 4 @@ -318,11 +310,7 @@ add_test( ########################## Primary Combined Unit Test Target add_test( test_unit - 18 - models/hymod/include/Reservoir_Test.cpp - models/hymod/include/Reservoir_Timeless_Test.cpp - models/tshirt/include/TshirtTest.cpp - realizations/catchments/Tshirt_C_Realization_Test.cpp + 14 geojson/JSONProperty_Test.cpp geojson/JSONGeometry_Test.cpp geojson/Feature_Test.cpp @@ -351,11 +339,7 @@ add_test( # All automated tests add_test( test_all - 17 - models/hymod/include/Reservoir_Test.cpp - models/hymod/include/Reservoir_Timeless_Test.cpp - models/tshirt/include/TshirtTest.cpp - realizations/catchments/Tshirt_C_Realization_Test.cpp + 13 geojson/JSONProperty_Test.cpp geojson/JSONGeometry_Test.cpp geojson/Feature_Test.cpp diff --git a/test/data/model/tshirt/expected_results.csv b/test/data/model/tshirt/expected_results.csv deleted file mode 100644 index 7a19faceea..0000000000 --- a/test/data/model/tshirt/expected_results.csv +++ /dev/null @@ -1,999 +0,0 @@ -2457357.50000000,10.000000,1.522619,0.091357,0.000000,191.106140,191.197497 -2457357.54166667,0.000000,0.000000,0.776536,0.000000,177.198214,177.974750 -2457357.58333333,0.000000,0.000000,0.426333,0.093168,165.163302,165.682804 -2457357.62500000,0.000000,0.000000,0.182714,0.180746,154.643488,155.006948 -2457357.66666667,0.000000,0.000000,0.045679,0.262986,145.367169,145.675833 -2457357.70833333,0.000000,0.000000,0.000000,0.340128,137.124396,137.464524 -2457357.75000000,0.000000,0.000000,0.000000,0.412405,129.750237,130.162643 -2457357.79166667,0.000000,0.000000,0.000000,0.480040,123.113277,123.593317 -2457357.83333333,0.000000,0.000000,0.000000,0.543245,117.107476,117.650721 -2457357.87500000,0.000000,0.000000,0.000000,0.602226,111.646304,112.248530 -2457357.91666667,0.000000,0.000000,0.000000,0.657179,106.658437,107.315616 -2457357.95833333,0.000000,0.000000,0.000000,0.708293,102.084540,102.792833 -2457358.00000000,0.000000,0.000000,0.000000,0.755749,97.874845,98.630593 -2457358.04166667,0.000000,0.000000,0.000000,0.799719,93.987286,94.787005 -2457358.08333333,0.000000,0.000000,0.000000,0.840372,90.386066,91.226438 -2457358.12500000,0.000000,0.000000,0.000000,0.877865,87.040521,87.918387 -2457358.16666667,0.000000,0.000000,0.000000,0.912353,83.924234,84.836587 -2457358.20833333,0.000000,0.000000,0.000000,0.943981,81.014320,81.958301 -2457358.25000000,0.000000,0.000000,0.000000,0.972891,78.290850,79.263741 -2457358.29166667,0.000000,0.000000,0.000000,0.999216,75.736395,76.735611 -2457358.33333333,0.000000,0.000000,0.000000,1.023086,73.335636,74.358722 -2457358.37500000,0.000000,0.000000,0.000000,1.044625,71.075061,72.119685 -2457358.41666667,0.000000,0.000000,0.000000,1.063950,68.942699,70.006649 -2457358.45833333,0.000000,0.000000,0.000000,1.081176,66.927909,68.009085 -2457358.50000000,0.000000,0.000000,0.000000,1.096411,65.021203,66.117613 -2457358.54166667,0.000000,0.000000,0.000000,1.109758,63.214085,64.323844 -2457358.58333333,0.000000,0.000000,0.000000,1.121318,61.498934,62.620253 -2457358.62500000,0.000000,0.000000,0.000000,1.131186,59.868887,61.000073 -2457358.66666667,0.000000,0.000000,0.000000,1.139452,58.317750,59.457202 -2457358.70833333,0.000000,0.000000,0.000000,1.146205,56.839916,57.986121 -2457358.75000000,0.000000,0.000000,0.000000,1.151526,55.430302,56.581828 -2457358.79166667,0.000000,0.000000,0.000000,1.155497,54.084282,55.239779 -2457358.83333333,0.000000,0.000000,0.000000,1.158193,52.797641,53.955835 -2457358.87500000,0.200000,0.000536,0.000032,1.159688,51.566675,52.726395 -2457358.91666667,0.000000,0.000000,0.000273,1.160050,50.387566,51.547890 -2457358.95833333,0.000000,0.000000,0.000150,1.159521,49.257235,50.416907 -2457359.00000000,0.000000,0.000000,0.000064,1.157980,48.172715,49.330759 -2457359.04166667,0.000000,0.000000,0.000016,1.155487,47.131273,48.286776 -2457359.08333333,0.000000,0.000000,0.000000,1.152102,46.130393,47.282495 -2457359.12500000,0.500000,0.003337,0.000200,1.147880,45.168069,46.316150 -2457359.16666667,0.800000,0.008510,0.002213,1.142874,44.242002,45.387089 -2457359.20833333,0.200000,0.000536,0.005307,1.137569,43.349635,44.492510 -2457359.25000000,0.000000,0.000000,0.003057,1.132244,42.489388,43.624688 -2457359.29166667,0.000000,0.000000,0.001271,1.126388,41.659675,42.787334 -2457359.33333333,0.000000,0.000000,0.000320,1.119865,40.858899,41.979084 -2457359.37500000,0.000000,0.000000,0.000016,1.112719,40.085574,41.198309 -2457359.41666667,0.000000,0.000000,0.000000,1.104994,39.338314,40.443308 -2457359.45833333,0.000000,0.000000,0.000000,1.096731,38.615824,39.712555 -2457359.50000000,0.000000,0.000000,0.000000,1.087970,37.916894,39.004864 -2457359.54166667,0.000000,0.000000,0.000000,1.078747,37.240393,38.319139 -2457359.58333333,0.000000,0.000000,0.000000,1.069098,36.585259,37.654357 -2457359.62500000,0.000000,0.000000,0.000000,1.059057,35.950498,37.009555 -2457359.66666667,0.000000,0.000000,0.000000,1.048657,35.335176,36.383833 -2457359.70833333,0.000000,0.000000,0.000000,1.037927,34.738416,35.776344 -2457359.75000000,0.000000,0.000000,0.000000,1.026898,34.159393,35.186291 -2457359.79166667,0.000000,0.000000,0.000000,1.015596,33.597329,34.612925 -2457359.83333333,0.000000,0.000000,0.000000,1.004048,33.051492,34.055540 -2457359.87500000,0.000000,0.000000,0.000000,0.992279,32.521191,33.513470 -2457359.91666667,0.000000,0.000000,0.000000,0.980312,32.005776,32.986088 -2457359.95833333,0.000000,0.000000,0.000000,0.968170,31.504629,32.472799 -2457360.00000000,0.000000,0.000000,0.000000,0.955875,31.017168,31.973043 -2457360.04166667,0.000000,0.000000,0.000000,0.943446,30.542844,31.486290 -2457360.08333333,0.000000,0.000000,0.000000,0.930902,30.081133,31.012035 -2457360.12500000,0.000000,0.000000,0.000000,0.918262,29.631541,30.549803 -2457360.16666667,0.000000,0.000000,0.000000,0.905543,29.193600,30.099142 -2457360.20833333,0.000000,0.000000,0.000000,0.892760,28.766864,29.659624 -2457360.25000000,0.000000,0.000000,0.000000,0.879929,28.350911,29.230840 -2457360.29166667,0.000000,0.000000,0.000000,0.867064,27.945340,28.812404 -2457360.33333333,0.000000,0.000000,0.000000,0.854180,27.549768,28.403948 -2457360.37500000,0.000000,0.000000,0.000000,0.841288,27.163832,28.005120 -2457360.41666667,0.000000,0.000000,0.000000,0.828401,26.787187,27.615588 -2457360.45833333,0.000000,0.000000,0.000000,0.815530,26.419502,27.235032 -2457360.50000000,0.000000,0.000000,0.000000,0.802686,26.060464,26.863150 -2457360.54166667,0.000000,0.000000,0.000000,0.789878,25.709773,26.499651 -2457360.58333333,0.000000,0.000000,0.000000,0.777117,25.367144,26.144260 -2457360.62500000,0.000000,0.000000,0.000000,0.764410,25.032303,25.796712 -2457360.66666667,0.000000,0.000000,0.000000,0.751766,24.704989,25.456755 -2457360.70833333,0.000000,0.000000,0.000000,0.739193,24.384954,25.124147 -2457360.75000000,0.000000,0.000000,0.000000,0.726698,24.071959,24.798657 -2457360.79166667,0.000000,0.000000,0.000000,0.714287,23.765776,24.480063 -2457360.83333333,0.000000,0.000000,0.000000,0.701966,23.466187,24.168153 -2457360.87500000,0.000000,0.000000,0.000000,0.689743,23.172982,23.862725 -2457360.91666667,0.000000,0.000000,0.000000,0.677620,22.885962,23.563583 -2457360.95833333,0.000000,0.000000,0.000000,0.665605,22.604935,23.270540 -2457361.00000000,0.000000,0.000000,0.000000,0.653700,22.329716,22.983416 -2457361.04166667,0.000000,0.000000,0.000000,0.641911,22.060129,22.702039 -2457361.08333333,0.000000,0.000000,0.000000,0.630240,21.796004,22.426244 -2457361.12500000,0.000000,0.000000,0.000000,0.618693,21.537178,22.155871 -2457361.16666667,0.000000,0.000000,0.000000,0.607270,21.283496,21.890766 -2457361.20833333,0.000000,0.000000,0.000000,0.595977,21.034806,21.630783 -2457361.25000000,0.000000,0.000000,0.000000,0.584814,20.790965,21.375779 -2457361.29166667,0.000000,0.000000,0.000000,0.573785,20.551833,21.125618 -2457361.33333333,0.000000,0.000000,0.000000,0.562891,20.317277,20.880168 -2457361.37500000,0.000000,0.000000,0.000000,0.552135,20.087168,20.639303 -2457361.41666667,0.000000,0.000000,0.000000,0.541517,19.861382,20.402899 -2457361.45833333,0.000000,0.000000,0.000000,0.531039,19.639800,20.170840 -2457361.50000000,0.000000,0.000000,0.000000,0.520703,19.422307,19.943010 -2457361.54166667,0.000000,0.000000,0.000000,0.510509,19.208793,19.719301 -2457361.58333333,0.000000,0.000000,0.000000,0.500458,18.999149,19.499607 -2457361.62500000,0.000000,0.000000,0.000000,0.490550,18.793273,19.283824 -2457361.66666667,0.000000,0.000000,0.000000,0.480787,18.591067,19.071854 -2457361.70833333,0.000000,0.000000,0.000000,0.471168,18.392433,18.863601 -2457361.75000000,0.000000,0.000000,0.000000,0.461693,18.197280,18.658973 -2457361.79166667,0.000000,0.000000,0.000000,0.452363,18.005517,18.457880 -2457361.83333333,0.000000,0.000000,0.000000,0.443177,17.817059,18.260236 -2457361.87500000,0.000000,0.000000,0.000000,0.434135,17.631822,18.065956 -2457361.91666667,0.000000,0.000000,0.000000,0.425237,17.449725,17.874961 -2457361.95833333,0.000000,0.000000,0.000000,0.416481,17.270690,17.687171 -2457362.00000000,0.000000,0.000000,0.000000,0.407869,17.094641,17.502510 -2457362.04166667,0.000000,0.000000,0.000000,0.399398,16.921507,17.320905 -2457362.08333333,0.000000,0.000000,0.000000,0.391069,16.751216,17.142285 -2457362.12500000,0.000000,0.000000,0.000000,0.382880,16.583699,16.966579 -2457362.16666667,0.000000,0.000000,0.000000,0.374830,16.418891,16.793721 -2457362.20833333,0.000000,0.000000,0.000000,0.366919,16.256728,16.623646 -2457362.25000000,0.000000,0.000000,0.000000,0.359145,16.097146,16.456291 -2457362.29166667,0.000000,0.000000,0.000000,0.351507,15.940087,16.291594 -2457362.33333333,0.000000,0.000000,0.000000,0.344004,15.785492,16.129496 -2457362.37500000,0.000000,0.000000,0.000000,0.336635,15.633304,15.969939 -2457362.41666667,0.000000,0.000000,0.000000,0.329399,15.483469,15.812868 -2457362.45833333,0.000000,0.000000,0.000000,0.322294,15.335933,15.658226 -2457362.50000000,0.000000,0.000000,0.000000,0.315318,15.190645,15.505963 -2457362.54166667,0.000000,0.000000,0.000000,0.308471,15.047554,15.356026 -2457362.58333333,0.000000,0.000000,0.000000,0.301751,14.906613,15.208365 -2457362.62500000,0.000000,0.000000,0.000000,0.295157,14.767774,15.062931 -2457362.66666667,0.000000,0.000000,0.000000,0.288687,14.630991,14.919678 -2457362.70833333,0.000000,0.000000,0.000000,0.282339,14.496221,14.778560 -2457362.75000000,0.000000,0.000000,0.000000,0.276113,14.363418,14.639531 -2457362.79166667,0.000000,0.000000,0.000000,0.270006,14.232543,14.502548 -2457362.83333333,0.000000,0.000000,0.000000,0.264016,14.103553,14.367570 -2457362.87500000,0.000000,0.000000,0.000000,0.258143,13.976410,14.234554 -2457362.91666667,0.000000,0.000000,0.000000,0.252385,13.851075,14.103460 -2457362.95833333,0.000000,0.000000,0.000000,0.246740,13.727510,13.974250 -2457363.00000000,0.000000,0.000000,0.000000,0.241207,13.605679,13.846886 -2457363.04166667,0.000000,0.000000,0.000000,0.235784,13.485546,13.721330 -2457363.08333333,0.000000,0.000000,0.000000,0.230468,13.367078,13.597546 -2457363.12500000,0.000000,0.000000,0.000000,0.225260,13.250240,13.475500 -2457363.16666667,0.000000,0.000000,0.000000,0.220157,13.135000,13.355157 -2457363.20833333,0.000000,0.000000,0.000000,0.215157,13.021326,13.236483 -2457363.25000000,0.000000,0.000000,0.000000,0.210259,12.909187,13.119445 -2457363.29166667,0.000000,0.000000,0.000000,0.205461,12.798553,13.004014 -2457363.33333333,0.000000,0.000000,0.000000,0.200762,12.689394,12.890156 -2457363.37500000,0.000000,0.000000,0.000000,0.196160,12.581682,12.777842 -2457363.41666667,0.000000,0.000000,0.000000,0.191653,12.475390,12.667043 -2457363.45833333,0.000000,0.000000,0.000000,0.187240,12.370490,12.557730 -2457363.50000000,0.000000,0.000000,0.000000,0.182920,12.266955,12.449875 -2457363.54166667,0.000000,0.000000,0.000000,0.178690,12.164760,12.343450 -2457363.58333333,0.000000,0.000000,0.000000,0.174549,12.063879,12.238428 -2457363.62500000,0.000000,0.000000,0.000000,0.170496,11.964288,12.134785 -2457363.66666667,0.000000,0.000000,0.000000,0.166529,11.865964,12.032493 -2457363.70833333,0.000000,0.000000,0.000000,0.162647,11.768882,11.931529 -2457363.75000000,0.000000,0.000000,0.000000,0.158848,11.673021,11.831868 -2457363.79166667,0.000000,0.000000,0.000000,0.155130,11.578357,11.733487 -2457363.83333333,0.000000,0.000000,0.000000,0.151492,11.484869,11.636362 -2457363.87500000,0.000000,0.000000,0.000000,0.147933,11.392537,11.540470 -2457363.91666667,0.000000,0.000000,0.000000,0.144451,11.301339,11.445790 -2457363.95833333,0.000000,0.000000,0.000000,0.141045,11.211254,11.352300 -2457364.00000000,0.000000,0.000000,0.000000,0.137714,11.122265,11.259979 -2457364.04166667,0.000000,0.000000,0.000000,0.134455,11.034351,11.168806 -2457364.08333333,0.000000,0.000000,0.000000,0.131268,10.947493,11.078761 -2457364.12500000,0.000000,0.000000,0.000000,0.128151,10.861673,10.989824 -2457364.16666667,0.000000,0.000000,0.000000,0.125103,10.776874,10.901977 -2457364.20833333,0.000000,0.000000,0.000000,0.122122,10.693077,10.815199 -2457364.25000000,0.000000,0.000000,0.000000,0.119208,10.610266,10.729474 -2457364.29166667,0.000000,0.000000,0.000000,0.116359,10.528424,10.644782 -2457364.33333333,0.000000,0.000000,0.000000,0.113573,10.447534,10.561107 -2457364.37500000,0.000000,0.000000,0.000000,0.110850,10.367581,10.478431 -2457364.41666667,0.000000,0.000000,0.000000,0.108188,10.288549,10.396736 -2457364.45833333,0.000000,0.000000,0.000000,0.105586,10.210422,10.316008 -2457364.50000000,0.000000,0.000000,0.000000,0.103043,10.133186,10.236229 -2457364.54166667,0.000000,0.000000,0.000000,0.100557,10.056826,10.157383 -2457364.58333333,0.000000,0.000000,0.000000,0.098127,9.981328,10.079456 -2457364.62500000,0.000000,0.000000,0.000000,0.095753,9.906678,10.002431 -2457364.66666667,0.000000,0.000000,0.000000,0.093433,9.832862,9.926295 -2457364.70833333,0.000000,0.000000,0.000000,0.091166,9.759866,9.851033 -2457364.75000000,0.000000,0.000000,0.000000,0.088951,9.687678,9.776630 -2457364.79166667,0.000000,0.000000,0.000000,0.086787,9.616285,9.703072 -2457364.83333333,0.000000,0.000000,0.000000,0.084673,9.545674,9.630347 -2457364.87500000,0.000000,0.000000,0.000000,0.082607,9.475833,9.558440 -2457364.91666667,0.000000,0.000000,0.000000,0.080589,9.406750,9.487339 -2457364.95833333,0.000000,0.000000,0.000000,0.078618,9.338412,9.417030 -2457365.00000000,0.000000,0.000000,0.000000,0.076693,9.270809,9.347502 -2457365.04166667,0.000000,0.000000,0.000000,0.074812,9.203930,9.278742 -2457365.08333333,0.000000,0.000000,0.000000,0.072975,9.137762,9.210737 -2457365.12500000,0.000000,0.000000,0.000000,0.071181,9.072296,9.143477 -2457365.16666667,0.000000,0.000000,0.000000,0.069429,9.007520,9.076949 -2457365.20833333,0.000000,0.000000,0.000000,0.067718,8.943424,9.011143 -2457365.25000000,0.000000,0.000000,0.000000,0.066048,8.879999,8.946046 -2457365.29166667,0.000000,0.000000,0.000000,0.064416,8.817233,8.881649 -2457365.33333333,0.000000,0.000000,0.000000,0.062823,8.755117,8.817940 -2457365.37500000,0.000000,0.000000,0.000000,0.061268,8.693642,8.754909 -2457365.41666667,0.000000,0.000000,0.000000,0.059749,8.632798,8.692547 -2457365.45833333,0.000000,0.000000,0.000000,0.058266,8.572575,8.630842 -2457365.50000000,0.000000,0.000000,0.000000,0.056819,8.512965,8.569784 -2457365.54166667,0.000000,0.000000,0.000000,0.055406,8.453960,8.509365 -2457365.58333333,0.000000,0.000000,0.000000,0.054026,8.395549,8.449575 -2457365.62500000,0.000000,0.000000,0.000000,0.052680,8.337725,8.390405 -2457365.66666667,0.000000,0.000000,0.000000,0.051366,8.280479,8.331844 -2457365.70833333,0.000000,0.000000,0.000000,0.050083,8.223802,8.273885 -2457365.75000000,0.000000,0.000000,0.000000,0.048831,8.167688,8.216519 -2457365.79166667,0.000000,0.000000,0.000000,0.047608,8.112128,8.159736 -2457365.83333333,0.000000,0.000000,0.000000,0.046416,8.057114,8.103529 -2457365.87500000,0.000000,0.000000,0.000000,0.045252,8.002638,8.047890 -2457365.91666667,0.000000,0.000000,0.000000,0.044116,7.948693,7.992809 -2457365.95833333,0.000000,0.000000,0.000000,0.043007,7.895272,7.938279 -2457366.00000000,0.000000,0.000000,0.000000,0.041926,7.842367,7.884293 -2457366.04166667,0.000000,0.000000,0.000000,0.040870,7.789972,7.830842 -2457366.08333333,0.000000,0.000000,0.000000,0.039840,7.738079,7.777919 -2457366.12500000,0.000000,0.000000,0.000000,0.038835,7.686682,7.725517 -2457366.16666667,0.000000,0.000000,0.000000,0.037855,7.635773,7.673628 -2457366.20833333,0.000000,0.000000,0.000000,0.036898,7.585347,7.622245 -2457366.25000000,0.000000,0.000000,0.000000,0.035965,7.535396,7.571361 -2457366.29166667,0.000000,0.000000,0.000000,0.035054,7.485915,7.520970 -2457366.33333333,0.000000,0.000000,0.000000,0.034166,7.436897,7.471063 -2457366.37500000,0.000000,0.000000,0.000000,0.033300,7.388336,7.421636 -2457366.41666667,0.000000,0.000000,0.000000,0.032455,7.340226,7.372681 -2457366.45833333,0.000000,0.000000,0.000000,0.031630,7.292562,7.324192 -2457366.50000000,0.000000,0.000000,0.000000,0.030826,7.245336,7.276162 -2457366.54166667,0.000000,0.000000,0.000000,0.030041,7.198544,7.228585 -2457366.58333333,0.000000,0.000000,0.000000,0.029276,7.152180,7.181456 -2457366.62500000,0.000000,0.000000,0.000000,0.028530,7.106238,7.134768 -2457366.66666667,0.000000,0.000000,0.000000,0.027802,7.060714,7.088516 -2457366.70833333,0.000000,0.000000,0.000000,0.027092,7.015601,7.042693 -2457366.75000000,0.000000,0.000000,0.000000,0.026400,6.970894,6.997294 -2457366.79166667,0.000000,0.000000,0.000000,0.025725,6.926588,6.952314 -2457366.83333333,0.000000,0.000000,0.000000,0.025067,6.882679,6.907746 -2457366.87500000,0.000000,0.000000,0.000000,0.024425,6.839161,6.863586 -2457366.91666667,0.000000,0.000000,0.000000,0.023799,6.796029,6.819828 -2457366.95833333,0.000000,0.000000,0.000000,0.023188,6.753279,6.776467 -2457367.00000000,0.000000,0.000000,0.000000,0.022593,6.710905,6.733498 -2457367.04166667,0.000000,0.000000,0.000000,0.022012,6.668904,6.690916 -2457367.08333333,0.000000,0.000000,0.000000,0.021447,6.627270,6.648716 -2457367.12500000,0.000000,0.000000,0.000000,0.020895,6.585998,6.606893 -2457367.16666667,0.000000,0.000000,0.000000,0.020357,6.545086,6.565442 -2457367.20833333,0.000000,0.000000,0.000000,0.019832,6.504527,6.524359 -2457367.25000000,0.000000,0.000000,0.000000,0.019321,6.464318,6.483639 -2457367.29166667,0.000000,0.000000,0.000000,0.018822,6.424455,6.443278 -2457367.33333333,0.000000,0.000000,0.000000,0.018336,6.384933,6.403270 -2457367.37500000,0.000000,0.000000,0.000000,0.017863,6.345749,6.363612 -2457367.41666667,0.000000,0.000000,0.000000,0.017401,6.306898,6.324299 -2457367.45833333,0.000000,0.000000,0.000000,0.016951,6.268376,6.285327 -2457367.50000000,0.000000,0.000000,0.000000,0.016512,6.230180,6.246692 -2457367.54166667,0.000000,0.000000,0.000000,0.016084,6.192306,6.208390 -2457367.58333333,0.000000,0.000000,0.000000,0.015667,6.154749,6.170416 -2457367.62500000,0.000000,0.000000,0.000000,0.015260,6.117506,6.132767 -2457367.66666667,0.000000,0.000000,0.000000,0.014864,6.080574,6.095438 -2457367.70833333,0.000000,0.000000,0.000000,0.014478,6.043949,6.058427 -2457367.75000000,0.000000,0.000000,0.000000,0.014102,6.007626,6.021728 -2457367.79166667,0.000000,0.000000,0.000000,0.013735,5.971604,5.985339 -2457367.83333333,0.000000,0.000000,0.000000,0.013378,5.935878,5.949256 -2457367.87500000,0.000000,0.000000,0.000000,0.013029,5.900445,5.913474 -2457367.91666667,0.000000,0.000000,0.000000,0.012690,5.865302,5.877991 -2457367.95833333,0.000000,0.000000,0.000000,0.012359,5.830444,5.842803 -2457368.00000000,0.000000,0.000000,0.000000,0.012036,5.795870,5.807907 -2457368.04166667,0.000000,0.000000,0.000000,0.011722,5.761576,5.773298 -2457368.08333333,0.000000,0.000000,0.000000,0.011416,5.727558,5.738974 -2457368.12500000,0.000000,0.000000,0.000000,0.011118,5.693814,5.704932 -2457368.16666667,0.000000,0.000000,0.000000,0.010827,5.660341,5.671168 -2457368.20833333,0.000000,0.000000,0.000000,0.010544,5.627135,5.637679 -2457368.25000000,0.000000,0.000000,0.000000,0.010268,5.594194,5.604461 -2457368.29166667,0.000000,0.000000,0.000000,0.009999,5.561514,5.571513 -2457368.33333333,0.000000,0.000000,0.000000,0.009736,5.529094,5.538830 -2457368.37500000,0.000000,0.000000,0.000000,0.009481,5.496929,5.506410 -2457368.41666667,0.000000,0.000000,0.000000,0.009232,5.465017,5.474249 -2457368.45833333,0.000000,0.000000,0.000000,0.008990,5.433356,5.442346 -2457368.50000000,0.000000,0.000000,0.000000,0.008754,5.401942,5.410696 -2457368.54166667,0.000000,0.000000,0.000000,0.008524,5.370774,5.379297 -2457368.58333333,0.000000,0.000000,0.000000,0.008299,5.339848,5.348147 -2457368.62500000,0.000000,0.000000,0.000000,0.008081,5.309161,5.317242 -2457368.66666667,0.000000,0.000000,0.000000,0.007868,5.278712,5.286580 -2457368.70833333,0.000000,0.000000,0.000000,0.007661,5.248497,5.256158 -2457368.75000000,0.000000,0.000000,0.000000,0.007459,5.218515,5.225974 -2457368.79166667,0.000000,0.000000,0.000000,0.007263,5.188762,5.196025 -2457368.83333333,0.000000,0.000000,0.000000,0.007071,5.159237,5.166308 -2457368.87500000,0.000000,0.000000,0.000000,0.006884,5.129937,5.136821 -2457368.91666667,0.000000,0.000000,0.000000,0.006702,5.100859,5.107561 -2457368.95833333,0.000000,0.000000,0.000000,0.006525,5.072001,5.078526 -2457369.00000000,0.000000,0.000000,0.000000,0.006353,5.043361,5.049714 -2457369.04166667,0.000000,0.000000,0.000000,0.006185,5.014937,5.021122 -2457369.08333333,0.000000,0.000000,0.000000,0.006021,4.986727,4.992748 -2457369.12500000,0.000000,0.000000,0.000000,0.005862,4.958727,4.964589 -2457369.16666667,0.000000,0.000000,0.000000,0.005707,4.930937,4.936644 -2457369.20833333,0.000000,0.000000,0.000000,0.005555,4.903354,4.908909 -2457369.25000000,0.000000,0.000000,0.000000,0.005408,4.875975,4.881383 -2457369.29166667,0.000000,0.000000,0.000000,0.005265,4.848800,4.854064 -2457369.33333333,0.000000,0.000000,0.000000,0.005125,4.821825,4.826949 -2457369.37500000,0.000000,0.000000,0.000000,0.004989,4.795048,4.800037 -2457369.41666667,0.000000,0.000000,0.000000,0.004856,4.768468,4.773325 -2457369.45833333,0.000000,0.000000,0.000000,0.004727,4.742083,4.746811 -2457369.50000000,0.000000,0.000000,0.000000,0.004602,4.715891,4.720493 -2457369.54166667,0.000000,0.000000,0.000000,0.004479,4.689890,4.694369 -2457369.58333333,0.000000,0.000000,0.000000,0.004360,4.664077,4.668437 -2457369.62500000,0.000000,0.000000,0.000000,0.004244,4.638452,4.642696 -2457369.66666667,0.000000,0.000000,0.000000,0.004131,4.613011,4.617142 -2457369.70833333,0.000000,0.000000,0.000000,0.004021,4.587755,4.591775 -2457369.75000000,0.000000,0.000000,0.000000,0.003914,4.562679,4.566593 -2457369.79166667,0.000000,0.000000,0.000000,0.003809,4.537784,4.541593 -2457369.83333333,0.000000,0.000000,0.000000,0.003708,4.513066,4.516774 -2457369.87500000,0.000000,0.000000,0.000000,0.003609,4.488525,4.492134 -2457369.91666667,0.000000,0.000000,0.000000,0.003512,4.464159,4.467671 -2457369.95833333,0.000000,0.000000,0.000000,0.003419,4.439965,4.443384 -2457370.00000000,0.000000,0.000000,0.000000,0.003327,4.415943,4.419270 -2457370.04166667,0.000000,0.000000,0.000000,0.003238,4.392090,4.395328 -2457370.08333333,0.000000,0.000000,0.000000,0.003152,4.368405,4.371557 -2457370.12500000,0.000000,0.000000,0.000000,0.003067,4.344887,4.347954 -2457370.16666667,0.000000,0.000000,0.000000,0.002985,4.321533,4.324518 -2457370.20833333,0.000000,0.000000,0.000000,0.002905,4.298343,4.301248 -2457370.25000000,0.000000,0.000000,0.000000,0.002827,4.275314,4.278142 -2457370.29166667,0.000000,0.000000,0.000000,0.002752,4.252446,4.255197 -2457370.33333333,0.000000,0.000000,0.000000,0.002678,4.229736,4.232414 -2457370.37500000,0.000000,0.000000,0.000000,0.002606,4.207183,4.209789 -2457370.41666667,0.000000,0.000000,0.000000,0.002536,4.184786,4.187323 -2457370.45833333,0.000000,0.000000,0.000000,0.002468,4.162544,4.165012 -2457370.50000000,0.000000,0.000000,0.000000,0.002402,4.140454,4.142856 -2457370.54166667,0.000000,0.000000,0.000000,0.002337,4.118516,4.120853 -2457370.58333333,0.000000,0.000000,0.000000,0.002274,4.096727,4.099002 -2457370.62500000,0.000000,0.000000,0.000000,0.002213,4.075087,4.077301 -2457370.66666667,0.000000,0.000000,0.000000,0.002154,4.053595,4.055749 -2457370.70833333,0.000000,0.000000,0.000000,0.002096,4.032248,4.034344 -2457370.75000000,0.000000,0.000000,0.000000,0.002039,4.011046,4.013086 -2457370.79166667,0.000000,0.000000,0.000000,0.001985,3.989987,3.991972 -2457370.83333333,0.000000,0.000000,0.000000,0.001931,3.969071,3.971002 -2457370.87500000,0.000000,0.000000,0.000000,0.001879,3.948295,3.950174 -2457370.91666667,0.000000,0.000000,0.000000,0.001828,3.927658,3.929486 -2457370.95833333,0.000000,0.000000,0.000000,0.001779,3.907159,3.908938 -2457371.00000000,0.000000,0.000000,0.000000,0.001731,3.886798,3.888529 -2457371.04166667,0.000000,0.000000,0.000000,0.001684,3.866572,3.868256 -2457371.08333333,0.000000,0.000000,0.000000,0.001639,3.846480,3.848119 -2457371.12500000,1.500000,0.029641,0.001778,0.001595,3.826759,3.830133 -2457371.16666667,0.000000,0.000000,0.015117,0.001552,3.806932,3.823600 -2457371.20833333,0.000000,0.000000,0.008299,0.002792,3.787235,3.798326 -2457371.25000000,0.000000,0.000000,0.003557,0.003956,3.767668,3.775181 -2457371.29166667,0.000000,0.000000,0.000889,0.005048,3.748230,3.754167 -2457371.33333333,0.100000,0.000134,0.000008,0.006071,3.728935,3.735014 -2457371.37500000,4.300000,0.234917,0.014163,0.007028,3.710402,3.731593 -2457371.41666667,0.000000,0.000000,0.119845,0.008009,3.691338,3.819193 -2457371.45833333,0.000000,0.000000,0.065793,0.012470,3.672399,3.750662 -2457371.50000000,0.000000,0.000000,0.028194,0.016656,3.653583,3.698433 -2457371.54166667,0.000000,0.000000,0.007048,0.020579,3.634889,3.662516 -2457371.58333333,0.000000,0.000000,0.000000,0.024252,3.616316,3.640569 -2457371.62500000,0.000000,0.000000,0.000000,0.027686,3.597863,3.625550 -2457371.66666667,0.000000,0.000000,0.000000,0.030893,3.579530,3.610422 -2457371.70833333,0.000000,0.000000,0.000000,0.033882,3.561314,3.595195 -2457371.75000000,0.000000,0.000000,0.000000,0.036664,3.543215,3.579878 -2457371.79166667,0.000000,0.000000,0.000000,0.039248,3.525232,3.564480 -2457371.83333333,0.000000,0.000000,0.000000,0.041644,3.507364,3.549008 -2457371.87500000,0.000000,0.000000,0.000000,0.043862,3.489610,3.533471 -2457371.91666667,0.000000,0.000000,0.000000,0.045908,3.471969,3.517877 -2457371.95833333,0.000000,0.000000,0.000000,0.047793,3.454440,3.502233 -2457372.00000000,0.000000,0.000000,0.000000,0.049523,3.437022,3.486545 -2457372.04166667,0.000000,0.000000,0.000000,0.051106,3.419714,3.470820 -2457372.08333333,0.000000,0.000000,0.000000,0.052550,3.402516,3.455066 -2457372.12500000,0.000000,0.000000,0.000000,0.053861,3.385426,3.439287 -2457372.16666667,0.000000,0.000000,0.000000,0.055046,3.368444,3.423490 -2457372.20833333,0.000000,0.000000,0.000000,0.056111,3.351568,3.407679 -2457372.25000000,0.000000,0.000000,0.000000,0.057064,3.334798,3.391861 -2457372.29166667,0.000000,0.000000,0.000000,0.057908,3.318132,3.376040 -2457372.33333333,0.000000,0.000000,0.000000,0.058650,3.301571,3.360221 -2457372.37500000,0.000000,0.000000,0.000000,0.059296,3.285112,3.344409 -2457372.41666667,0.000000,0.000000,0.000000,0.059850,3.268756,3.328607 -2457372.45833333,0.000000,0.000000,0.000000,0.060318,3.252502,3.312820 -2457372.50000000,0.000000,0.000000,0.000000,0.060704,3.236348,3.297051 -2457372.54166667,0.000000,0.000000,0.000000,0.061012,3.220293,3.281305 -2457372.58333333,0.000000,0.000000,0.000000,0.061247,3.204338,3.265585 -2457372.62500000,0.000000,0.000000,0.000000,0.061413,3.188481,3.249894 -2457372.66666667,0.000000,0.000000,0.000000,0.061514,3.172721,3.234235 -2457372.70833333,0.000000,0.000000,0.000000,0.061554,3.157058,3.218612 -2457372.75000000,0.000000,0.000000,0.000000,0.061536,3.141490,3.203026 -2457372.79166667,0.000000,0.000000,0.000000,0.061463,3.126018,3.187481 -2457372.83333333,0.000000,0.000000,0.000000,0.061340,3.110640,3.171980 -2457372.87500000,0.000000,0.000000,0.000000,0.061169,3.095356,3.156524 -2457372.91666667,0.000000,0.000000,0.000000,0.060952,3.080164,3.141116 -2457372.95833333,0.000000,0.000000,0.000000,0.060694,3.065064,3.125758 -2457373.00000000,0.000000,0.000000,0.000000,0.060396,3.050056,3.110452 -2457373.04166667,0.000000,0.000000,0.000000,0.060062,3.035138,3.095200 -2457373.08333333,0.000000,0.000000,0.000000,0.059693,3.020310,3.080003 -2457373.12500000,0.000000,0.000000,0.000000,0.059293,3.005572,3.064864 -2457373.16666667,0.000000,0.000000,0.000000,0.058862,2.990921,3.049784 -2457373.20833333,0.000000,0.000000,0.000000,0.058404,2.976359,3.034763 -2457373.25000000,0.000000,0.000000,0.000000,0.057921,2.961884,3.019805 -2457373.29166667,0.000000,0.000000,0.000000,0.057414,2.947495,3.004909 -2457373.33333333,0.000000,0.000000,0.000000,0.056885,2.933192,2.990077 -2457373.37500000,0.000000,0.000000,0.000000,0.056337,2.918974,2.975310 -2457373.41666667,0.000000,0.000000,0.000000,0.055770,2.904840,2.960610 -2457373.45833333,0.000000,0.000000,0.000000,0.055186,2.890790,2.945976 -2457373.50000000,0.000000,0.000000,0.000000,0.054587,2.876824,2.931411 -2457373.54166667,0.000000,0.000000,0.000000,0.053975,2.862940,2.916914 -2457373.58333333,0.000000,0.000000,0.000000,0.053350,2.849137,2.902487 -2457373.62500000,0.000000,0.000000,0.000000,0.052714,2.835416,2.888130 -2457373.66666667,0.500000,0.003337,0.000200,0.052068,2.821850,2.874119 -2457373.70833333,1.000000,0.013261,0.002498,0.051414,2.808437,2.862348 -2457373.75000000,3.000000,0.116266,0.014674,0.051185,2.795386,2.861244 -2457373.79166667,1.000000,0.013261,0.064205,0.051783,2.782127,2.898115 -2457373.83333333,2.200000,0.063178,0.044800,0.054815,2.769117,2.868732 -2457373.87500000,0.000000,0.000000,0.050284,0.058478,2.755864,2.864626 -2457373.91666667,0.000000,0.000000,0.022769,0.063736,2.742688,2.829193 -2457373.95833333,0.000000,0.000000,0.007979,0.068625,2.729589,2.806193 -2457374.00000000,3.000000,0.116266,0.008871,0.073164,2.716994,2.799029 -2457374.04166667,4.400000,0.245658,0.074035,0.077368,2.704660,2.856063 -2457374.08333333,4.400000,0.245658,0.172580,0.083769,2.692397,2.948745 -2457374.12500000,1.600000,0.033680,0.210043,0.093339,2.679820,2.983202 -2457374.16666667,2.300000,0.068962,0.123066,0.105881,2.667414,2.896362 -2457374.20833333,3.800000,0.184634,0.092528,0.118953,2.655284,2.866765 -2457374.25000000,0.000000,0.000000,0.124884,0.133091,2.642689,2.900664 -2457374.29166667,0.000000,0.000000,0.060983,0.149426,2.630166,2.840575 -2457374.33333333,0.000000,0.000000,0.024225,0.164661,2.617715,2.806601 -2457374.37500000,0.000000,0.000000,0.005539,0.178847,2.605335,2.789721 -2457374.41666667,0.000000,0.000000,0.000000,0.192034,2.593026,2.785060 -2457374.45833333,0.000000,0.000000,0.000000,0.204269,2.580786,2.785055 -2457374.50000000,0.000000,0.000000,0.000000,0.215597,2.568616,2.784214 -2457374.54166667,0.000000,0.000000,0.000000,0.226062,2.556516,2.782577 -2457374.58333333,0.000000,0.000000,0.000000,0.235704,2.544484,2.780188 -2457374.62500000,0.000000,0.000000,0.000000,0.244565,2.532520,2.777084 -2457374.66666667,0.000000,0.000000,0.000000,0.252681,2.520623,2.773305 -2457374.70833333,0.000000,0.000000,0.000000,0.260091,2.508794,2.768885 -2457374.75000000,0.000000,0.000000,0.000000,0.266828,2.497031,2.763860 -2457374.79166667,0.000000,0.000000,0.000000,0.272928,2.485335,2.758263 -2457374.83333333,0.000000,0.000000,0.000000,0.278421,2.473704,2.752125 -2457374.87500000,0.000000,0.000000,0.000000,0.283339,2.462138,2.745477 -2457374.91666667,0.000000,0.000000,0.000000,0.287711,2.450637,2.738348 -2457374.95833333,0.000000,0.000000,0.000000,0.291566,2.439200,2.730767 -2457375.00000000,0.000000,0.000000,0.000000,0.294931,2.427827,2.722758 -2457375.04166667,0.000000,0.000000,0.000000,0.297831,2.416518,2.714349 -2457375.08333333,0.000000,0.000000,0.000000,0.300292,2.405271,2.705563 -2457375.12500000,0.000000,0.000000,0.000000,0.302338,2.394087,2.696424 -2457375.16666667,0.000000,0.000000,0.000000,0.303990,2.382965,2.686955 -2457375.20833333,0.000000,0.000000,0.000000,0.305271,2.371904,2.677175 -2457375.25000000,0.000000,0.000000,0.000000,0.306202,2.360905,2.667106 -2457375.29166667,0.000000,0.000000,0.000000,0.306802,2.349966,2.656767 -2457375.33333333,0.000000,0.000000,0.000000,0.307090,2.339087,2.646177 -2457375.37500000,0.000000,0.000000,0.000000,0.307085,2.328269,2.635354 -2457375.41666667,0.000000,0.000000,0.000000,0.306804,2.317510,2.624314 -2457375.45833333,0.000000,0.000000,0.000000,0.306263,2.306810,2.613073 -2457375.50000000,0.000000,0.000000,0.000000,0.305479,2.296168,2.601647 -2457375.54166667,0.000000,0.000000,0.000000,0.304466,2.285585,2.590051 -2457375.58333333,0.000000,0.000000,0.000000,0.303239,2.275060,2.578299 -2457375.62500000,0.000000,0.000000,0.000000,0.301812,2.264592,2.566403 -2457375.66666667,0.000000,0.000000,0.000000,0.300197,2.254181,2.554378 -2457375.70833333,0.000000,0.000000,0.000000,0.298408,2.243826,2.542234 -2457375.75000000,0.000000,0.000000,0.000000,0.296455,2.233528,2.529984 -2457375.79166667,0.000000,0.000000,0.000000,0.294352,2.223286,2.517638 -2457375.83333333,0.000000,0.000000,0.000000,0.292107,2.213099,2.505207 -2457375.87500000,0.000000,0.000000,0.000000,0.289733,2.202968,2.492701 -2457375.91666667,0.000000,0.000000,0.000000,0.287238,2.192891,2.480129 -2457375.95833333,0.000000,0.000000,0.000000,0.284632,2.182868,2.467500 -2457376.00000000,0.000000,0.000000,0.000000,0.281924,2.172900,2.454824 -2457376.04166667,0.000000,0.000000,0.000000,0.279122,2.162985,2.442107 -2457376.08333333,0.000000,0.000000,0.000000,0.276234,2.153123,2.429358 -2457376.12500000,0.000000,0.000000,0.000000,0.273269,2.143315,2.416584 -2457376.16666667,0.000000,0.000000,0.000000,0.270233,2.133558,2.403791 -2457376.20833333,0.000000,0.000000,0.000000,0.267133,2.123854,2.390987 -2457376.25000000,0.000000,0.000000,0.000000,0.263976,2.114202,2.378178 -2457376.29166667,0.000000,0.000000,0.000000,0.260768,2.104602,2.365369 -2457376.33333333,0.000000,0.000000,0.000000,0.257515,2.095052,2.352567 -2457376.37500000,0.000000,0.000000,0.000000,0.254222,2.085554,2.339775 -2457376.41666667,0.000000,0.000000,0.000000,0.250895,2.076105,2.327000 -2457376.45833333,0.000000,0.000000,0.000000,0.247539,2.066707,2.314246 -2457376.50000000,0.000000,0.000000,0.000000,0.244158,2.057359,2.301517 -2457376.54166667,0.000000,0.000000,0.000000,0.240758,2.048060,2.288818 -2457376.58333333,0.000000,0.000000,0.000000,0.237342,2.038811,2.276152 -2457376.62500000,0.000000,0.000000,0.000000,0.233914,2.029610,2.263524 -2457376.66666667,0.000000,0.000000,0.000000,0.230478,2.020457,2.250935 -2457376.70833333,0.000000,0.000000,0.000000,0.227038,2.011353,2.238391 -2457376.75000000,0.000000,0.000000,0.000000,0.223597,2.002297,2.225894 -2457376.79166667,0.000000,0.000000,0.000000,0.220157,1.993288,2.213446 -2457376.83333333,0.000000,0.000000,0.000000,0.216723,1.984327,2.201051 -2457376.87500000,0.000000,0.000000,0.000000,0.213297,1.975413,2.188710 -2457376.91666667,0.000000,0.000000,0.000000,0.209882,1.966545,2.176427 -2457376.95833333,0.000000,0.000000,0.000000,0.206479,1.957723,2.164202 -2457377.00000000,0.000000,0.000000,0.000000,0.203092,1.948948,2.152040 -2457377.04166667,0.000000,0.000000,0.000000,0.199722,1.940218,2.139940 -2457377.08333333,0.000000,0.000000,0.000000,0.196371,1.931534,2.127905 -2457377.12500000,0.000000,0.000000,0.000000,0.193042,1.922894,2.115937 -2457377.16666667,0.000000,0.000000,0.000000,0.189736,1.914300,2.104036 -2457377.20833333,0.000000,0.000000,0.000000,0.186454,1.905750,2.092205 -2457377.25000000,0.000000,0.000000,0.000000,0.183199,1.897245,2.080444 -2457377.29166667,0.000000,0.000000,0.000000,0.179971,1.888783,2.068754 -2457377.33333333,0.000000,0.000000,0.000000,0.176772,1.880366,2.057138 -2457377.37500000,0.000000,0.000000,0.000000,0.173603,1.871991,2.045594 -2457377.41666667,0.000000,0.000000,0.000000,0.170465,1.863660,2.034125 -2457377.45833333,0.000000,0.000000,0.000000,0.167359,1.855372,2.022730 -2457377.50000000,0.000000,0.000000,0.000000,0.164285,1.847126,2.011412 -2457377.54166667,0.000000,0.000000,0.000000,0.161246,1.838923,2.000169 -2457377.58333333,0.000000,0.000000,0.000000,0.158241,1.830762,1.989003 -2457377.62500000,0.000000,0.000000,0.000000,0.155272,1.822642,1.977914 -2457377.66666667,0.000000,0.000000,0.000000,0.152338,1.814564,1.966902 -2457377.70833333,0.000000,0.000000,0.000000,0.149440,1.806528,1.955968 -2457377.75000000,0.000000,0.000000,0.000000,0.146579,1.798532,1.945111 -2457377.79166667,0.000000,0.000000,0.000000,0.143755,1.790577,1.934333 -2457377.83333333,0.000000,0.000000,0.000000,0.140969,1.782663,1.923632 -2457377.87500000,0.000000,0.000000,0.000000,0.138221,1.774789,1.913010 -2457377.91666667,0.000000,0.000000,0.000000,0.135510,1.766955,1.902465 -2457377.95833333,0.000000,0.000000,0.000000,0.132838,1.759160,1.891999 -2457378.00000000,0.000000,0.000000,0.000000,0.130205,1.751406,1.881610 -2457378.04166667,0.000000,0.000000,0.000000,0.127609,1.743690,1.871299 -2457378.08333333,0.000000,0.000000,0.000000,0.125052,1.736014,1.861066 -2457378.12500000,0.000000,0.000000,0.000000,0.122534,1.728376,1.850910 -2457378.16666667,0.000000,0.000000,0.000000,0.120055,1.720777,1.840831 -2457378.20833333,0.000000,0.000000,0.000000,0.117613,1.713216,1.830829 -2457378.25000000,0.000000,0.000000,0.000000,0.115211,1.705693,1.820904 -2457378.29166667,0.000000,0.000000,0.000000,0.112846,1.698208,1.811054 -2457378.33333333,0.000000,0.000000,0.000000,0.110520,1.690761,1.801281 -2457378.37500000,0.000000,0.000000,0.000000,0.108232,1.683351,1.791582 -2457378.41666667,0.000000,0.000000,0.000000,0.105981,1.675978,1.781959 -2457378.45833333,0.000000,0.000000,0.000000,0.103768,1.668642,1.772410 -2457378.50000000,0.000000,0.000000,0.000000,0.101593,1.661343,1.762936 -2457378.54166667,0.000000,0.000000,0.000000,0.099455,1.654080,1.753534 -2457378.58333333,0.000000,0.000000,0.000000,0.097353,1.646853,1.744206 -2457378.62500000,0.000000,0.000000,0.000000,0.095288,1.639663,1.734951 -2457378.66666667,0.000000,0.000000,0.000000,0.093260,1.632508,1.725768 -2457378.70833333,0.000000,0.000000,0.000000,0.091267,1.625389,1.716656 -2457378.75000000,1.300000,0.022322,0.001339,0.089310,1.618478,1.709128 -2457378.79166667,6.200000,0.476871,0.039997,0.087388,1.612205,1.739590 -2457378.83333333,6.500000,0.522193,0.280786,0.086616,1.605997,1.973399 -2457378.87500000,7.300000,0.652195,0.441653,0.090801,1.599910,2.132364 -2457378.91666667,37.700001,12.677781,1.297394,0.099870,1.596340,2.993605 -2457378.95833333,9.700000,1.118671,6.792372,0.114110,1.590562,8.497044 -2457379.00000000,2.500000,0.081265,4.219106,0.149224,1.583978,5.952307 -2457379.04166667,1.300000,0.022322,1.896912,0.189611,1.577272,3.663795 -2457379.08333333,0.800000,0.008510,0.549223,0.229550,1.570533,2.349306 -2457379.12500000,0.000000,0.000000,0.053902,0.268036,1.563721,1.885659 -2457379.16666667,0.000000,0.000000,0.007499,0.304696,1.556942,1.869137 -2457379.20833333,0.000000,0.000000,0.001691,0.338916,1.550196,1.890803 -2457379.25000000,0.000000,0.000000,0.000255,0.370808,1.543484,1.914547 -2457379.29166667,0.000000,0.000000,0.000000,0.400482,1.536804,1.937286 -2457379.33333333,0.000000,0.000000,0.000000,0.428041,1.530157,1.958198 -2457379.37500000,0.000000,0.000000,0.000000,0.453586,1.523543,1.977130 -2457379.41666667,0.000000,0.000000,0.000000,0.477214,1.516961,1.994175 -2457379.45833333,0.000000,0.000000,0.000000,0.499015,1.510412,2.009427 -2457379.50000000,0.000000,0.000000,0.000000,0.519079,1.503894,2.022973 -2457379.54166667,0.000000,0.000000,0.000000,0.537490,1.497408,2.034898 -2457379.58333333,0.000000,0.000000,0.000000,0.554330,1.490954,2.045283 -2457379.62500000,0.000000,0.000000,0.000000,0.569675,1.484531,2.054206 -2457379.66666667,0.000000,0.000000,0.000000,0.583600,1.478139,2.061739 -2457379.70833333,0.000000,0.000000,0.000000,0.596177,1.471779,2.067956 -2457379.75000000,0.000000,0.000000,0.000000,0.607475,1.465449,2.072924 -2457379.79166667,0.000000,0.000000,0.000000,0.617558,1.459150,2.076708 -2457379.83333333,0.000000,0.000000,0.000000,0.626489,1.452881,2.079370 -2457379.87500000,0.000000,0.000000,0.000000,0.634329,1.446643,2.080972 -2457379.91666667,0.000000,0.000000,0.000000,0.641134,1.440435,2.081569 -2457379.95833333,0.000000,0.000000,0.000000,0.646960,1.434257,2.081217 -2457380.00000000,0.000000,0.000000,0.000000,0.651860,1.428109,2.079969 -2457380.04166667,3.900000,0.194231,0.011654,0.655883,1.422485,2.090022 -2457380.08333333,1.200000,0.019045,0.100201,0.659078,1.416551,2.175830 -2457380.12500000,11.600000,1.564492,0.157967,0.664723,1.411825,2.234515 -2457380.16666667,1.800000,0.042515,0.829082,0.670466,1.406019,2.905567 -2457380.20833333,0.900000,0.010756,0.468498,0.684016,1.400125,2.552639 -2457380.25000000,0.000000,0.000000,0.205700,0.697682,1.394141,2.297523 -2457380.29166667,10.000000,1.184711,0.126131,0.710689,1.389358,2.226177 -2457380.33333333,1.800000,0.042515,0.609320,0.722287,1.383659,2.715266 -2457380.37500000,0.600000,0.004799,0.354012,0.740236,1.377833,2.472081 -2457380.41666667,2.200000,0.063178,0.160308,0.757990,1.372240,2.290538 -2457380.45833333,0.700000,0.006524,0.074599,0.774532,1.366481,2.215613 -2457380.50000000,3.800000,0.184634,0.033946,0.791263,1.361138,2.186347 -2457380.54166667,1.400000,0.025855,0.105267,0.806898,1.355522,2.267687 -2457380.58333333,2.200000,0.063178,0.071352,0.824034,1.350035,2.245421 -2457380.62500000,1.600000,0.033680,0.063833,0.840615,1.344498,2.248945 -2457380.66666667,1.600000,0.033680,0.045529,0.857322,1.338987,2.241838 -2457380.70833333,0.000000,0.000000,0.034964,0.873636,1.333295,2.241895 -2457380.75000000,0.000000,0.000000,0.015367,0.889566,1.327630,2.232563 -2457380.79166667,0.000000,0.000000,0.005052,0.903753,1.321991,2.230797 -2457380.83333333,0.000000,0.000000,0.001010,0.916289,1.316380,2.233679 -2457380.87500000,0.000000,0.000000,0.000000,0.927259,1.310795,2.238054 -2457380.91666667,0.000000,0.000000,0.000000,0.936746,1.305237,2.241983 -2457380.95833333,0.000000,0.000000,0.000000,0.944830,1.299704,2.244534 -2457381.00000000,0.000000,0.000000,0.000000,0.951585,1.294198,2.245783 -2457381.04166667,0.000000,0.000000,0.000000,0.957085,1.288718,2.245803 -2457381.08333333,0.000000,0.000000,0.000000,0.961398,1.283264,2.244662 -2457381.12500000,0.000000,0.000000,0.000000,0.964591,1.277836,2.242427 -2457381.16666667,0.000000,0.000000,0.000000,0.966727,1.272433,2.239160 -2457381.20833333,0.000000,0.000000,0.000000,0.967867,1.267055,2.234923 -2457381.25000000,0.000000,0.000000,0.000000,0.968069,1.261703,2.229772 -2457381.29166667,0.000000,0.000000,0.000000,0.967387,1.256376,2.223763 -2457381.33333333,0.000000,0.000000,0.000000,0.965875,1.251074,2.216949 -2457381.37500000,0.000000,0.000000,0.000000,0.963583,1.245797,2.209380 -2457381.41666667,0.000000,0.000000,0.000000,0.960559,1.240544,2.201104 -2457381.45833333,0.600000,0.004799,0.000288,0.956850,1.235394,2.192532 -2457381.50000000,1.500000,0.029641,0.004226,0.952498,1.230383,2.187108 -2457381.54166667,0.100000,0.000134,0.016469,0.948066,1.225216,2.189751 -2457381.58333333,0.000000,0.000000,0.008944,0.944324,1.220060,2.173328 -2457381.62500000,0.000000,0.000000,0.003738,0.940041,1.214927,2.158707 -2457381.66666667,0.400000,0.002139,0.001034,0.935164,1.209871,2.146069 -2457381.70833333,0.000000,0.000000,0.001095,0.929735,1.204786,2.135616 -2457381.75000000,0.200000,0.000536,0.000631,0.924136,1.199751,2.124518 -2457381.79166667,0.000000,0.000000,0.000530,0.918036,1.194714,2.113280 -2457381.83333333,0.200000,0.000536,0.000246,0.911645,1.189725,2.101617 -2457381.87500000,0.000000,0.000000,0.000338,0.904811,1.184734,2.089883 -2457381.91666667,0.000000,0.000000,0.000166,0.897740,1.179766,2.077673 -2457381.95833333,0.000000,0.000000,0.000064,0.890280,1.174821,2.065165 -2457382.00000000,0.000000,0.000000,0.000016,0.882459,1.169899,2.052374 -2457382.04166667,0.000000,0.000000,0.000000,0.874306,1.165000,2.039306 -2457382.08333333,0.000000,0.000000,0.000000,0.865848,1.160123,2.025971 -2457382.12500000,0.000000,0.000000,0.000000,0.857110,1.155269,2.012379 -2457382.16666667,0.000000,0.000000,0.000000,0.848118,1.150437,1.998555 -2457382.20833333,0.000000,0.000000,0.000000,0.838893,1.145628,1.984522 -2457382.25000000,0.000000,0.000000,0.000000,0.829459,1.140841,1.970300 -2457382.29166667,0.000000,0.000000,0.000000,0.819836,1.136076,1.955912 -2457382.33333333,0.000000,0.000000,0.000000,0.810044,1.131332,1.941376 -2457382.37500000,0.000000,0.000000,0.000000,0.800101,1.126611,1.926712 -2457382.41666667,0.000000,0.000000,0.000000,0.790026,1.121911,1.911938 -2457382.45833333,0.000000,0.000000,0.000000,0.779836,1.117233,1.897069 -2457382.50000000,0.000000,0.000000,0.000000,0.769546,1.112576,1.882122 -2457382.54166667,0.000000,0.000000,0.000000,0.759171,1.107941,1.867112 -2457382.58333333,1.500000,0.029641,0.001778,0.748726,1.103517,1.854022 -2457382.62500000,0.000000,0.000000,0.015117,0.738225,1.098923,1.852265 -2457382.66666667,0.000000,0.000000,0.008299,0.728962,1.094351,1.831612 -2457382.70833333,0.000000,0.000000,0.003557,0.719590,1.089799,1.812946 -2457382.75000000,0.000000,0.000000,0.000889,0.710125,1.085267,1.796281 -2457382.79166667,0.000000,0.000000,0.000000,0.700580,1.080757,1.781337 -2457382.83333333,0.000000,0.000000,0.000000,0.690968,1.076267,1.767235 -2457382.87500000,0.000000,0.000000,0.000000,0.681303,1.071797,1.753100 -2457382.91666667,0.000000,0.000000,0.000000,0.671596,1.067348,1.738944 -2457382.95833333,0.000000,0.000000,0.000000,0.661858,1.062919,1.724777 -2457383.00000000,0.000000,0.000000,0.000000,0.652099,1.058511,1.710610 -2457383.04166667,0.000000,0.000000,0.000000,0.642331,1.054122,1.696453 -2457383.08333333,0.000000,0.000000,0.000000,0.632562,1.049753,1.682315 -2457383.12500000,0.000000,0.000000,0.000000,0.622801,1.045404,1.668205 -2457383.16666667,0.000000,0.000000,0.000000,0.613056,1.041075,1.654131 -2457383.20833333,0.000000,0.000000,0.000000,0.603336,1.036765,1.640101 -2457383.25000000,0.000000,0.000000,0.000000,0.593646,1.032475,1.626122 -2457383.29166667,0.000000,0.000000,0.000000,0.583996,1.028205,1.612200 -2457383.33333333,0.000000,0.000000,0.000000,0.574390,1.023953,1.598343 -2457383.37500000,0.000000,0.000000,0.000000,0.564834,1.019721,1.584555 -2457383.41666667,0.000000,0.000000,0.000000,0.555336,1.015508,1.570843 -2457383.45833333,0.000000,0.000000,0.000000,0.545898,1.011314,1.557212 -2457383.50000000,0.000000,0.000000,0.000000,0.536528,1.007139,1.543666 -2457383.54166667,0.000000,0.000000,0.000000,0.527228,1.002982,1.530210 -2457383.58333333,0.000000,0.000000,0.000000,0.518003,0.998845,1.516847 -2457383.62500000,0.000000,0.000000,0.000000,0.508857,0.994726,1.503583 -2457383.66666667,0.000000,0.000000,0.000000,0.499794,0.990625,1.490419 -2457383.70833333,0.000000,0.000000,0.000000,0.490816,0.986543,1.477359 -2457383.75000000,0.000000,0.000000,0.000000,0.481928,0.982479,1.464407 -2457383.79166667,0.000000,0.000000,0.000000,0.473131,0.978434,1.451564 -2457383.83333333,0.000000,0.000000,0.000000,0.464428,0.974406,1.438834 -2457383.87500000,0.000000,0.000000,0.000000,0.455821,0.970397,1.426218 -2457383.91666667,0.000000,0.000000,0.000000,0.447313,0.966406,1.413719 -2457383.95833333,0.000000,0.000000,0.000000,0.438905,0.962432,1.401337 -2457384.00000000,0.000000,0.000000,0.000000,0.430599,0.958476,1.389075 -2457384.04166667,0.000000,0.000000,0.000000,0.422396,0.954538,1.376935 -2457384.08333333,1.500000,0.029641,0.001778,0.414298,0.950806,1.366883 -2457384.12500000,4.300000,0.234917,0.029212,0.406306,0.947422,1.382940 -2457384.16666667,0.000000,0.000000,0.128107,0.399703,0.943533,1.471342 -2457384.20833333,0.000000,0.000000,0.069334,0.396674,0.939661,1.405669 -2457384.25000000,0.000000,0.000000,0.029079,0.393468,0.935807,1.358354 -2457384.29166667,0.000000,0.000000,0.007048,0.390097,0.931970,1.329115 -2457384.33333333,0.000000,0.000000,0.000000,0.386574,0.928150,1.314724 -2457384.37500000,0.000000,0.000000,0.000000,0.382911,0.924347,1.307259 -2457384.41666667,0.000000,0.000000,0.000000,0.379121,0.920561,1.299682 -2457384.45833333,0.000000,0.000000,0.000000,0.375213,0.916792,1.292005 -2457384.50000000,0.000000,0.000000,0.000000,0.371198,0.913040,1.284238 -2457384.54166667,0.000000,0.000000,0.000000,0.367087,0.909304,1.276390 -2457384.58333333,0.000000,0.000000,0.000000,0.362888,0.905584,1.268472 -2457384.62500000,0.000000,0.000000,0.000000,0.358610,0.901882,1.260492 -2457384.66666667,0.000000,0.000000,0.000000,0.354263,0.898195,1.252458 -2457384.70833333,0.000000,0.000000,0.000000,0.349854,0.894525,1.244379 -2457384.75000000,0.000000,0.000000,0.000000,0.345390,0.890871,1.236261 -2457384.79166667,0.000000,0.000000,0.000000,0.340880,0.887233,1.228113 -2457384.83333333,0.000000,0.000000,0.000000,0.336329,0.883612,1.219941 -2457384.87500000,0.000000,0.000000,0.000000,0.331745,0.880006,1.211750 -2457384.91666667,0.000000,0.000000,0.000000,0.327133,0.876416,1.203549 -2457384.95833333,0.000000,0.000000,0.000000,0.322499,0.872842,1.195341 -2457385.00000000,0.000000,0.000000,0.000000,0.317848,0.869284,1.187132 -2457385.04166667,0.000000,0.000000,0.000000,0.313187,0.865741,1.178928 -2457385.08333333,0.000000,0.000000,0.000000,0.308519,0.862214,1.170733 -2457385.12500000,0.000000,0.000000,0.000000,0.303849,0.858703,1.162552 -2457385.16666667,0.000000,0.000000,0.000000,0.299182,0.855206,1.154388 -2457385.20833333,0.000000,0.000000,0.000000,0.294521,0.851726,1.146247 -2457385.25000000,0.000000,0.000000,0.000000,0.289871,0.848260,1.138131 -2457385.29166667,0.000000,0.000000,0.000000,0.285235,0.844810,1.130045 -2457385.33333333,0.000000,0.000000,0.000000,0.280616,0.841375,1.121990 -2457385.37500000,0.000000,0.000000,0.000000,0.276017,0.837955,1.113972 -2457385.41666667,0.000000,0.000000,0.000000,0.271442,0.834550,1.105991 -2457385.45833333,0.000000,0.000000,0.000000,0.266893,0.831159,1.098052 -2457385.50000000,0.000000,0.000000,0.000000,0.262372,0.827784,1.090156 -2457385.54166667,0.000000,0.000000,0.000000,0.257883,0.824423,1.082306 -2457385.58333333,0.000000,0.000000,0.000000,0.253426,0.821078,1.074504 -2457385.62500000,0.000000,0.000000,0.000000,0.249005,0.817746,1.066752 -2457385.66666667,0.400000,0.002139,0.000128,0.244622,0.814480,1.059230 -2457385.70833333,0.100000,0.000134,0.001099,0.240277,0.811190,1.052565 -2457385.75000000,2.200000,0.063178,0.004458,0.236319,0.808171,1.048949 -2457385.79166667,2.300000,0.068962,0.036653,0.232470,0.805178,1.074301 -2457385.83333333,5.100000,0.327134,0.072569,0.230503,0.802518,1.105589 -2457385.87500000,5.900000,0.433452,0.219740,0.230389,0.799957,1.250086 -2457385.91666667,1.500000,0.029641,0.324607,0.234237,0.796903,1.355747 -2457385.95833333,1.100000,0.016025,0.178770,0.242413,0.793813,1.214996 -2457386.00000000,0.000000,0.000000,0.078300,0.251170,0.790601,1.120071 -2457386.04166667,0.000000,0.000000,0.021047,0.260129,0.787402,1.068578 -2457386.08333333,0.000000,0.000000,0.002812,0.268324,0.784217,1.055353 -2457386.12500000,0.000000,0.000000,0.000481,0.275793,0.781046,1.057320 -2457386.16666667,0.000000,0.000000,0.000000,0.282573,0.777889,1.060462 -2457386.20833333,0.000000,0.000000,0.000000,0.288697,0.774745,1.063443 -2457386.25000000,0.000000,0.000000,0.000000,0.294200,0.771615,1.065816 -2457386.29166667,0.000000,0.000000,0.000000,0.299113,0.768499,1.067612 -2457386.33333333,0.000000,0.000000,0.000000,0.303466,0.765396,1.068862 -2457386.37500000,0.000000,0.000000,0.000000,0.307289,0.762307,1.069596 -2457386.41666667,0.000000,0.000000,0.000000,0.310609,0.759231,1.069840 -2457386.45833333,0.000000,0.000000,0.000000,0.313454,0.756168,1.069622 -2457386.50000000,0.000000,0.000000,0.000000,0.315848,0.753118,1.068966 -2457386.54166667,0.000000,0.000000,0.000000,0.317817,0.750082,1.067898 -2457386.58333333,0.000000,0.000000,0.000000,0.319383,0.747058,1.066441 -2457386.62500000,0.000000,0.000000,0.000000,0.320569,0.744048,1.064617 -2457386.66666667,0.000000,0.000000,0.000000,0.321397,0.741051,1.062447 -2457386.70833333,0.000000,0.000000,0.000000,0.321886,0.738066,1.059952 -2457386.75000000,0.000000,0.000000,0.000000,0.322057,0.735095,1.057151 -2457386.79166667,0.000000,0.000000,0.000000,0.321928,0.732136,1.054064 -2457386.83333333,0.000000,0.000000,0.000000,0.321517,0.729190,1.050706 -2457386.87500000,0.000000,0.000000,0.000000,0.320840,0.726256,1.047096 -2457386.91666667,0.700000,0.006524,0.000391,0.319915,0.723422,1.043729 -2457386.95833333,1.400000,0.025855,0.004878,0.318757,0.720685,1.044321 -2457387.00000000,0.000000,0.000000,0.015012,0.317985,0.717788,1.050786 -2457387.04166667,0.000000,0.000000,0.008022,0.318170,0.714904,1.041096 -2457387.08333333,0.000000,0.000000,0.003298,0.318059,0.712032,1.033389 -2457387.12500000,0.100000,0.000134,0.000784,0.317667,0.709184,1.027635 -2457387.16666667,8.400001,0.852084,0.051193,0.317013,0.707280,1.075487 -2457387.20833333,34.799999,11.089222,1.099954,0.316199,0.707405,2.123558 -2457387.25000000,7.300000,0.652195,5.933234,0.321731,0.705396,6.960361 -2457387.29166667,4.700000,0.279236,3.556610,0.347321,0.703117,4.607048 -2457387.33333333,3.700000,0.175268,1.691810,0.376883,0.700736,2.769429 -2457387.37500000,4.200000,0.224403,0.591977,0.408213,0.698421,1.698611 -2457387.41666667,8.300000,0.832923,0.266570,0.440398,0.696552,1.403520 -2457387.45833333,3.300000,0.140139,0.525441,0.473751,0.694153,1.693345 -2457387.50000000,0.000000,0.000000,0.336875,0.511217,0.691370,1.539462 -2457387.54166667,0.000000,0.000000,0.145922,0.548764,0.688598,1.383284 -2457387.58333333,0.000000,0.000000,0.041804,0.583599,0.685839,1.311242 -2457387.62500000,0.000000,0.000000,0.004204,0.615849,0.683091,1.303144 -2457387.66666667,0.000000,0.000000,0.000000,0.645640,0.680354,1.325994 -2457387.70833333,0.000000,0.000000,0.000000,0.673088,0.677630,1.350718 -2457387.75000000,0.000000,0.000000,0.000000,0.698309,0.674917,1.373226 -2457387.79166667,0.000000,0.000000,0.000000,0.721410,0.672216,1.393626 -2457387.83333333,0.000000,0.000000,0.000000,0.742497,0.669526,1.412023 -2457387.87500000,0.000000,0.000000,0.000000,0.761670,0.666847,1.428517 -2457387.91666667,0.000000,0.000000,0.000000,0.779024,0.664180,1.443204 -2457387.95833333,0.000000,0.000000,0.000000,0.794651,0.661524,1.456175 -2457388.00000000,0.000000,0.000000,0.000000,0.808639,0.658880,1.467519 -2457388.04166667,0.000000,0.000000,0.000000,0.821073,0.656246,1.477320 -2457388.08333333,0.000000,0.000000,0.000000,0.832033,0.653624,1.485658 -2457388.12500000,0.000000,0.000000,0.000000,0.841597,0.651013,1.492610 -2457388.16666667,0.000000,0.000000,0.000000,0.849838,0.648413,1.498251 -2457388.20833333,0.000000,0.000000,0.000000,0.856827,0.645825,1.502652 -2457388.25000000,0.000000,0.000000,0.000000,0.862632,0.643247,1.505878 -2457388.29166667,0.000000,0.000000,0.000000,0.867317,0.640680,1.507997 -2457388.33333333,0.000000,0.000000,0.000000,0.870945,0.638123,1.509068 -2457388.37500000,0.000000,0.000000,0.000000,0.873575,0.635578,1.509153 -2457388.41666667,0.000000,0.000000,0.000000,0.875263,0.633044,1.508306 -2457388.45833333,0.000000,0.000000,0.000000,0.876063,0.630520,1.506583 -2457388.50000000,0.000000,0.000000,0.000000,0.876028,0.628006,1.504034 -2457388.54166667,0.000000,0.000000,0.000000,0.875206,0.625504,1.500710 -2457388.58333333,0.000000,0.000000,0.000000,0.873646,0.623012,1.496657 -2457388.62500000,0.000000,0.000000,0.000000,0.871391,0.620530,1.491921 -2457388.66666667,0.000000,0.000000,0.000000,0.868485,0.618059,1.486544 -2457388.70833333,0.000000,0.000000,0.000000,0.864969,0.615598,1.480568 -2457388.75000000,0.000000,0.000000,0.000000,0.860883,0.613148,1.474031 -2457388.79166667,0.000000,0.000000,0.000000,0.856263,0.610708,1.466971 -2457388.83333333,0.000000,0.000000,0.000000,0.851146,0.608278,1.459424 -2457388.87500000,0.000000,0.000000,0.000000,0.845565,0.605859,1.451424 -2457388.91666667,0.000000,0.000000,0.000000,0.839553,0.603449,1.443003 -2457388.95833333,0.000000,0.000000,0.000000,0.833141,0.601050,1.434191 -2457389.00000000,0.000000,0.000000,0.000000,0.826358,0.598661,1.425019 -2457389.04166667,0.000000,0.000000,0.000000,0.819232,0.596282,1.415514 -2457389.08333333,0.000000,0.000000,0.000000,0.811789,0.593913,1.405702 -2457389.12500000,0.000000,0.000000,0.000000,0.804056,0.591554,1.395610 -2457389.16666667,0.000000,0.000000,0.000000,0.796057,0.589204,1.385261 -2457389.20833333,0.000000,0.000000,0.000000,0.787814,0.586865,1.374678 -2457389.25000000,0.000000,0.000000,0.000000,0.779348,0.584535,1.363884 -2457389.29166667,0.000000,0.000000,0.000000,0.770682,0.582215,1.352898 -2457389.33333333,0.000000,0.000000,0.000000,0.761834,0.579905,1.341740 -2457389.37500000,0.000000,0.000000,0.000000,0.752824,0.577605,1.330429 -2457389.41666667,0.000000,0.000000,0.000000,0.743669,0.575314,1.318982 -2457389.45833333,0.000000,0.000000,0.000000,0.734385,0.573033,1.307417 -2457389.50000000,0.000000,0.000000,0.000000,0.724989,0.570761,1.295750 -2457389.54166667,0.000000,0.000000,0.000000,0.715496,0.568498,1.283994 -2457389.58333333,0.000000,0.000000,0.000000,0.705920,0.566246,1.272165 -2457389.62500000,0.000000,0.000000,0.000000,0.696274,0.564002,1.260276 -2457389.66666667,0.000000,0.000000,0.000000,0.686572,0.561768,1.248340 -2457389.70833333,0.000000,0.000000,0.000000,0.676826,0.559543,1.236369 -2457389.75000000,0.000000,0.000000,0.000000,0.667046,0.557328,1.224374 -2457389.79166667,0.000000,0.000000,0.000000,0.657244,0.555122,1.212366 -2457389.83333333,0.000000,0.000000,0.000000,0.647430,0.552925,1.200354 -2457389.87500000,0.000000,0.000000,0.000000,0.637613,0.550737,1.188350 -2457389.91666667,0.000000,0.000000,0.000000,0.627802,0.548558,1.176360 -2457389.95833333,0.000000,0.000000,0.000000,0.618006,0.546388,1.164395 -2457390.00000000,0.000000,0.000000,0.000000,0.608233,0.544228,1.152461 -2457390.04166667,0.000000,0.000000,0.000000,0.598490,0.542076,1.140566 -2457390.08333333,0.000000,0.000000,0.000000,0.588785,0.539933,1.128718 -2457390.12500000,0.000000,0.000000,0.000000,0.579123,0.537799,1.116922 -2457390.16666667,0.000000,0.000000,0.000000,0.569510,0.535674,1.105185 -2457390.20833333,0.000000,0.000000,0.000000,0.559954,0.533558,1.093512 -2457390.25000000,0.000000,0.000000,0.000000,0.550458,0.531451,1.081909 -2457390.29166667,0.000000,0.000000,0.000000,0.541028,0.529352,1.070380 -2457390.33333333,0.000000,0.000000,0.000000,0.531668,0.527262,1.058930 -2457390.37500000,0.000000,0.000000,0.000000,0.522383,0.525181,1.047564 -2457390.41666667,0.000000,0.000000,0.000000,0.513177,0.523108,1.036285 -2457390.45833333,0.000000,0.000000,0.000000,0.504052,0.521044,1.025097 -2457390.50000000,0.000000,0.000000,0.000000,0.495014,0.518989,1.014002 -2457390.54166667,0.000000,0.000000,0.000000,0.486064,0.516942,1.003006 -2457390.58333333,0.000000,0.000000,0.000000,0.477205,0.514903,0.992109 -2457390.62500000,0.000000,0.000000,0.000000,0.468441,0.512873,0.981314 -2457390.66666667,0.000000,0.000000,0.000000,0.459773,0.510851,0.970624 -2457390.70833333,0.000000,0.000000,0.000000,0.451204,0.508838,0.960041 -2457390.75000000,0.000000,0.000000,0.000000,0.442734,0.506833,0.949567 -2457390.79166667,0.000000,0.000000,0.000000,0.434367,0.504836,0.939203 -2457390.83333333,0.000000,0.000000,0.000000,0.426104,0.502848,0.928951 -2457390.87500000,0.000000,0.000000,0.000000,0.417945,0.500867,0.918812 -2457390.91666667,0.000000,0.000000,0.000000,0.409893,0.498895,0.908788 -2457390.95833333,0.000000,0.000000,0.000000,0.401947,0.496931,0.898878 -2457391.00000000,0.000000,0.000000,0.000000,0.394109,0.494975,0.889084 -2457391.04166667,0.000000,0.000000,0.000000,0.386380,0.493027,0.879407 -2457391.08333333,0.000000,0.000000,0.000000,0.378760,0.491088,0.869847 -2457391.12500000,0.000000,0.000000,0.000000,0.371249,0.489156,0.860405 -2457391.16666667,0.000000,0.000000,0.000000,0.363848,0.487232,0.851080 -2457391.20833333,0.000000,0.000000,0.000000,0.356557,0.485316,0.841873 -2457391.25000000,0.000000,0.000000,0.000000,0.349376,0.483408,0.832784 -2457391.29166667,0.000000,0.000000,0.000000,0.342305,0.481507,0.823812 -2457391.33333333,0.000000,0.000000,0.000000,0.335344,0.479615,0.814959 -2457391.37500000,0.000000,0.000000,0.000000,0.328492,0.477730,0.806222 -2457391.41666667,0.000000,0.000000,0.000000,0.321750,0.475854,0.797603 -2457391.45833333,0.000000,0.000000,0.000000,0.315116,0.473984,0.789100 -2457391.50000000,0.000000,0.000000,0.000000,0.308591,0.472123,0.780714 -2457391.54166667,0.000000,0.000000,0.000000,0.302174,0.470269,0.772443 -2457391.58333333,0.000000,0.000000,0.000000,0.295864,0.468423,0.764287 -2457391.62500000,0.000000,0.000000,0.000000,0.289661,0.466584,0.756245 -2457391.66666667,0.000000,0.000000,0.000000,0.283564,0.464753,0.748316 -2457391.70833333,0.000000,0.000000,0.000000,0.277571,0.462929,0.740501 -2457391.75000000,0.000000,0.000000,0.000000,0.271684,0.461113,0.732797 -2457391.79166667,0.000000,0.000000,0.000000,0.265899,0.459304,0.725203 -2457391.83333333,0.000000,0.000000,0.000000,0.260217,0.457503,0.717720 -2457391.87500000,0.000000,0.000000,0.000000,0.254637,0.455709,0.710346 -2457391.91666667,0.000000,0.000000,0.000000,0.249158,0.453922,0.703080 -2457391.95833333,0.000000,0.000000,0.000000,0.243778,0.452143,0.695921 -2457392.00000000,0.000000,0.000000,0.000000,0.238496,0.450371,0.688867 -2457392.04166667,0.000000,0.000000,0.000000,0.233312,0.448606,0.681918 -2457392.08333333,0.000000,0.000000,0.000000,0.228225,0.446848,0.675073 -2457392.12500000,0.000000,0.000000,0.000000,0.223232,0.445098,0.668330 -2457392.16666667,0.000000,0.000000,0.000000,0.218334,0.443355,0.661689 -2457392.20833333,0.000000,0.000000,0.000000,0.213529,0.441619,0.655148 -2457392.25000000,0.000000,0.000000,0.000000,0.208816,0.439890,0.648705 -2457392.29166667,0.000000,0.000000,0.000000,0.204193,0.438168,0.642360 -2457392.33333333,0.000000,0.000000,0.000000,0.199660,0.436453,0.636112 -2457392.37500000,0.000000,0.000000,0.000000,0.195215,0.434745,0.629959 -2457392.41666667,0.000000,0.000000,0.000000,0.190856,0.433044,0.623900 -2457392.45833333,0.000000,0.000000,0.000000,0.186584,0.431349,0.617934 -2457392.50000000,0.000000,0.000000,0.000000,0.182396,0.429662,0.612059 -2457392.54166667,0.000000,0.000000,0.000000,0.178292,0.427982,0.606274 -2457392.58333333,0.000000,0.000000,0.000000,0.174270,0.426308,0.600578 -2457392.62500000,0.000000,0.000000,0.000000,0.170328,0.424642,0.594970 -2457392.66666667,0.000000,0.000000,0.000000,0.166467,0.422982,0.589449 -2457392.70833333,0.000000,0.000000,0.000000,0.162683,0.421329,0.584012 -2457392.75000000,0.000000,0.000000,0.000000,0.158977,0.419682,0.578659 -2457392.79166667,0.000000,0.000000,0.000000,0.155347,0.418043,0.573389 -2457392.83333333,0.000000,0.000000,0.000000,0.151791,0.416410,0.568201 -2457392.87500000,0.000000,0.000000,0.000000,0.148309,0.414783,0.563093 -2457392.91666667,0.000000,0.000000,0.000000,0.144900,0.413163,0.558063 -2457392.95833333,0.000000,0.000000,0.000000,0.141561,0.411550,0.553111 -2457393.00000000,0.000000,0.000000,0.000000,0.138293,0.409943,0.548236 -2457393.04166667,0.000000,0.000000,0.000000,0.135093,0.408343,0.543436 -2457393.08333333,0.000000,0.000000,0.000000,0.131960,0.406749,0.538710 -2457393.12500000,0.000000,0.000000,0.000000,0.128894,0.405162,0.534056 -2457393.16666667,0.000000,0.000000,0.000000,0.125893,0.403581,0.529475 -2457393.20833333,0.000000,0.000000,0.000000,0.122957,0.402007,0.524964 -2457393.25000000,0.000000,0.000000,0.000000,0.120083,0.400439,0.520522 -2457393.29166667,0.000000,0.000000,0.000000,0.117271,0.398877,0.516148 -2457393.33333333,0.000000,0.000000,0.000000,0.114519,0.397322,0.511841 -2457393.37500000,0.000000,0.000000,0.000000,0.111827,0.395773,0.507600 -2457393.41666667,0.000000,0.000000,0.000000,0.109194,0.394230,0.503424 -2457393.45833333,0.000000,0.000000,0.000000,0.106618,0.392694,0.499311 -2457393.50000000,0.000000,0.000000,0.000000,0.104098,0.391163,0.495261 -2457393.54166667,0.000000,0.000000,0.000000,0.101633,0.389639,0.491272 -2457393.58333333,0.000000,0.000000,0.000000,0.099223,0.388121,0.487344 -2457393.62500000,0.000000,0.000000,0.000000,0.096866,0.386609,0.483475 -2457393.66666667,0.000000,0.000000,0.000000,0.094560,0.385104,0.479664 -2457393.70833333,0.000000,0.000000,0.000000,0.092306,0.383604,0.475910 -2457393.75000000,0.000000,0.000000,0.000000,0.090103,0.382110,0.472213 -2457393.79166667,0.000000,0.000000,0.000000,0.087948,0.380623,0.468571 -2457393.83333333,0.000000,0.000000,0.000000,0.085841,0.379141,0.464983 -2457393.87500000,0.000000,0.000000,0.000000,0.083782,0.377666,0.461448 -2457393.91666667,0.000000,0.000000,0.000000,0.081769,0.376196,0.457965 -2457393.95833333,0.000000,0.000000,0.000000,0.079802,0.374732,0.454534 -2457394.00000000,0.000000,0.000000,0.000000,0.077878,0.373275,0.451153 -2457394.04166667,0.000000,0.000000,0.000000,0.075999,0.371823,0.447822 -2457394.08333333,0.000000,0.000000,0.000000,0.074162,0.370377,0.444539 -2457394.12500000,0.000000,0.000000,0.000000,0.072367,0.368936,0.441303 -2457394.16666667,0.000000,0.000000,0.000000,0.070613,0.367502,0.438115 -2457394.20833333,0.000000,0.000000,0.000000,0.068899,0.366073,0.434972 -2457394.25000000,0.000000,0.000000,0.000000,0.067224,0.364650,0.431874 -2457394.29166667,0.000000,0.000000,0.000000,0.065588,0.363233,0.428821 -2457394.33333333,0.000000,0.000000,0.000000,0.063989,0.361822,0.425811 -2457394.37500000,0.000000,0.000000,0.000000,0.062428,0.360416,0.422843 -2457394.41666667,0.000000,0.000000,0.000000,0.060902,0.359016,0.419918 -2457394.45833333,0.000000,0.000000,0.000000,0.059412,0.357621,0.417033 -2457394.50000000,0.000000,0.000000,0.000000,0.057956,0.356232,0.414188 -2457394.54166667,0.000000,0.000000,0.000000,0.056534,0.354849,0.411383 -2457394.58333333,0.000000,0.000000,0.000000,0.055146,0.353471,0.408617 -2457394.62500000,0.000000,0.000000,0.000000,0.053790,0.352099,0.405888 -2457394.66666667,0.000000,0.000000,0.000000,0.052465,0.350732,0.403197 -2457394.70833333,0.000000,0.000000,0.000000,0.051172,0.349371,0.400542 -2457394.75000000,0.000000,0.000000,0.000000,0.049909,0.348015,0.397924 -2457394.79166667,0.000000,0.000000,0.000000,0.048676,0.346664,0.395340 -2457394.83333333,0.000000,0.000000,0.000000,0.047471,0.345319,0.392791 -2457394.87500000,0.000000,0.000000,0.000000,0.046296,0.343980,0.390276 -2457394.91666667,0.000000,0.000000,0.000000,0.045148,0.342646,0.387793 -2457394.95833333,0.000000,0.000000,0.000000,0.044027,0.341317,0.385344 -2457395.00000000,0.000000,0.000000,0.000000,0.042933,0.339993,0.382926 -2457395.04166667,0.000000,0.000000,0.000000,0.041865,0.338675,0.380540 -2457395.08333333,0.000000,0.000000,0.000000,0.040822,0.337362,0.378184 -2457395.12500000,0.000000,0.000000,0.000000,0.039804,0.336054,0.375859 -2457395.16666667,0.000000,0.000000,0.000000,0.038811,0.334752,0.373563 -2457395.20833333,0.000000,0.000000,0.000000,0.037841,0.333455,0.371296 -2457395.25000000,0.000000,0.000000,0.000000,0.036895,0.332162,0.369057 -2457395.29166667,0.000000,0.000000,0.000000,0.035971,0.330876,0.366847 -2457395.33333333,0.000000,0.000000,0.000000,0.035070,0.329594,0.364664 -2457395.37500000,0.000000,0.000000,0.000000,0.034190,0.328317,0.362507 -2457395.41666667,0.000000,0.000000,0.000000,0.033332,0.327046,0.360377 -2457395.45833333,0.000000,0.000000,0.000000,0.032494,0.325779,0.358273 -2457395.50000000,0.000000,0.000000,0.000000,0.031676,0.324518,0.356194 -2457395.54166667,0.000000,0.000000,0.000000,0.030879,0.323261,0.354140 -2457395.58333333,0.000000,0.000000,0.000000,0.030100,0.322010,0.352110 -2457395.62500000,0.000000,0.000000,0.000000,0.029341,0.320764,0.350105 -2457395.66666667,0.000000,0.000000,0.000000,0.028600,0.319522,0.348122 -2457395.70833333,0.000000,0.000000,0.000000,0.027877,0.318286,0.346163 -2457395.75000000,0.000000,0.000000,0.000000,0.027172,0.317054,0.344226 -2457395.79166667,0.000000,0.000000,0.000000,0.026484,0.315828,0.342311 -2457395.83333333,0.000000,0.000000,0.000000,0.025812,0.314606,0.340419 -2457395.87500000,0.000000,0.000000,0.000000,0.025158,0.313389,0.338547 -2457395.91666667,0.000000,0.000000,0.000000,0.024519,0.312177,0.336696 -2457395.95833333,0.000000,0.000000,0.000000,0.023896,0.310970,0.334866 -2457396.00000000,0.000000,0.000000,0.000000,0.023288,0.309768,0.333056 -2457396.04166667,0.000000,0.000000,0.000000,0.022695,0.308570,0.331266 -2457396.08333333,0.000000,0.000000,0.000000,0.022117,0.307378,0.329495 -2457396.12500000,0.000000,0.000000,0.000000,0.021553,0.306190,0.327743 -2457396.16666667,0.000000,0.000000,0.000000,0.021003,0.305006,0.326010 -2457396.20833333,0.000000,0.000000,0.000000,0.020467,0.303828,0.324295 -2457396.25000000,0.000000,0.000000,0.000000,0.019944,0.302654,0.322598 -2457396.29166667,0.000000,0.000000,0.000000,0.019434,0.301485,0.320918 -2457396.33333333,0.000000,0.000000,0.000000,0.018936,0.300320,0.319256 -2457396.37500000,0.000000,0.000000,0.000000,0.018451,0.299160,0.317611 -2457396.41666667,0.000000,0.000000,0.000000,0.017978,0.298005,0.315983 -2457396.45833333,0.000000,0.000000,0.000000,0.017517,0.296854,0.314371 -2457396.50000000,0.000000,0.000000,0.000000,0.017067,0.295708,0.312775 -2457396.54166667,0.000000,0.000000,0.000000,0.016628,0.294566,0.311194 -2457396.58333333,0.000000,0.000000,0.000000,0.016201,0.293429,0.309630 -2457396.62500000,0.000000,0.000000,0.000000,0.015784,0.292296,0.308080 -2457396.66666667,0.000000,0.000000,0.000000,0.015377,0.291168,0.306546 -2457396.70833333,0.000000,0.000000,0.000000,0.014981,0.290045,0.305026 -2457396.75000000,0.000000,0.000000,0.000000,0.014595,0.288926,0.303520 -2457396.79166667,0.000000,0.000000,0.000000,0.014218,0.287811,0.302029 -2457396.83333333,0.000000,0.000000,0.000000,0.013851,0.286701,0.300551 -2457396.87500000,0.000000,0.000000,0.000000,0.013493,0.285595,0.299088 -2457396.91666667,0.000000,0.000000,0.000000,0.013144,0.284493,0.297637 -2457396.95833333,0.000000,0.000000,0.000000,0.012804,0.283396,0.296200 -2457397.00000000,0.000000,0.000000,0.000000,0.012472,0.282303,0.294776 -2457397.04166667,0.000000,0.000000,0.000000,0.012149,0.281215,0.293364 -2457397.08333333,0.000000,0.000000,0.000000,0.011834,0.280131,0.291965 -2457397.12500000,0.000000,0.000000,0.000000,0.011527,0.279051,0.290578 -2457397.16666667,0.000000,0.000000,0.000000,0.011228,0.277975,0.289203 -2457397.20833333,0.000000,0.000000,0.000000,0.010936,0.276904,0.287840 -2457397.25000000,0.000000,0.000000,0.000000,0.010652,0.275837,0.286488 -2457397.29166667,0.000000,0.000000,0.000000,0.010375,0.274774,0.285148 -2457397.33333333,0.000000,0.000000,0.000000,0.010104,0.273715,0.283820 -2457397.37500000,0.000000,0.000000,0.000000,0.009841,0.272661,0.282502 -2457397.41666667,0.000000,0.000000,0.000000,0.009585,0.271610,0.281195 -2457397.45833333,0.000000,0.000000,0.000000,0.009335,0.270564,0.279899 -2457397.50000000,0.000000,0.000000,0.000000,0.009091,0.269522,0.278614 -2457397.54166667,0.000000,0.000000,0.000000,0.008854,0.268484,0.277338 -2457397.58333333,0.000000,0.000000,0.000000,0.008623,0.267451,0.276073 -2457397.62500000,0.000000,0.000000,0.000000,0.008397,0.266421,0.274818 -2457397.66666667,0.000000,0.000000,0.000000,0.008178,0.265395,0.273573 -2457397.70833333,0.000000,0.000000,0.000000,0.007964,0.264374,0.272337 -2457397.75000000,0.000000,0.000000,0.000000,0.007755,0.263356,0.271111 -2457397.79166667,0.000000,0.000000,0.000000,0.007552,0.262343,0.269894 -2457397.83333333,0.000000,0.000000,0.000000,0.007354,0.261333,0.268687 -2457397.87500000,0.000000,0.000000,0.000000,0.007161,0.260327,0.267488 -2457397.91666667,0.000000,0.000000,0.000000,0.006973,0.259326,0.266299 -2457397.95833333,0.000000,0.000000,0.000000,0.006790,0.258328,0.265118 -2457398.00000000,0.000000,0.000000,0.000000,0.006611,0.257334,0.263946 -2457398.04166667,0.000000,0.000000,0.000000,0.006438,0.256345,0.262782 -2457398.08333333,0.000000,0.000000,0.000000,0.006268,0.255359,0.261627 -2457398.12500000,0.000000,0.000000,0.000000,0.006103,0.254377,0.260480 -2457398.16666667,0.000000,0.000000,0.000000,0.005943,0.253399,0.259341 -2457398.20833333,0.000000,0.000000,0.000000,0.005786,0.252424,0.258210 -2457398.25000000,0.000000,0.000000,0.000000,0.005634,0.251454,0.257088 -2457398.29166667,0.000000,0.000000,0.000000,0.005485,0.250487,0.255972 -2457398.33333333,0.000000,0.000000,0.000000,0.005340,0.249525,0.254865 -2457398.37500000,0.000000,0.000000,0.000000,0.005199,0.248566,0.253765 -2457398.41666667,0.000000,0.000000,0.000000,0.005062,0.247610,0.252672 -2457398.45833333,0.000000,0.000000,0.000000,0.004928,0.246659,0.251587 -2457398.50000000,0.000000,0.000000,0.000000,0.004798,0.245711,0.250509 -2457398.54166667,0.000000,0.000000,0.000000,0.004671,0.244767,0.249438 -2457398.58333333,0.000000,0.000000,0.000000,0.004547,0.243827,0.248374 -2457398.62500000,0.000000,0.000000,0.000000,0.004427,0.242890,0.247317 -2457398.66666667,0.000000,0.000000,0.000000,0.004310,0.241957,0.246267 -2457398.70833333,0.000000,0.000000,0.000000,0.004195,0.241028,0.245223 -2457398.75000000,0.000000,0.000000,0.000000,0.004084,0.240102,0.244187 -2457398.79166667,0.000000,0.000000,0.000000,0.003976,0.239180,0.243156 -2457398.83333333,0.000000,0.000000,0.000000,0.003870,0.238262,0.242132 -2457398.87500000,0.000000,0.000000,0.000000,0.003768,0.237347,0.241115 -2457398.91666667,0.000000,0.000000,0.000000,0.003668,0.236436,0.240104 -2457398.95833333,0.000000,0.000000,0.000000,0.003570,0.235529,0.239099 -2457399.00000000,0.000000,0.000000,0.000000,0.003475,0.234625,0.238100 -2457399.04166667,0.000000,0.000000,0.000000,0.003383,0.233724,0.237107 -2457399.08333333,0.000000,0.000000,0.000000,0.003293,0.232827,0.236120 diff --git a/test/models/tshirt/include/TshirtTest.cpp b/test/models/tshirt/include/TshirtTest.cpp deleted file mode 100644 index ddac01238b..0000000000 --- a/test/models/tshirt/include/TshirtTest.cpp +++ /dev/null @@ -1,63 +0,0 @@ -#include "gtest/gtest.h" -#include "tshirt/include/Tshirt.h" -#include "tshirt/include/tshirt_params.h" -#include "GIUH.hpp" - -class TshirtModelTest : public ::testing::Test { - -protected: - - TshirtModelTest() { - - } - - ~TshirtModelTest() override { - - } - - void SetUp() override; - - void TearDown() override; - - void setupArbitraryExampleCase(); - -}; - -void TshirtModelTest::SetUp() { - -} - -void TshirtModelTest::TearDown() { - -} - -// Simple test to make sure the run function executes and that the inherent mass-balance check returned by run is good. -TEST_F(TshirtModelTest, TestRun0) { - - double et_storage = 0.0; - - tshirt::tshirt_params params{1000.0, 1.0, 10.0, 0.1, 0.01, 3, 1.0, 1.0, 1.0, 1.0, 8, 1.0, 1.0, 100.0}; - double storage = 1.0; - double input_flux = 10.0; - - // Testing with implied 0.0's state - tshirt::tshirt_state state_1(1.0, 1.0); - tshirt::tshirt_model model_1(params, make_shared(state_1)); - shared_ptr et_params_1 = make_shared(pdm03_struct()); - int model_1_result = model_1.run(86400.0, input_flux, et_params_1); - - // Testing with explicit state of correct size - tshirt::tshirt_state state_2(1.0, 1.0, - vector { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0 }); - tshirt::tshirt_model model_2(params, make_shared(state_2)); - shared_ptr et_params_2 = make_shared(pdm03_struct()); - int model_2_result = model_2.run(86400.0, input_flux, et_params_2); - - // TODO: figure out how to test for bogus/mismatched nash_n and state nash vector size (without silent error) - - // Should return 0 if mass balance check was good - EXPECT_EQ(model_1_result, 0); - EXPECT_EQ(model_2_result, 0); - -} - diff --git a/test/realizations/Formulation_Manager_Test.cpp b/test/realizations/Formulation_Manager_Test.cpp index a2f6609943..114ca246d6 100644 --- a/test/realizations/Formulation_Manager_Test.cpp +++ b/test/realizations/Formulation_Manager_Test.cpp @@ -596,7 +596,6 @@ TEST_F(Formulation_Manager_Test, basic_run_3) { std::vector actual_results(expected_results.size()); for (int i = 0; i < expected_results.size(); i++) { - // Remember that for the Tshirt_C_Realization, the timestep sizes are implicit actual_results[i] = manager.get_formulation("cat-67")->get_response(i, 3600); } diff --git a/test/realizations/catchments/Tshirt_C_Realization_Test.cpp b/test/realizations/catchments/Tshirt_C_Realization_Test.cpp deleted file mode 100644 index 81e0be04fd..0000000000 --- a/test/realizations/catchments/Tshirt_C_Realization_Test.cpp +++ /dev/null @@ -1,1079 +0,0 @@ -#include "gtest/gtest.h" -#include "Tshirt_C_Realization.hpp" -#include "tshirt/include/tshirt_params.h" -#include "GIUH.hpp" -#include -#include -#include -#include -#include -#include "FileChecker.h" - -class Tshirt_C_Realization_Test : public ::testing::Test { - -protected: - - typedef boost::tokenizer> Tokenizer; - - Tshirt_C_Realization_Test() { - error_bound_percent = 0.01; - error_upper_bound_min = 0.001; - upper_bound_factor = 1.0; - lower_bound_factor = 1.0; - - is_giuh_ordinate_examples = false; - is_forcing_params_examples = false; - } - - ~Tshirt_C_Realization_Test() override { - - } - - void SetUp() override; - - void TearDown() override; - - void open_standalone_c_impl_data_stream(); - - void setup_standalone_c_impl_example_case(); - - /** The name of the data file for data generated from the standalone C-based implementation. */ - std::string standalone_output_data_file; - /** The input stream of the data file for data generated from the standalone C-based implementation. */ - std::ifstream standalone_data_ingest_stream; - /** Whether there was an input stream opened or not for the current test instance. */ - bool is_standalone_data_ingest_stream; - std::shared_ptr c_impl_ex_tshirt_params; - - // TODO: may need to convert these - //std::shared_ptr c_impl_ex_initial_state; - //std::shared_ptr c_impl_ex_pdm_et_data; - - // Don't think this is needed, as it is implied now. - //double c_impl_ex_timestep_size_s; - - // TODO: probably don't need this - /* - ** - * The fluxes from Fred's standalone code are all in units of meters per time step. Multiply them by the - * "c_impl_ex_catchment_area_km2" variable to convert them into cubic meters per time step. - */ - //double c_impl_ex_catchment_area_km2; - - double error_bound_percent; - double error_upper_bound_min; - - double upper_bound_factor; - double lower_bound_factor; - - std::vector> giuh_ordinate_examples; - bool is_giuh_ordinate_examples; - std::vector forcing_params_examples; - bool is_forcing_params_examples; - - void init_forcing_params_examples(); - - void init_giuh_ordinate_examples(); - - -}; - -void Tshirt_C_Realization_Test::SetUp() { - - // If changed in debugging, modify here. - error_bound_percent = 0.01; - error_upper_bound_min = 0.001; - upper_bound_factor = 1.0 + error_bound_percent; - lower_bound_factor = 1.0 - error_bound_percent; - - // By default, this is false - is_standalone_data_ingest_stream = false; - - std::vector path_options = { - "test/data/model/tshirt/", - "./test/data/model/tshirt/", - "../test/data/model/tshirt/", - "../../test/data/model/tshirt/", - }; - - std::vector name_options = {"expected_results.csv"}; - - // Build vector of names by building combinations of the path and basename options - std::vector data_file_names(path_options.size() * name_options.size()); - - // Build so that all path names are tried for given basename before trying a different basename option - for (auto & name_option : name_options) { - for (auto & path_option : path_options) { - std::string string_combo = path_option + name_option; - data_file_names.push_back(string_combo); - } - } - - // Then go through in order and find the fist existing combination - standalone_output_data_file = utils::FileChecker::find_first_readable(data_file_names); - - //c_impl_ex_timestep_size_s = 60 * 60; // 60 seconds * 60 minutes => 1 hour size - - init_giuh_ordinate_examples(); - init_forcing_params_examples(); -} - -void Tshirt_C_Realization_Test::TearDown() { - if (is_standalone_data_ingest_stream && standalone_data_ingest_stream.is_open()) { - standalone_data_ingest_stream.close(); - } -} - -void Tshirt_C_Realization_Test::open_standalone_c_impl_data_stream() { - ASSERT_FALSE(standalone_output_data_file.empty()); - - standalone_data_ingest_stream.open(standalone_output_data_file.c_str()); - - ASSERT_TRUE(standalone_data_ingest_stream.is_open()); - - is_standalone_data_ingest_stream = true; -} - -void Tshirt_C_Realization_Test::setup_standalone_c_impl_example_case() { - - // Note using NGen value instead of Fred's value of 0.0 (which will cancel out things like Klf) - //double mult = 1.0; - // Valid range for this is 10-10000 - // TODO: considering tinkering with the value for this in multiple samples - // Also called 'lksatfac' - double mult = 1000.0; - - double maxsmc = 0.439; - //double wltsmc = 0.055; - double wltsmc = 0.066; - double satdk = 3.38e-06; - double satpsi = 0.355; - // High slop (0.0-1.0) will lead to almost no lateral flow (consider < 0.05) - // TODO: considering tinkering with the value for this in multiple samples - double slop = 1.0; - double bb = 4.05; - double alpha_fc = 0.33; - double expon = 6.0; - - // From Freds code; (approx) the average blue line drainage density for CONUS - double drainage_density_km_per_km2 = 3.5; - - double assumed_near_channel_water_table_slope = 0.01; - - // NGen has 0.0000672 - // see lines 326-329 in Fred's code - double Klf = 2.0 * assumed_near_channel_water_table_slope * mult * satdk * 2.0 * drainage_density_km_per_km2; - Klf *= TSHIRT_C_FIXED_TIMESTEP_SIZE_S; // Use this to convert Klf from m/s to m/timestep (m/h) - - // NGen has 1.08 - double Cgw = 0.01; // NGen has 1.08 - - // Note that NGen uses 8 - int nash_n = 2; - - // Note that NGen uses 0.1 - double Kn = 0.03; - - double max_gw_storage = 16.0; - - //Define tshirt params - //{maxsmc, wltsmc, satdk, satpsi, slope, b, multiplier, aplha_fx, klf, kn, nash_n, Cgw, expon, max_gw_storage} - c_impl_ex_tshirt_params = std::make_shared(tshirt::tshirt_params{ - maxsmc, //maxsmc FWRFH - wltsmc, //wltsmc - satdk, //satdk FWRFH - satpsi, //satpsi - slop, //slope - bb, //b bexp? FWRFH - mult, //multipier FIXMME what should this value be - alpha_fc, //aplha_fc - Klf, //Klf F - Kn, //Kn Kn 0.001-0.03 F - nash_n, //nash_n - Cgw, //Cgw C? FWRFH - expon, //expon FWRFH - max_gw_storage //max_gw_storage Sgwmax FWRFH - }); -} - -void Tshirt_C_Realization_Test::init_giuh_ordinate_examples() { - if (!is_giuh_ordinate_examples) { - std::vector giuh_ordinates_1(5); - giuh_ordinates_1[0] = 0.06; // note these sum to 1.0. If we have N ordinates, we need a queue sized N+1 to perform - giuh_ordinates_1[1] = 0.51; // the convolution. - giuh_ordinates_1[2] = 0.28; - giuh_ordinates_1[3] = 0.12; - giuh_ordinates_1[4] = 0.03; - giuh_ordinate_examples.push_back(giuh_ordinates_1); - - is_giuh_ordinate_examples = true; - } -} - -void Tshirt_C_Realization_Test::init_forcing_params_examples() { - if (!is_forcing_params_examples) { - std::vector forcing_path_opts = { - "test/data/forcing/cat-89_2015-12-01 00_00_00_2015-12-30 23_00_00.csv", - "./test/data/forcing/cat-89_2015-12-01 00_00_00_2015-12-30 23_00_00.csv", - "../test/data/forcing/cat-89_2015-12-01 00_00_00_2015-12-30 23_00_00.csv", - "../../test/data/forcing/cat-89_2015-12-01 00_00_00_2015-12-30 23_00_00.csv" - }; - std::string path = utils::FileChecker::find_first_readable(forcing_path_opts); - forcing_params_examples.push_back(forcing_params(path, "legacy", "2015-12-01 00:00:00", "2015-12-01 23:00:00")); - - is_forcing_params_examples = true; - } -} - -// Simple test to make sure the run function executes and that the inherent mass-balance check returned by run is good. -TEST_F(Tshirt_C_Realization_Test, TestRun0) { - int example_index = 0; - - tshirt::tshirt_params params{1000.0, 1.0, 10.0, 0.1, 0.01, 3, 1.0, 1.0, 1.0, 1.0, 2, 1.0, 1.0, 100.0}; - - std::vector nash_storage(params.nash_n); - for (int i = 0; i < params.nash_n; i++) { - nash_storage[i] = 0.0; - } - - std::vector giuh_ordinates = giuh_ordinate_examples[example_index]; - - realization::Tshirt_C_Realization tshirt_c_real( - forcing_params_examples[example_index], - utils::StreamHandler(), - 0.667, - 0.5, - true, - "wat-88", - giuh_ordinates, - params, - nash_storage); - - int result = tshirt_c_real.run_formulation_for_timestep(0.0, 3600); - - // TODO: figure out how to test for bogus/mismatched nash_n and state nash vector size (without silent error) - - // Should return 0 if mass balance check was good - EXPECT_EQ(result, 0); -} - -/** Test function for getting values for time step in formatted output string, for first and last time step. */ -TEST_F(Tshirt_C_Realization_Test, TestGetOutputLineForTimestep1a) { - int example_index = 0; - - open_standalone_c_impl_data_stream(); - - setup_standalone_c_impl_example_case(); - - // init gw res as half full for test - double gw_storage_ratio = 0.5; - - // init soil reservoir as 2/3 full - double soil_storage_ratio = 0.667; - - std::vector nash_storage(c_impl_ex_tshirt_params->nash_n); - for (int i = 0; i < c_impl_ex_tshirt_params->nash_n; i++) { - nash_storage[i] = 0.0; - } - - std::vector giuh_ordinates = giuh_ordinate_examples[example_index]; - - realization::Tshirt_C_Realization tshirt_c_real( - forcing_params_examples[example_index], - utils::StreamHandler(), - soil_storage_ratio, - gw_storage_ratio, - true, - "wat-88", - giuh_ordinates, - *c_impl_ex_tshirt_params, - nash_storage); - - std::vector result_vector; - std::string line; - int timestep = 0; - - while (getline(standalone_data_ingest_stream, line)) { - Tokenizer tokenizer(line); - result_vector.assign(tokenizer.begin(), tokenizer.end()); - - double input_storage = std::stod(result_vector[1]) / 1000 / 3600; - - // Output the line essentially - //copy(result_vector.begin(), result_vector.end(), std::ostream_iterator(std::cout, "|")); - //std::cout << "\n"; - - tshirt_c_real.run_formulation_for_timestep(input_storage, 3600); - timestep++; - } - - std::string actual_first = tshirt_c_real.get_output_line_for_timestep(0); - std::string actual_last = tshirt_c_real.get_output_line_for_timestep(timestep - 1); - std::string outside_bounds = tshirt_c_real.get_output_line_for_timestep(timestep); - - EXPECT_EQ(actual_first, "0.010000,0.001523,0.000091,0.000000,0.191106,0.191197"); - ASSERT_EQ(actual_last, "0.000000,0.000000,0.000000,0.000003,0.000233,0.000236"); -} - -/** Test function for getting values for time step in formatted output string handles out-of-bounds case as expected. */ -TEST_F(Tshirt_C_Realization_Test, TestGetOutputLineForTimestep1b) { - int example_index = 0; - - open_standalone_c_impl_data_stream(); - - setup_standalone_c_impl_example_case(); - - // init gw res as half full for test - double gw_storage_ratio = 0.5; - - // init soil reservoir as 2/3 full - double soil_storage_ratio = 0.667; - - std::vector nash_storage(c_impl_ex_tshirt_params->nash_n); - for (int i = 0; i < c_impl_ex_tshirt_params->nash_n; i++) { - nash_storage[i] = 0.0; - } - - std::vector giuh_ordinates = giuh_ordinate_examples[example_index]; - - realization::Tshirt_C_Realization tshirt_c_real( - forcing_params_examples[example_index], - utils::StreamHandler(), - soil_storage_ratio, - gw_storage_ratio, - true, - "wat-88", - giuh_ordinates, - *c_impl_ex_tshirt_params, - nash_storage); - - std::vector result_vector; - std::string line; - - int timestep = 0; - - while (getline(standalone_data_ingest_stream, line)) { - Tokenizer tokenizer(line); - result_vector.assign(tokenizer.begin(), tokenizer.end()); - double input_storage = std::stod(result_vector[1]) / 1000; - tshirt_c_real.run_formulation_for_timestep(input_storage, 3600); - timestep++; - } - - std::string actual_last = tshirt_c_real.get_output_line_for_timestep(timestep - 1); - std::string outside_bounds = tshirt_c_real.get_output_line_for_timestep(timestep); - - EXPECT_EQ(actual_last, "0.000000,0.000000,0.000000,0.000049,0.000336,0.000385"); - ASSERT_EQ(outside_bounds, ""); -} - -/** - * Test function for getting values for time step in formatted output string, for first and last time step, when reading - * from forcing data. - */ -TEST_F(Tshirt_C_Realization_Test, TestGetOutputLineForTimestep2a) { - int example_index = 0; - - open_standalone_c_impl_data_stream(); - - setup_standalone_c_impl_example_case(); - - // init gw res as half full for test - double gw_storage_ratio = 0.5; - - // init soil reservoir as 2/3 full - double soil_storage_ratio = 0.667; - - std::vector nash_storage(c_impl_ex_tshirt_params->nash_n); - for (int i = 0; i < c_impl_ex_tshirt_params->nash_n; i++) { - nash_storage[i] = 0.0; - } - - std::vector giuh_ordinates = giuh_ordinate_examples[example_index]; - - realization::Tshirt_C_Realization tshirt_c_real( - forcing_params_examples[example_index], - utils::StreamHandler(), - soil_storage_ratio, - gw_storage_ratio, - true, - "wat-88", - giuh_ordinates, - *c_impl_ex_tshirt_params, - nash_storage); - - std::vector result_vector; - std::string line; - - for (int timestep = 0; timestep < 5; ++timestep) { - tshirt_c_real.get_response(timestep, 3600); - } - - std::string actual_first = tshirt_c_real.get_output_line_for_timestep(0); - std::string actual_last = tshirt_c_real.get_output_line_for_timestep(4); - - EXPECT_EQ(actual_first, "0.000000,0.000000,0.000000,0.000000,0.191086,0.191086"); - ASSERT_EQ(actual_last, "0.000000,0.000000,0.000000,0.000242,0.145356,0.145598"); -} - -/** - * Test function for getting values for time step in formatted output string handles out-of-bounds case as expected, - * when reading from forcing data. - */ -TEST_F(Tshirt_C_Realization_Test, TestGetOutputLineForTimestep2b) { - int example_index = 0; - - open_standalone_c_impl_data_stream(); - - setup_standalone_c_impl_example_case(); - - // init gw res as half full for test - double gw_storage_ratio = 0.5; - - // init soil reservoir as 2/3 full - double soil_storage_ratio = 0.667; - - std::vector nash_storage(c_impl_ex_tshirt_params->nash_n); - for (int i = 0; i < c_impl_ex_tshirt_params->nash_n; i++) { - nash_storage[i] = 0.0; - } - - std::vector giuh_ordinates = giuh_ordinate_examples[example_index]; - - realization::Tshirt_C_Realization tshirt_c_real( - forcing_params_examples[example_index], - utils::StreamHandler(), - soil_storage_ratio, - gw_storage_ratio, - true, - "wat-88", - giuh_ordinates, - *c_impl_ex_tshirt_params, - nash_storage); - - std::vector result_vector; - std::string line; - - for (int timestep = 0; timestep < 5; ++timestep) { - tshirt_c_real.get_response(timestep, 3600); - } - - std::string actual_last = tshirt_c_real.get_output_line_for_timestep(4); - std::string outside_bounds = tshirt_c_real.get_output_line_for_timestep(5); - - EXPECT_EQ(actual_last, "0.000000,0.000000,0.000000,0.000242,0.145356,0.145598"); - ASSERT_EQ(outside_bounds, ""); -} - -/** Test function for getting number of output variables for realization type. */ -TEST_F(Tshirt_C_Realization_Test, TestGetOutputItemCount1a) { - int example_index = 0; - - open_standalone_c_impl_data_stream(); - - setup_standalone_c_impl_example_case(); - - // init gw res as half full for test - double gw_storage_ratio = 0.5; - - // init soil reservoir as 2/3 full - double soil_storage_ratio = 0.667; - - std::vector nash_storage(c_impl_ex_tshirt_params->nash_n); - for (int i = 0; i < c_impl_ex_tshirt_params->nash_n; i++) { - nash_storage[i] = 0.0; - } - - std::vector giuh_ordinates = giuh_ordinate_examples[example_index]; - - realization::Tshirt_C_Realization tshirt_c_real( - forcing_params_examples[example_index], - utils::StreamHandler(), - soil_storage_ratio, - gw_storage_ratio, - true, - "wat-88", - giuh_ordinates, - *c_impl_ex_tshirt_params, - nash_storage); - - ASSERT_EQ(tshirt_c_real.get_output_item_count(), tshirt_c_real.get_output_var_names().size()); -} - -/** Test header fields are reasonable (same size collection) for output variables. */ -TEST_F(Tshirt_C_Realization_Test, TestGetOutputHeaderFields1a) { - int example_index = 0; - - open_standalone_c_impl_data_stream(); - - setup_standalone_c_impl_example_case(); - - // init gw res as half full for test - double gw_storage_ratio = 0.5; - - // init soil reservoir as 2/3 full - double soil_storage_ratio = 0.667; - - std::vector nash_storage(c_impl_ex_tshirt_params->nash_n); - for (int i = 0; i < c_impl_ex_tshirt_params->nash_n; i++) { - nash_storage[i] = 0.0; - } - - std::vector giuh_ordinates = giuh_ordinate_examples[example_index]; - - realization::Tshirt_C_Realization tshirt_c_real( - forcing_params_examples[example_index], - utils::StreamHandler(), - soil_storage_ratio, - gw_storage_ratio, - true, - "wat-88", - giuh_ordinates, - *c_impl_ex_tshirt_params, - nash_storage); - - ASSERT_EQ(tshirt_c_real.get_output_header_fields().size(), tshirt_c_real.get_output_var_names().size()); -} - -/** Test function for getting the output variable names for realization type. */ -TEST_F(Tshirt_C_Realization_Test, TestGetOutputVariableNames1a) { - int example_index = 0; - - open_standalone_c_impl_data_stream(); - - setup_standalone_c_impl_example_case(); - - // init gw res as half full for test - double gw_storage_ratio = 0.5; - - // init soil reservoir as 2/3 full - double soil_storage_ratio = 0.667; - - std::vector nash_storage(c_impl_ex_tshirt_params->nash_n); - for (int i = 0; i < c_impl_ex_tshirt_params->nash_n; i++) { - nash_storage[i] = 0.0; - } - - std::vector giuh_ordinates = giuh_ordinate_examples[example_index]; - - realization::Tshirt_C_Realization tshirt_c_real( - forcing_params_examples[example_index], - utils::StreamHandler(), - soil_storage_ratio, - gw_storage_ratio, - true, - "wat-88", - giuh_ordinates, - *c_impl_ex_tshirt_params, - nash_storage); - - std::vector var_names = { - OUT_VAR_RAINFALL, - OUT_VAR_SURFACE_RUNOFF, - OUT_VAR_GIUH_RUNOFF, - OUT_VAR_LATERAL_FLOW, - OUT_VAR_BASE_FLOW, - OUT_VAR_TOTAL_DISCHARGE - }; - - ASSERT_EQ(tshirt_c_real.get_output_var_names(), var_names); -} - -/** Test values getting function for surface runoff. */ -TEST_F(Tshirt_C_Realization_Test, TestGetValue1a) { - int example_index = 0; - - open_standalone_c_impl_data_stream(); - - setup_standalone_c_impl_example_case(); - - // init gw res as half full for test - double gw_storage_ratio = 0.5; - - // init soil reservoir as 2/3 full - double soil_storage_ratio = 0.667; - - std::vector nash_storage(c_impl_ex_tshirt_params->nash_n); - for (int i = 0; i < c_impl_ex_tshirt_params->nash_n; i++) { - nash_storage[i] = 0.0; - } - - std::vector giuh_ordinates = giuh_ordinate_examples[example_index]; - - realization::Tshirt_C_Realization tshirt_c_real( - forcing_params_examples[example_index], - utils::StreamHandler(), - soil_storage_ratio, - gw_storage_ratio, - true, - "wat-88", - giuh_ordinates, - *c_impl_ex_tshirt_params, - nash_storage); - - std::vector result_vector; - std::vector values_vector; - std::string line; - - while (getline(standalone_data_ingest_stream, line)) { - Tokenizer tokenizer(line); - result_vector.assign(tokenizer.begin(), tokenizer.end()); - double input_storage = std::stod(result_vector[1]); - tshirt_c_real.run_formulation_for_timestep(input_storage, 3600); - values_vector.emplace_back(tshirt_c_real.get_latest_flux_surface_runoff()); - } - - ASSERT_EQ(tshirt_c_real.get_value(OUT_VAR_SURFACE_RUNOFF), values_vector); -} - -/** Test values getting function for total discharge. */ -TEST_F(Tshirt_C_Realization_Test, TestGetValue1b) { - int example_index = 0; - - open_standalone_c_impl_data_stream(); - - setup_standalone_c_impl_example_case(); - - // init gw res as half full for test - double gw_storage_ratio = 0.5; - - // init soil reservoir as 2/3 full - double soil_storage_ratio = 0.667; - - std::vector nash_storage(c_impl_ex_tshirt_params->nash_n); - for (int i = 0; i < c_impl_ex_tshirt_params->nash_n; i++) { - nash_storage[i] = 0.0; - } - - std::vector giuh_ordinates = giuh_ordinate_examples[example_index]; - - realization::Tshirt_C_Realization tshirt_c_real( - forcing_params_examples[example_index], - utils::StreamHandler(), - soil_storage_ratio, - gw_storage_ratio, - true, - "wat-88", - giuh_ordinates, - *c_impl_ex_tshirt_params, - nash_storage); - - std::vector result_vector; - std::vector values_vector; - std::string line; - - while (getline(standalone_data_ingest_stream, line)) { - Tokenizer tokenizer(line); - result_vector.assign(tokenizer.begin(), tokenizer.end()); - double input_storage = std::stod(result_vector[1]); - tshirt_c_real.run_formulation_for_timestep(input_storage, 3600); - values_vector.emplace_back(tshirt_c_real.get_latest_flux_total_discharge()); - } - - ASSERT_EQ(tshirt_c_real.get_value(OUT_VAR_TOTAL_DISCHARGE), values_vector); -} - -/** Test direct surface runoff flux calculations for first example, within bounds. */ -TEST_F(Tshirt_C_Realization_Test, TestSurfaceRunoffCalc1a) { - int example_index = 0; - - open_standalone_c_impl_data_stream(); - - setup_standalone_c_impl_example_case(); - - // init gw res as half full for test - double gw_storage_ratio = 0.5; - - // init soil reservoir as 2/3 full - double soil_storage_ratio = 0.667; - - std::vector nash_storage(c_impl_ex_tshirt_params->nash_n); - for (int i = 0; i < c_impl_ex_tshirt_params->nash_n; i++) { - nash_storage[i] = 0.0; - } - - std::vector giuh_ordinates = giuh_ordinate_examples[example_index]; - - realization::Tshirt_C_Realization tshirt_c_real( - forcing_params_examples[example_index], - utils::StreamHandler(), - soil_storage_ratio, - gw_storage_ratio, - true, - "wat-88", - giuh_ordinates, - *c_impl_ex_tshirt_params, - nash_storage); - - std::vector result_vector; - std::string line; - - while (getline(standalone_data_ingest_stream, line)) { - Tokenizer tokenizer(line); - result_vector.assign(tokenizer.begin(), tokenizer.end()); - - double input_storage = std::stod(result_vector[1]); - // The fluxes from Fred's code are all in units of meters per time step. Multiply them by the "c_impl_ex_catchment_area_km2" - // variable to convert them into cubic meters per time step. - //input_storage *= c_impl_ex_catchment_area_km2; - - // However, the data starts in mm/h, so first convert mm to m ... - input_storage /= 1000; - // ... then convert m / h to m / s - input_storage /= 3600; - - // runoff is index 2 - double expected = std::stod(result_vector[2]); - - // Also convert expected from mm / h to m / h - expected /= 1000; - - // Output the line essentially - copy(result_vector.begin(), result_vector.end(), std::ostream_iterator(std::cout, "|")); - std::cout << "\n"; - - tshirt_c_real.run_formulation_for_timestep(input_storage, 3600); - - double actual = tshirt_c_real.get_latest_flux_surface_runoff(); - - // Note that, for non-zero values, having to work within a reasonable upper and lower bounds to allow for - // precision and rounding errors both on the calculation and sample-data-recording side. - if (expected == 0.0) { - EXPECT_EQ(actual, expected); - } - else { - double upper_bound = expected * upper_bound_factor; - double lower_bound = expected * lower_bound_factor; - EXPECT_LE(actual, upper_bound); - EXPECT_GE(actual, lower_bound); - } - - } -} - -/** Test GIUH runoff calculations. */ -TEST_F(Tshirt_C_Realization_Test, TestGiuhRunoffCalc1a) { - int example_index = 0; - - open_standalone_c_impl_data_stream(); - - setup_standalone_c_impl_example_case(); - - // init gw res as half full for test - double gw_storage_ratio = 0.5; - - // init soil reservoir as 2/3 full - double soil_storage_ratio = 0.667; - - std::vector nash_storage(c_impl_ex_tshirt_params->nash_n); - for (int i = 0; i < c_impl_ex_tshirt_params->nash_n; i++) { - nash_storage[i] = 0.0; - } - - std::vector giuh_ordinates = giuh_ordinate_examples[example_index]; - - realization::Tshirt_C_Realization tshirt_c_real( - forcing_params_examples[example_index], - utils::StreamHandler(), - soil_storage_ratio, - gw_storage_ratio, - true, - "wat-88", - giuh_ordinates, - *c_impl_ex_tshirt_params, - nash_storage); - - std::vector result_vector; - std::string line; - - while (getline(standalone_data_ingest_stream, line)) { - Tokenizer tokenizer(line); - result_vector.assign(tokenizer.begin(), tokenizer.end()); - - double input_storage = std::stod(result_vector[1]); - // The fluxes from Fred's code are all in units of meters per time step. Multiply them by the "c_impl_ex_catchment_area_km2" - // variable to convert them into cubic meters per time step. - //input_storage *= c_impl_ex_catchment_area_km2; - - // However, the data starts in mm/h, so first convert mm to m ... - input_storage /= 1000; - // ... then convert m / h to m / s - input_storage /= 3600; - - // giuh runoff is index 3 - double expected = std::stod(result_vector[3]); - - // Also convert expected from mm / h to m / h - expected /= 1000; - - // Output the line essentially - copy(result_vector.begin(), result_vector.end(), std::ostream_iterator(std::cout, "|")); - std::cout << "\n"; - - tshirt_c_real.run_formulation_for_timestep(input_storage, 3600); - double actual = tshirt_c_real.get_latest_flux_giuh_runoff(); - - // Note that, for non-zero values, having to work within a reasonable upper and lower bounds to allow for - // precision and rounding errors both on the calculation and sample-data-recording side. - if (expected == 0.0) { - // Might be acceptable if actual value is less than 10^-9 - //EXPECT_LT(actual - expected, 0.000000001); - // ... but for now, do this - EXPECT_EQ(actual, expected); - } - else { - double upper_bound = expected * upper_bound_factor; - double lower_bound = expected * lower_bound_factor; - // TODO: apply this logic to the other tests also. - // Keep in mind that the sample data only has so much precision (for brevity/formatting purposes). As such, - // upper and lower bounds for this test also need to be limited by what is representable. - // This is the smallest value that can be written to test data with current formatting (then divided by - // 1000, since in millimeters) - double representable_bound_limit = 0.000001 / 1000; - double rep_upper_bound = expected + representable_bound_limit; - double rep_lower_bound = expected - representable_bound_limit; - upper_bound = max(upper_bound, rep_upper_bound); - lower_bound = min(lower_bound, rep_lower_bound); - - EXPECT_LE(actual, upper_bound); - EXPECT_GE(actual, lower_bound); - } - } -} - -/** Test soil lateral flow calculations. */ -TEST_F(Tshirt_C_Realization_Test, TestLateralFlowCalc1a) { - int example_index = 0; - - open_standalone_c_impl_data_stream(); - - setup_standalone_c_impl_example_case(); - - // init gw res as half full for test - double gw_storage_ratio = 0.5; - - // init soil reservoir as 2/3 full - double soil_storage_ratio = 0.667; - - std::vector nash_storage(c_impl_ex_tshirt_params->nash_n); - for (int i = 0; i < c_impl_ex_tshirt_params->nash_n; i++) { - nash_storage[i] = 0.0; - } - - std::vector giuh_ordinates = giuh_ordinate_examples[example_index]; - - realization::Tshirt_C_Realization tshirt_c_real( - forcing_params_examples[example_index], - utils::StreamHandler(), - soil_storage_ratio, - gw_storage_ratio, - true, - "wat-88", - giuh_ordinates, - *c_impl_ex_tshirt_params, - nash_storage); - - std::vector result_vector; - std::string line; - - while (getline(standalone_data_ingest_stream, line)) { - Tokenizer tokenizer(line); - result_vector.assign(tokenizer.begin(), tokenizer.end()); - - double input_storage = std::stod(result_vector[1]); - // The fluxes from Fred's code are all in units of meters per time step. Multiply them by the "c_impl_ex_catchment_area_km2" - // variable to convert them into cubic meters per time step. - //input_storage *= c_impl_ex_catchment_area_km2; - - // However, the data starts in mm/h, so first convert mm to m ... - input_storage /= 1000; - // ... then convert m / h to m / s - input_storage /= 3600; - - // lateral flow is index 4 - double expected = std::stod(result_vector[4]); - - // Also convert expected from mm / h to m / h - expected /= 1000; - - // Output the line essentially - copy(result_vector.begin(), result_vector.end(), std::ostream_iterator(std::cout, "|")); - std::cout << "\n"; - - tshirt_c_real.run_formulation_for_timestep(input_storage, 3600); - double actual = tshirt_c_real.get_latest_flux_lateral_flow(); - - // Note that, for non-zero values, having to work within a reasonable upper and lower bounds to allow for - // precision and rounding errors both on the calculation and sample-data-recording side. - if (expected == 0.0) { - // Might be acceptable if actual value is less than 10^-9 - //EXPECT_LT(actual - expected, 0.000000001); - // ... but for now, do this - EXPECT_EQ(actual, expected); - } - else { - double upper_bound = expected * upper_bound_factor; - double lower_bound = expected * lower_bound_factor; - EXPECT_LE(actual, upper_bound); - EXPECT_GE(actual, lower_bound); - } - } -} - -/** Test base flow calculations. */ -TEST_F(Tshirt_C_Realization_Test, TestBaseFlowCalc1a) { - int example_index = 0; - - open_standalone_c_impl_data_stream(); - - setup_standalone_c_impl_example_case(); - - // init gw res as half full for test - double gw_storage_ratio = 0.5; - - // init soil reservoir as 2/3 full - double soil_storage_ratio = 0.667; - - std::vector nash_storage(c_impl_ex_tshirt_params->nash_n); - for (int i = 0; i < c_impl_ex_tshirt_params->nash_n; i++) { - nash_storage[i] = 0.0; - } - - std::vector giuh_ordinates = giuh_ordinate_examples[example_index]; - - realization::Tshirt_C_Realization tshirt_c_real( - forcing_params_examples[example_index], - utils::StreamHandler(), - soil_storage_ratio, - gw_storage_ratio, - true, - "wat-88", - giuh_ordinates, - *c_impl_ex_tshirt_params, - nash_storage); - - std::vector result_vector; - std::string line; - - while (getline(standalone_data_ingest_stream, line)) { - Tokenizer tokenizer(line); - result_vector.assign(tokenizer.begin(), tokenizer.end()); - - double input_storage = std::stod(result_vector[1]); - // The fluxes from Fred's code are all in units of meters per time step. Multiply them by the "c_impl_ex_catchment_area_km2" - // variable to convert them into cubic meters per time step. - //input_storage *= c_impl_ex_catchment_area_km2; - - // However, the data starts in mm/h, so first convert mm to m ... - input_storage /= 1000; - // ... then convert m / h to m / s - input_storage /= 3600; - - // base flow is index 5 - double expected = std::stod(result_vector[5]); - - // Also convert expected from mm / h to m / h - expected /= 1000; - - // Output the line essentially - copy(result_vector.begin(), result_vector.end(), std::ostream_iterator(std::cout, "|")); - std::cout << "\n"; - - tshirt_c_real.run_formulation_for_timestep(input_storage, 3600); - double actual = tshirt_c_real.get_latest_flux_base_flow(); - - // Note that, for non-zero values, having to work within a reasonable upper and lower bounds to allow for - // precision and rounding errors both on the calculation and sample-data-recording side. - if (expected == 0.0) { - // Might be acceptable if actual value is less than 10^-9 - //EXPECT_LT(actual - expected, 0.000000001); - // ... but for now, do this - EXPECT_EQ(actual, expected); - } - else { - double upper_bound = expected * upper_bound_factor; - double lower_bound = expected * lower_bound_factor; - EXPECT_LE(actual, upper_bound); - EXPECT_GE(actual, lower_bound); - } - } -} - -/** Test total discharge output calculations. */ -TEST_F(Tshirt_C_Realization_Test, TestTotalDischargeOutputCalc1a) { - int example_index = 0; - - open_standalone_c_impl_data_stream(); - - setup_standalone_c_impl_example_case(); - - // init gw res as half full for test - double gw_storage_ratio = 0.5; - - // init soil reservoir as 2/3 full - double soil_storage_ratio = 0.667; - - std::vector nash_storage(c_impl_ex_tshirt_params->nash_n); - for (int i = 0; i < c_impl_ex_tshirt_params->nash_n; i++) { - nash_storage[i] = 0.0; - } - - std::vector giuh_ordinates = giuh_ordinate_examples[example_index]; - - realization::Tshirt_C_Realization tshirt_c_real( - forcing_params_examples[example_index], - utils::StreamHandler(), - soil_storage_ratio, - gw_storage_ratio, - true, - "wat-88", - giuh_ordinates, - *c_impl_ex_tshirt_params, - nash_storage); - - std::vector result_vector; - std::string line; - - while (getline(standalone_data_ingest_stream, line)) { - Tokenizer tokenizer(line); - result_vector.assign(tokenizer.begin(), tokenizer.end()); - - double input_storage = std::stod(result_vector[1]); - // The fluxes from Fred's code are all in units of meters per time step. Multiply them by the "c_impl_ex_catchment_area_km2" - // variable to convert them into cubic meters per time step. - //input_storage *= c_impl_ex_catchment_area_km2; - - // However, the data starts in mm/h, so first convert mm to m ... - input_storage /= 1000; - // ... then convert m / h to m / s - input_storage /= 3600; - - // output is index 5 - double expected = std::stod(result_vector[6]); - - // Also convert expected from mm / h to m / h - expected /= 1000; - - // Output the line essentially - copy(result_vector.begin(), result_vector.end(), std::ostream_iterator(std::cout, "|")); - std::cout << "\n"; - - tshirt_c_real.run_formulation_for_timestep(input_storage, 3600); - double actual = tshirt_c_real.get_latest_flux_total_discharge(); - - // Note that, for non-zero values, having to work within a reasonable upper and lower bounds to allow for - // precision and rounding errors both on the calculation and sample-data-recording side. - if (expected == 0.0) { - // Might be acceptable if actual value is less than 10^-9 - //EXPECT_LT(actual - expected, 0.000000001); - // ... but for now, do this - EXPECT_EQ(actual, expected); - } - else { - double upper_bound = expected * upper_bound_factor; - double lower_bound = expected * lower_bound_factor; - EXPECT_LE(actual, upper_bound); - EXPECT_GE(actual, lower_bound); - } - } -}