562 lines
22 KiB
Go
562 lines
22 KiB
Go
package planet
|
|
|
|
import (
|
|
"encoding/json"
|
|
"fmt"
|
|
"math"
|
|
"os"
|
|
"path/filepath"
|
|
"sort"
|
|
"time"
|
|
|
|
"salty/terrain/internal/manifest"
|
|
"salty/terrain/internal/overlay"
|
|
)
|
|
|
|
// Plan reads a template and works out what baking it would involve, without eroding anything.
|
|
//
|
|
// It exists because two decisions can wreck an hour-long bake and both are settled before the first erosion
|
|
// step: how the legend read the painting, and how the planet was cut into regions. Looking at them costs
|
|
// about a minute here and an hour if the bake has to be thrown away.
|
|
func Plan(m *manifest.Manifest, outDir string, mapWidth int, log func(string, ...any)) (*Inputs, error) {
|
|
return PlanPainting(m, nil, outDir, mapWidth, log)
|
|
}
|
|
|
|
// WriteMaps draws the four diagnostic maps. Separate from Plan because a caller that already has an Inputs
|
|
// may want to redraw them without preparing again: the studio does, when only the legend's numbers changed
|
|
// and so only the colouring-in can differ.
|
|
func WriteMaps(outDir string, in *Inputs, mapWidth int) error {
|
|
if err := os.MkdirAll(outDir, 0o755); err != nil {
|
|
return err
|
|
}
|
|
for _, w := range []func(string, *Inputs, int) error{
|
|
WriteClassMap, WriteRegionMap, WriteUpliftMap, WriteErodibilityMap, WriteOverlayMap, WritePlateMap,
|
|
} {
|
|
if err := w(outDir, in, mapWidth); err != nil {
|
|
return err
|
|
}
|
|
}
|
|
return nil
|
|
}
|
|
|
|
// MapNames is what WriteMaps wrote, which is what the studio lists as buttons. The overlay map is there only
|
|
// when there is an overlay, because a map of nothing is a map nobody should be offered.
|
|
func (in *Inputs) MapNames() []string {
|
|
out := []string{"map_class", "map_uplift", "map_regions", "map_erodibility"}
|
|
if in.OverlayRaster != nil {
|
|
out = append(out, "map_overlay")
|
|
}
|
|
if in.Plates != nil {
|
|
out = append(out, "map_plates")
|
|
}
|
|
return out
|
|
}
|
|
|
|
// PlanPainting is Plan over paintings already in memory. See PrepareWith.
|
|
func PlanPainting(m *manifest.Manifest, art *Painting, outDir string, mapWidth int,
|
|
log func(string, ...any)) (*Inputs, error) {
|
|
|
|
in, err := PrepareWith(m, art, log)
|
|
if err != nil {
|
|
return nil, err
|
|
}
|
|
|
|
if err := WriteMaps(outDir, in, mapWidth); err != nil {
|
|
return nil, err
|
|
}
|
|
if in.OverlayDoc != nil {
|
|
if err := in.OverlayDoc.WriteJSON(outDir); err != nil {
|
|
return nil, err
|
|
}
|
|
}
|
|
|
|
rep := in.Report()
|
|
data, err := json.MarshalIndent(rep, "", " ")
|
|
if err != nil {
|
|
return nil, err
|
|
}
|
|
if err := os.WriteFile(filepath.Join(outDir, "plan.json"), append(data, '\n'), 0o644); err != nil {
|
|
return nil, err
|
|
}
|
|
return in, nil
|
|
}
|
|
|
|
// Report is the machine-readable half of a plan, written to plan.json beside the maps.
|
|
type Report struct {
|
|
Manifest string `json:"manifest"`
|
|
Template string `json:"template"`
|
|
Legend string `json:"legend"`
|
|
When time.Time `json:"when"`
|
|
|
|
PaintW int `json:"paint_w"`
|
|
PaintH int `json:"paint_h"`
|
|
|
|
CircumferenceKm float64 `json:"circumference_km"`
|
|
HeightKm float64 `json:"height_km"`
|
|
AreaKm2 float64 `json:"area_km2"`
|
|
CellM float64 `json:"cell_m"`
|
|
GridW int `json:"grid_w"`
|
|
GridH int `json:"grid_h"`
|
|
PadRows int `json:"pad_rows"`
|
|
MarginCells int `json:"margin_cells"`
|
|
TalusDeg float64 `json:"talus_deg"`
|
|
ClampCeilMmYr float64 `json:"clamp_ceiling_mm_yr"`
|
|
ElevationMinM float64 `json:"elevation_min_m"`
|
|
ElevationMaxM float64 `json:"elevation_max_m"`
|
|
Seed int64 `json:"seed"`
|
|
|
|
// MassifWavelengthKm is the upland fabric's size after rounding to a whole number of lattice cells, and
|
|
// zero when no class asked for one.
|
|
MassifWavelengthKm float64 `json:"massif_wavelength_km,omitempty"`
|
|
|
|
// The two fields a seed re-rolls that the painting does not fix: the rock provinces and the fault set.
|
|
// Zero wavelength means the planet asks for none of that field at all.
|
|
LithologyKm float64 `json:"lithology_wavelength_km,omitempty"`
|
|
LithologyTypes int `json:"lithology_types,omitempty"`
|
|
FaultGrainKm float64 `json:"fault_grain_km,omitempty"`
|
|
FaultCount int `json:"faults,omitempty"`
|
|
|
|
MatchFar int `json:"match_far_px"`
|
|
MatchWorst float64 `json:"match_worst_distance"`
|
|
MatchWorstAt [2]int `json:"match_worst_at"`
|
|
EdgeRescued int `json:"stroke_rescued_at_poles_px"`
|
|
Dissolved int `json:"stroke_dissolved_px"`
|
|
WrapRows int `json:"wrap_rows"`
|
|
WrapDiffer int `json:"wrap_differ"`
|
|
WrapLandSea int `json:"wrap_land_against_sea"`
|
|
WrapFarEdge int `json:"wrap_far_edge_px"`
|
|
|
|
Classes []ClassShare `json:"classes"`
|
|
Regions []RegionPlan `json:"regions"`
|
|
|
|
// Overlay is the annotation layer's share of the plan, or nil when the planet has none. The full
|
|
// document - every feature in world metres - goes to overlay.json beside the maps; this is the summary.
|
|
Overlay *OverlayShare `json:"overlay,omitempty"`
|
|
|
|
LandCells int `json:"land_cells"`
|
|
SolveCells int `json:"solve_cells"`
|
|
DroppedRegions int `json:"dropped_regions"`
|
|
DroppedCells int `json:"dropped_cells"`
|
|
EstimateMin float64 `json:"estimate_minutes"`
|
|
EstimatePeakGB float64 `json:"estimate_peak_gb"`
|
|
PrepareSeconds float64 `json:"prepare_seconds"`
|
|
}
|
|
|
|
type ClassShare struct {
|
|
Name string `json:"name"`
|
|
Sea bool `json:"sea"`
|
|
Cells int `json:"cells"`
|
|
Share float64 `json:"share"`
|
|
UpliftMmYr float64 `json:"uplift_mm_yr,omitempty"`
|
|
KMult float64 `json:"k_mult,omitempty"`
|
|
DepthM float64 `json:"depth_m,omitempty"`
|
|
|
|
// DivideDeg is the hillslope angle this class's numbers imply at a divide, and Clamped says whether
|
|
// that is past the angle of repose. See the note on divideAngle.
|
|
DivideDeg float64 `json:"divide_deg,omitempty"`
|
|
Clamped bool `json:"clamped,omitempty"`
|
|
|
|
// MedianDeg and P90Deg are what the ground actually comes out as: the median slope over the class and
|
|
// the ninetieth percentile. See typicalFromDivide - a divide is the *steepest* place in a catchment and
|
|
// there are very few of them, so the divide angle is about three times the ground, and an author reading
|
|
// it as the landscape sets every rate they own two or three times too hot.
|
|
MedianDeg float64 `json:"median_deg,omitempty"`
|
|
P90Deg float64 `json:"p90_deg,omitempty"`
|
|
|
|
// ReadsAs names the ground that angle makes. It is here because an uplift rate does not look like
|
|
// anything, and reading one as terrain is the mistake the massif field exists to undo. It is taken from
|
|
// the median rather than from the divide, because "what does this read as" is a question about the
|
|
// ground somebody is standing on.
|
|
ReadsAs string `json:"reads_as,omitempty"`
|
|
|
|
// The massif block, when this class has one: the plain between the massifs, the angle *it* makes, and
|
|
// how much of the class stands above the midpoint of the two. Absent for a class that is one rate all
|
|
// over, which is what every class was before D-55.
|
|
FloorMmYr float64 `json:"floor_mm_yr,omitempty"`
|
|
FloorDeg float64 `json:"floor_divide_deg,omitempty"`
|
|
FloorMedianDeg float64 `json:"floor_median_deg,omitempty"`
|
|
FloorReadsAs string `json:"floor_reads_as,omitempty"`
|
|
MassifFraction float64 `json:"massif_fraction,omitempty"`
|
|
|
|
// Faults is how many traces landed in this class's ground, and FaultThrowM the range they were drawn
|
|
// from. Zero when the class asked for none.
|
|
Faults int `json:"faults,omitempty"`
|
|
FaultThrowM [2]float64 `json:"fault_throw_m,omitempty"`
|
|
|
|
// LithologyMix is how much of the planet's rock field this class lets through. Reported even at 1, which
|
|
// is the default, because the useful reading is the column rather than one entry in it.
|
|
LithologyMix float64 `json:"lithology_mix,omitempty"`
|
|
}
|
|
|
|
// OverlayShare is how much of the world the annotation layer covers and what it asked for.
|
|
type OverlayShare struct {
|
|
Legend string `json:"legend"`
|
|
PaintedPx int `json:"painted_px"`
|
|
FarPx int `json:"far_px"`
|
|
Features int `json:"features"`
|
|
Marks []overlay.MarkShare `json:"marks"`
|
|
}
|
|
|
|
type RegionPlan struct {
|
|
ID int `json:"id"`
|
|
X0 int `json:"x0"`
|
|
Y0 int `json:"y0"`
|
|
W int `json:"w"`
|
|
H int `json:"h"`
|
|
WidthKm float64 `json:"width_km"`
|
|
HeightKm float64 `json:"height_km"`
|
|
Cells int `json:"cells"`
|
|
LandCells int `json:"land_cells"`
|
|
Seam bool `json:"seam"`
|
|
EstimateMin float64 `json:"estimate_minutes"`
|
|
EstimateGB float64 `json:"estimate_gb"`
|
|
}
|
|
|
|
// Report gathers everything the plan knows.
|
|
func (in *Inputs) Report() *Report {
|
|
perClass, land, total := in.Map.Counts()
|
|
r := &Report{
|
|
Manifest: in.M.Path, Template: in.M.Planet.Template, Legend: in.M.Planet.Legend,
|
|
When: time.Now().UTC().Truncate(time.Second),
|
|
PaintW: in.PaintW, PaintH: in.PaintH,
|
|
|
|
CircumferenceKm: in.P.CircumferenceM() / 1000,
|
|
HeightKm: in.P.HeightM() / 1000,
|
|
AreaKm2: in.P.CircumferenceM() * in.P.HeightM() / 1e6,
|
|
CellM: in.P.CellM,
|
|
GridW: in.P.W, GridH: in.P.PaintH(),
|
|
PadRows: in.P.PadY,
|
|
MarginCells: in.MarginCells,
|
|
TalusDeg: in.M.Pipeline.Thermal.TalusDeg,
|
|
ClampCeilMmYr: clampCeiling(in.M.Pipeline.Fluvial.K, in.P.CellM,
|
|
in.M.Pipeline.Fluvial.M, in.M.Pipeline.Thermal.TalusDeg),
|
|
ElevationMinM: in.M.ElevationM.Min,
|
|
ElevationMaxM: in.M.ElevationM.Max,
|
|
Seed: in.M.Source.Seed,
|
|
|
|
MatchFar: in.Match.Far, MatchWorst: in.Match.MaxDist, MatchWorstAt: in.Match.MaxAt,
|
|
EdgeRescued: in.EdgeRewritten, Dissolved: in.Dissolved,
|
|
WrapRows: in.Match.WrapRows, WrapDiffer: in.Match.WrapDiffer,
|
|
WrapLandSea: in.Match.WrapLandSea, WrapFarEdge: in.Match.WrapFarEdge,
|
|
|
|
LandCells: land,
|
|
SolveCells: in.SolveCells(),
|
|
DroppedRegions: in.Part.DroppedRegions,
|
|
DroppedCells: in.Part.DroppedCells,
|
|
PrepareSeconds: in.Elapsed.Seconds(),
|
|
}
|
|
|
|
if in.Legend.HasMassifs() {
|
|
r.MassifWavelengthKm = in.M.Planet.MassifWavelengthRoundedKm()
|
|
}
|
|
|
|
for i, c := range in.Legend.Classes {
|
|
if perClass[i] == 0 {
|
|
continue
|
|
}
|
|
cs := ClassShare{Name: c.Name, Sea: c.Sea, Cells: perClass[i],
|
|
Share: float64(perClass[i]) / float64(total)}
|
|
if c.Land() {
|
|
cs.UpliftMmYr, cs.KMult = c.UpliftMmYr, c.K()
|
|
cs.DivideDeg = divideAngle(c.RateMYr(), in.M.Pipeline.Fluvial.K*c.K(),
|
|
in.P.CellM, in.M.Pipeline.Fluvial.M)
|
|
cs.Clamped = cs.DivideDeg >= in.M.Pipeline.Thermal.TalusDeg
|
|
cs.MedianDeg, cs.P90Deg = typicalFromDivide(cs.DivideDeg)
|
|
cs.ReadsAs = readsAs(cs.MedianDeg)
|
|
if fr := c.MassifFraction(); fr > 0 {
|
|
cs.MassifFraction = fr
|
|
cs.FloorMmYr = c.Massif.FloorMmYr
|
|
cs.FloorDeg = divideAngle(c.MassifFloorMYr(), in.M.Pipeline.Fluvial.K*c.K(),
|
|
in.P.CellM, in.M.Pipeline.Fluvial.M)
|
|
cs.FloorMedianDeg, _ = typicalFromDivide(cs.FloorDeg)
|
|
cs.FloorReadsAs = readsAs(cs.FloorMedianDeg)
|
|
}
|
|
cs.LithologyMix = c.LithMix()
|
|
if c.Faults != nil {
|
|
cs.FaultThrowM = c.ThrowM()
|
|
for _, f := range in.Faults {
|
|
if f.Class == i {
|
|
cs.Faults++
|
|
}
|
|
}
|
|
}
|
|
} else {
|
|
cs.DepthM = c.DepthM
|
|
}
|
|
r.Classes = append(r.Classes, cs)
|
|
}
|
|
|
|
r.FaultGrainKm = in.M.Planet.FaultGrainKm
|
|
r.FaultCount = len(in.Faults)
|
|
if cells := in.M.Planet.LithologyCells(); cells > 0 && in.Legend.HasLithology() {
|
|
r.LithologyKm = in.M.Planet.NoisePeriodKm / float64(cells)
|
|
r.LithologyTypes = len(in.M.Pipeline.Lithology.KMultipliers)
|
|
}
|
|
|
|
if d := in.OverlayDoc; d != nil {
|
|
r.Overlay = &OverlayShare{
|
|
Legend: in.M.Planet.OverlayLegend, PaintedPx: in.OverlayMatch.Total - in.OverlayMatch.Blank,
|
|
FarPx: in.OverlayMatch.Far, Features: len(d.Features), Marks: d.Marks,
|
|
}
|
|
}
|
|
|
|
peak := 0.0
|
|
for _, rg := range in.Part.Regions {
|
|
mins := in.EstimateSeconds(rg.Cells()) / 60
|
|
gb := float64(in.EstimateBytes(rg.Cells())) / (1 << 30)
|
|
if gb > peak {
|
|
peak = gb
|
|
}
|
|
r.EstimateMin += mins
|
|
r.Regions = append(r.Regions, RegionPlan{
|
|
ID: rg.ID, X0: rg.Frame.X0, Y0: rg.Frame.Y0, W: rg.Frame.W, H: rg.Frame.H,
|
|
WidthKm: float64(rg.Frame.W) * in.P.CellM / 1000,
|
|
HeightKm: float64(rg.Frame.H) * in.P.CellM / 1000,
|
|
Cells: rg.Cells(), LandCells: rg.LandCells, Seam: rg.Seam,
|
|
EstimateMin: mins, EstimateGB: gb,
|
|
})
|
|
}
|
|
r.EstimatePeakGB = peak
|
|
sort.Slice(r.Regions, func(a, b int) bool { return r.Regions[a].Cells > r.Regions[b].Cells })
|
|
return r
|
|
}
|
|
|
|
// Print is the human half: the two tables worth reading before spending an hour.
|
|
func (r *Report) Print(w *os.File) {
|
|
p := func(format string, a ...any) { fmt.Fprintf(w, format+"\n", a...) }
|
|
|
|
p("")
|
|
p(" planet %.1f x %.1f km, %.0f km2 - %d x %d cells of %.1f m (+%d rows of polar pad)",
|
|
r.CircumferenceKm, r.HeightKm, r.AreaKm2, r.GridW, r.GridH, r.CellM, r.PadRows)
|
|
p(" template %d x %d px, %.2f m a pixel - the paint is %s than the grid",
|
|
r.PaintW, r.PaintH, r.CircumferenceKm*1000/float64(r.PaintW),
|
|
coarserOrFiner(r.CircumferenceKm*1000/float64(r.PaintW), r.CellM))
|
|
p(" classify %d px further than the warn distance from any class (worst %.0f at %d,%d)",
|
|
r.MatchFar, r.MatchWorst, r.MatchWorstAt[0], r.MatchWorstAt[1])
|
|
p(" %d px rescued as map-edge class, %d px of stroke dissolved",
|
|
r.EdgeRescued, r.Dissolved)
|
|
if r.WrapRows > 0 {
|
|
pct := 100 * float64(r.WrapDiffer) / float64(r.WrapRows)
|
|
p(" wrap the left and right edges are the same meridian: they disagree on %d of %d rows (%.1f%%),",
|
|
r.WrapDiffer, r.WrapRows, pct)
|
|
p(" %d of those land against water, and %d px in the outermost columns match no class.",
|
|
r.WrapLandSea, r.WrapFarEdge)
|
|
if r.WrapLandSea > r.WrapRows/50 {
|
|
p(" THAT IS A VISIBLE SEAM. The generator wraps; the painting has to as well.")
|
|
}
|
|
}
|
|
p("")
|
|
p(" class share cells uplift mm/yr K depth m divide typical")
|
|
clamped := 0
|
|
for _, c := range r.Classes {
|
|
if c.Sea {
|
|
p(" %-12s %5.1f%% %10d - - %7.0f", c.Name, 100*c.Share, c.Cells, c.DepthM)
|
|
continue
|
|
}
|
|
note := ""
|
|
if c.Clamped {
|
|
note = " CLAMPED"
|
|
clamped++
|
|
}
|
|
p(" %-12s %5.1f%% %10d %8.3f %4.2f - %5.1f deg %6.1f deg %s%s",
|
|
c.Name, 100*c.Share, c.Cells, c.UpliftMmYr, c.KMult, c.DivideDeg, c.MedianDeg, c.ReadsAs, note)
|
|
if c.MassifFraction > 0 {
|
|
p(" %-12s %s", "",
|
|
fmt.Sprintf("massif over %.0f%% of it; the other %.0f%% is %.3f mm/yr, %.1f deg at a divide "+
|
|
"and %.1f typical - %s",
|
|
100*c.MassifFraction, 100*(1-c.MassifFraction), c.FloorMmYr, c.FloorDeg,
|
|
c.FloorMedianDeg, c.FloorReadsAs))
|
|
}
|
|
}
|
|
if clamped > 0 {
|
|
p("")
|
|
p(" %d class(es) sit past the %.0f degree angle of repose at a divide, so the repose clamp shapes",
|
|
clamped, r.TalusDeg)
|
|
p(" them rather than erosion does, and the ground comes out as flat polygonal facets cut along the")
|
|
p(" eight D8 directions. Steady state is S = U/(K*A^m) applied down to a single cell, so at a %.0f m",
|
|
r.CellM)
|
|
p(" cell the ceiling is U = tan(talus)*K*cell = %.3f mm/yr at K x1. Above it, relief and steepness", r.ClampCeilMmYr)
|
|
p(" are the same knob and you get talus, not mountains. See Terrain-Next 4.B1 and 4.D.3.")
|
|
}
|
|
p("")
|
|
p(" `divide` is the steepest ground a rate can make and `typical` is the median over the class, which")
|
|
p(" is about a third of it: a divide is the top of a catchment and there are very few of them. Read the")
|
|
p(" second column. Reading the first as the landscape is how a legend ends up two or three times too")
|
|
p(" hot everywhere, which is the defect the massif block was added to undo one size up.")
|
|
if r.FaultCount > 0 || r.LithologyKm > 0 {
|
|
p("")
|
|
p(" what the seed re-rolls, and the painting does not")
|
|
if r.LithologyKm > 0 {
|
|
p(" lithology %d rock types over provinces of %.1f km, cut at quantiles of the *planet* so "+
|
|
"every", r.LithologyTypes, r.LithologyKm)
|
|
p(" region agrees; it multiplies each class's own k_mult by its lithology_mix")
|
|
}
|
|
if r.FaultCount > 0 {
|
|
p(" faults %d traces, strike from a %.0f km grain field. A trace is a rate difference "+
|
|
"across a", r.FaultCount, r.FaultGrainKm)
|
|
p(" line, steep one side and gentle the other, which erosion carves into a scarp")
|
|
for _, c := range r.Classes {
|
|
if c.Faults > 0 {
|
|
word := "traces"
|
|
if c.Faults == 1 {
|
|
word = "trace"
|
|
}
|
|
p(" %-12s %4d %-7s throw %.0f..%.0f m over the run",
|
|
c.Name, c.Faults, word+",", c.FaultThrowM[0], c.FaultThrowM[1])
|
|
}
|
|
}
|
|
}
|
|
p(" change source.seed, or pass --seed, and all of it moves while the painting stays put")
|
|
}
|
|
if r.MassifWavelengthKm > 0 {
|
|
p("")
|
|
p(" the massif fabric is %.1f km and is one field for the whole planet, so a highland belt and the",
|
|
r.MassifWavelengthKm)
|
|
p(" hills in the lowland beside it are high and low parts of the same structure. A fraction is a")
|
|
p(" share of the planet's surface, so it is only the *expected* share of any one island: a small one")
|
|
p(" may get all of a massif or none, which is the point of not normalising it per landmass.")
|
|
}
|
|
if o := r.Overlay; o != nil {
|
|
p("")
|
|
p(" overlay %s: %d px painted, %d features", o.Legend, o.PaintedPx, o.Features)
|
|
if o.FarPx > 0 {
|
|
p(" %d px are painted but match no mark and were dropped", o.FarPx)
|
|
}
|
|
for _, m := range o.Marks {
|
|
line := fmt.Sprintf(" %-12s %8.2f km2 %4d %s", m.Name, m.AreaKm2, m.Pieces,
|
|
pieces(m.Pieces))
|
|
if m.HasJitter {
|
|
if m.Jitter == 0 {
|
|
line += " coastline pinned as drawn"
|
|
} else {
|
|
line += fmt.Sprintf(" coast jitter x%.2g", m.Jitter)
|
|
}
|
|
}
|
|
if m.Kind == overlay.KindPath && m.WidthM > 0 {
|
|
line += fmt.Sprintf(" %.0f m wide", m.WidthM)
|
|
}
|
|
p("%s", line)
|
|
}
|
|
}
|
|
|
|
p("")
|
|
p(" %d regions, %d cells to solve against %d cells of painted land; margin %d cells",
|
|
len(r.Regions), r.SolveCells, r.LandCells, r.MarginCells)
|
|
if r.DroppedRegions > 0 {
|
|
p(" %d specks dropped, %d land cells, below the minimum", r.DroppedRegions, r.DroppedCells)
|
|
}
|
|
p("")
|
|
p(" id rect km cells land est min est GB")
|
|
shown := r.Regions
|
|
if len(shown) > 12 {
|
|
shown = shown[:12]
|
|
}
|
|
for _, rg := range shown {
|
|
seam := " "
|
|
if rg.Seam {
|
|
seam = "*"
|
|
}
|
|
p(" %3d%s %6.1f x %6.1f %10d %10d %8.1f %6.2f",
|
|
rg.ID, seam, rg.WidthKm, rg.HeightKm, rg.Cells, rg.LandCells, rg.EstimateMin, rg.EstimateGB)
|
|
}
|
|
if len(r.Regions) > len(shown) {
|
|
p(" ... and %d smaller", len(r.Regions)-len(shown))
|
|
}
|
|
p("")
|
|
p(" estimate %.0f min of solve in total, %.2f GB at the largest region, both scaled from one measured",
|
|
r.EstimateMin, r.EstimatePeakGB)
|
|
p(" lowland region and to be read as a floor: steep ground costs about five times what a plain")
|
|
p(" does per cell, because it drives the hillslope law to its full sub-step budget every step.")
|
|
p(" prepared in %.1f s", r.PrepareSeconds)
|
|
p("")
|
|
}
|
|
|
|
func coarserOrFiner(paintM, cellM float64) string {
|
|
if paintM > cellM {
|
|
return "coarser"
|
|
}
|
|
return "finer"
|
|
}
|
|
|
|
// divideAngle is the hillslope angle a class's numbers imply at a drainage divide, in degrees.
|
|
//
|
|
// Steady state is S = U/(K*A^m), and with critical_area_m2 at 0 that law is applied down to a single cell, so
|
|
// at a divide A is one cell squared and A^m is just the cell size. For n = 1 the uplift rate alone therefore
|
|
// fixes the hillslope angle - that is D-49, and it is the most useful number in the whole legend, because it
|
|
// decides whether the ground is shaped by erosion or by landsliding.
|
|
func divideAngle(rateMYr, k, cellM, m float64) float64 {
|
|
if k <= 0 || cellM <= 0 {
|
|
return 0
|
|
}
|
|
s := rateMYr / (k * math.Pow(cellM*cellM, m))
|
|
return math.Atan(s) * 180 / math.Pi
|
|
}
|
|
|
|
// typicalMedianFrac and typicalP90Frac turn a divide angle into the ground underneath it.
|
|
//
|
|
// The divide angle is exact and it is not the landscape. S = U/(K*A^m) is largest where A is smallest, which
|
|
// is the top of a catchment; slope falls away downstream from there, and almost none of a map is divide. So
|
|
// the number the legend hands an author is the steepest place in their world and they read it as the world.
|
|
//
|
|
// Measured rather than derived, on a 600 x 600 grid of 8 m cells with one coast, the manifest's own
|
|
// constants, 1000 steps, and a uniform rate:
|
|
//
|
|
// U mm/yr divide median P90 over 3 deg
|
|
// 0.012 1.7 0.58 1.00 4 %
|
|
// 0.045 6.4 2.12 2.85 8 %
|
|
// 0.080 11.3 3.72 4.85 75 %
|
|
// 0.250 32.0 11.13 14.03 99 %
|
|
//
|
|
// In tangent the median/divide ratio is 0.34, 0.33, 0.33 and 0.32 - flat enough over a factor of twenty in
|
|
// rate to be worth quoting as one number. The P90 ratio drifts from 0.59 to 0.40 as the ground steepens,
|
|
// because the tail of the slope distribution is the part the repose clamp eventually binds; 0.45 is the
|
|
// middle of it and it is the weaker of the two.
|
|
//
|
|
// Both are fractions of the *tangent*, not of the angle, because the steady-state law is about slope.
|
|
const (
|
|
typicalMedianFrac = 0.33
|
|
typicalP90Frac = 0.45
|
|
)
|
|
|
|
// typicalFromDivide is the median and P90 slope, in degrees, for a class whose divide angle is deg.
|
|
func typicalFromDivide(deg float64) (median, p90 float64) {
|
|
t := math.Tan(deg * math.Pi / 180)
|
|
return math.Atan(t*typicalMedianFrac) * 180 / math.Pi,
|
|
math.Atan(t*typicalP90Frac) * 180 / math.Pi
|
|
}
|
|
|
|
func pieces(n int) string {
|
|
if n == 1 {
|
|
return "piece"
|
|
}
|
|
return "pieces"
|
|
}
|
|
|
|
// readsAs names the ground a divide angle makes.
|
|
//
|
|
// It exists because an uplift rate does not look like anything, and the one number an author has to hand is
|
|
// therefore the one they cannot picture. The boundaries are angles rather than rates deliberately: the mistake
|
|
// this is here to stop is reading internal/stats' "plain below 0.1 mm/yr" as a description of terrain. It is
|
|
// not - it is a reporting bucket calibrated for the procedural path's intraplate rates - and 0.1 mm/yr is a 14
|
|
// degree hillslope on every divide of the map, which is hill country wherever it is painted.
|
|
func readsAs(deg float64) string {
|
|
switch {
|
|
case deg < 3:
|
|
return "plain"
|
|
case deg < 8:
|
|
return "rolling"
|
|
case deg < 16:
|
|
return "hill country"
|
|
case deg < 28:
|
|
return "mountain"
|
|
default:
|
|
return "alpine"
|
|
}
|
|
}
|
|
|
|
// clampCeiling is the uplift rate, in mm/yr, at which a divide reaches the angle of repose at K x1. Above it
|
|
// the repose clamp does the shaping and the terrain comes out as flat polygonal facets.
|
|
func clampCeiling(k, cellM, m, talusDeg float64) float64 {
|
|
return math.Tan(talusDeg*math.Pi/180) * k * math.Pow(cellM*cellM, m) * 1000
|
|
}
|