distribution/power_law.h
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /// @file | ||
| 2 | /// @brief Power-law size-frequency distributions and their segments. | ||
| 3 | #pragma once | ||
| 4 | |||
| 5 | #include "terraingen/distribution/size_frequency.h" | ||
| 6 | |||
| 7 | #include <cmath> | ||
| 8 | |||
| 9 | namespace terraingen { | ||
| 10 | |||
| 11 | /// A power law N(>=D) = `coefficient` * D^(-`exponent`) from `start` metres, alone or as | ||
| 12 | /// one segment of a `PiecewisePowerLawSfd`. | ||
| 13 | struct PowerLawSegment { | ||
| 14 | 8 | double start = 0.0; ///< First diameter, metres. | |
| 15 | 8 | double coefficient = 1.0; ///< N(>=1 m), per m^2 (or per m^2 per Gyr); positive. | |
| 16 | 8 | double exponent = 2.0; ///< Cumulative slope; positive. | |
| 17 | |||
| 18 | /// N(>=D) = coefficient * D^(-exponent), with D = `diameter` in metres. | ||
| 19 | 26014543 | double Cumulative(double diameter) const { | |
| 20 | 26014543 | return coefficient * std::pow(diameter, -exponent); | |
| 21 | } | ||
| 22 | |||
| 23 | /// -dN/dD = coefficient * exponent * D^(-exponent - 1), with D in metres. | ||
| 24 | 2285088 | double Incremental(double diameter) const { | |
| 25 | 2285088 | return coefficient * exponent * std::pow(diameter, -exponent - 1.0); | |
| 26 | } | ||
| 27 | |||
| 28 | /// The diameter in metres where N(>=D) is a positive `target`. | ||
| 29 | 1579674 | double Invert(double target) const { | |
| 30 | 1579674 | return std::pow(target / coefficient, -1.0 / exponent); | |
| 31 | } | ||
| 32 | |||
| 33 | 8 | bool operator==(const PowerLawSegment&) const = default; | |
| 34 | }; | ||
| 35 | |||
| 36 | /// N(>=D) = `coefficient` * D^(-`exponent`), with D in metres and N per m^2 (or per m^2 per | ||
| 37 | /// Gyr). | ||
| 38 | class PowerLawSfd final : public SizeFrequencyDistribution { | ||
| 39 | public: | ||
| 40 | /// The power law with the given coefficient and exponent. | ||
| 41 | /// @param coefficient N(>=1 m), per m^2 (or per m^2 per Gyr). | ||
| 42 | /// @param range The valid range, in metres. | ||
| 43 | /// @throws std::invalid_argument unless `coefficient` and `exponent` are positive and finite | ||
| 44 | /// and the range is valid. | ||
| 45 | PowerLawSfd(double coefficient, double exponent, SfdKind kind = SfdKind::Density, | ||
| 46 | DiameterRange range = {}); | ||
| 47 | |||
| 48 | /// N(>=1 m), per m^2 (or per m^2 per Gyr). | ||
| 49 | 2 | double Coefficient() const { | |
| 50 | 2 | return law_.coefficient; | |
| 51 | } | ||
| 52 | |||
| 53 | /// The cumulative slope. | ||
| 54 | 2 | double Exponent() const { | |
| 55 | 2 | return law_.exponent; | |
| 56 | } | ||
| 57 | |||
| 58 | protected: | ||
| 59 | double CumulativeAt(double diameter) const override; | ||
| 60 | double IncrementalAt(double diameter) const override; | ||
| 61 | double InvertCumulative(double target, double dMin, double dMax) const override; | ||
| 62 | |||
| 63 | private: | ||
| 64 | PowerLawSegment law_; | ||
| 65 | }; | ||
| 66 | |||
| 67 | /// Converts a power-law coefficient published for D in km and N per km^2 to D in metres and | ||
| 68 | /// N per m^2. | ||
| 69 | double PowerLawCoefficientFromKilometres(double coefficient, double exponent); | ||
| 70 | |||
| 71 | } // namespace terraingen | ||
| 72 |