layers/craters/profiles/complex.cpp

Line	Branch	Exec	Source
1			#include "terraingen/layers/craters/profiles/complex.h"
2
3			#include "terraingen/util/checks.h"
4
5			#include <algorithm>
6			#include <cmath>
7			#include <numbers>
8			#include <stdexcept>
9			#include <utility>
10
11			namespace terraingen {
12
13		99798	double KilometrePowerLaw::At(double diameter, bool volume) const {
14		99798	const double kilometres = coefficient * std::pow(diameter / 1000.0, exponent);
15	2/2 ✓ Branch 0 taken 12687 times. ✓ Branch 1 taken 87111 times.	99798	return volume ? kilometres * 1e9 : kilometres * 1000.0;
16			}
17
18		19	void ValidateComplexCraterShape(const ComplexCraterShape& shape) {
19	2/2 ✓ Branch 1 taken 1 time. ✓ Branch 2 taken 18 times.	19	if (!IsPositive(shape.transitionDiameter)) {
20	1/1 ✓ Branch 2 taken 1 time.	1	throw std::invalid_argument("The complex crater transition diameter must be positive");
21			}
22	2/2 ✓ Branch 1 taken 4 times. ✓ Branch 2 taken 14 times.	18	const KilometrePowerLaw peakWidth = shape.peakDiameter ? *shape.peakDiameter : shape.peakVolume;
23	2/2 ✓ Branch 2 taken 54 times. ✓ Branch 3 taken 16 times.	70	for (const KilometrePowerLaw& law : {shape.floorDiameter, shape.peakHeight, peakWidth}) {
24	5/6 ✓ Branch 1 taken 52 times. ✓ Branch 2 taken 2 times. ✗ Branch 4 not taken. ✓ Branch 5 taken 52 times. ✓ Branch 6 taken 2 times. ✓ Branch 7 taken 52 times.	54	if (!(IsPositive(law.coefficient) && std::isfinite(law.exponent))) {
25	1/1 ✓ Branch 2 taken 2 times.	2	throw std::invalid_argument("Complex crater power laws need a positive coefficient");
26			}
27			}
28	2/2 ✓ Branch 1 taken 1 time. ✓ Branch 2 taken 15 times.	16	if (!IsNonNegative(shape.peakFitMaxDiameter)) {
29	1/1 ✓ Branch 2 taken 1 time.	1	throw std::invalid_argument("The peak fit limit must not be negative");
30			}
31	4/6 ✓ Branch 1 taken 15 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 15 times. ✗ Branch 5 not taken. ✓ Branch 6 taken 1 time. ✓ Branch 7 taken 14 times.	30	if (!(IsPositive(shape.ejectaExponent) && std::isfinite(shape.reachRadii) &&
32	2/2 ✓ Branch 0 taken 1 time. ✓ Branch 1 taken 14 times.	15	shape.reachRadii > 1.0)) {
33			throw std::invalid_argument(
34	1/1 ✓ Branch 2 taken 1 time.	1	"Complex craters need a positive ejecta exponent and a reach over 1 radius");
35			}
36		14	}
37
38		14	ComplexCraterProfile::ComplexCraterProfile(std::shared_ptr<const CraterProfile> simple,
39		14	ComplexCraterShape shape)
40		14	: simple_(std::move(simple)), shape_(shape) {
41	2/2 ✓ Branch 1 taken 1 time. ✓ Branch 2 taken 13 times.	14	if (!simple_) {
42	1/1 ✓ Branch 2 taken 1 time.	1	throw std::invalid_argument("A complex crater profile needs a simple one");
43			}
44	1/1 ✓ Branch 1 taken 13 times.	13	ValidateComplexCraterShape(shape_);
45		15	}
46
47		264620	double ComplexCraterProfile::ReachRadii() const {
48		264620	return std::fmax(simple_->ReachRadii(), shape_.reachRadii);
49			}
50
51		104113	double ComplexCraterProfile::Height(double distance, double diameter,
52			const CraterState& state) const {
53	2/2 ✓ Branch 1 taken 53223 times. ✓ Branch 2 taken 50890 times.	104113	if (!IsComplex(diameter)) {
54		53223	return simple_->Height(distance, diameter, state);
55			}
56		50890	const double radius = 0.5 * diameter;
57		50890	const double x = distance / radius;
58		50890	const double depth = state.current.depthRatio * diameter;
59		50890	const double rim = state.current.rimHeightRatio * diameter;
60	2/2 ✓ Branch 0 taken 30540 times. ✓ Branch 1 taken 20350 times.	50890	if (x >= 1.0) {
61		30540	const double edge = std::pow(shape_.reachRadii, -shape_.ejectaExponent);
62	2/2 ✓ Branch 0 taken 1 time. ✓ Branch 1 taken 30539 times.	30540	if (x >= shape_.reachRadii) {
63		1	return 0.0;
64			}
65		30539	return rim * (std::pow(x, -shape_.ejectaExponent) - edge) / (1.0 - edge);
66			}
67		20350	const double floor = std::min(0.9, shape_.floorDiameter.At(diameter) / diameter);
68		20350	double height = rim - depth;
69	2/2 ✓ Branch 0 taken 12912 times. ✓ Branch 1 taken 7438 times.	20350	if (x > floor) {
70		12912	const double t = (x - floor) / (1.0 - floor);
71		12912	height += depth * t * t;
72			}
73		20350	const double fitted =
74	1/2 ✓ Branch 0 taken 20350 times. ✗ Branch 1 not taken.	20350	shape_.peakFitMaxDiameter > 0.0 ? std::fmin(diameter, shape_.peakFitMaxDiameter) : diameter;
75		20350	const double peakHeight = shape_.peakHeight.At(fitted);
76	2/2 ✓ Branch 1 taken 7664 times. ✓ Branch 2 taken 12686 times.	20350	const double base = (shape_.peakDiameter ? 0.5 * shape_.peakDiameter->At(fitted)
77		12686	: std::sqrt(3.0 * shape_.peakVolume.At(fitted, true) /
78		12686	(std::numbers::pi * peakHeight))) *
79		20350	(diameter / fitted);
80	2/2 ✓ Branch 0 taken 4613 times. ✓ Branch 1 taken 15737 times.	20350	if (distance < base) {
81		4613	const double fresh = state.fresh.depthRatio;
82	1/2 ✓ Branch 0 taken 4613 times. ✗ Branch 1 not taken.	4613	const double kept = fresh > 0.0 ? state.current.depthRatio / fresh : 1.0;
83		4613	const double u = 1.0 - (distance / base) * (distance / base);
84		4613	height += kept * peakHeight * u * u;
85			}
86		20350	return height;
87			}
88
89			} // namespace terraingen
90

TerrainGen coverage