// Package noise is the small toolkit the uplift field is built from: value noise, fBm, domain warping and
// Worley crest lines. Ported from Scripts/Authoring/heightmap_noise.py, and the tuning is the point of the
// port, not the shapes. Read the constants in Docs/Terrain.md before changing any of them; every one of them
// is a round of measurement that has already been paid for.
//
// The one change of intent from the numpy original (D-47): what comes out of here is no longer the terrain.
// It is an uplift rate field that the fluvial pass integrates. Amplitudes are therefore far smaller and the
// shapes matter more than the heights.
package noise

import (
	"math"
	"math/rand/v2"

	"salty/terrain/internal/field"
)

// Source is a per-pass random source. Never the global one: determinism is cross-cutting rule 12 and a
// shared source makes the result depend on which pass drew first.
type Source struct{ r *rand.Rand }

// NewSource seeds from the run seed and the pass index, so inserting a pass does not reshuffle the ones
// before it.
func NewSource(seed int64, pass uint64) *Source {
	return &Source{r: rand.New(rand.NewPCG(uint64(seed), pass))}
}

func (s *Source) Float() float64               { return s.r.Float64() }
func (s *Source) Range(lo, hi float64) float64 { return lo + (hi-lo)*s.r.Float64() }
func (s *Source) IntN(n int) int               { return s.r.IntN(n) }

func Smoothstep(t float64) float64 { return t * t * (3 - 2*t) }

func smoothstep32(t float32) float32 { return t * t * (3 - 2*t) }

// Lattice is one octave's random grid. Periodic: sampling outside it wraps, so warped or stretched
// coordinates never run off an edge.
type Lattice struct {
	Cells int
	Data  []float32
}

// NewLattice draws cells*cells uniforms in index order, which is what makes the octave reproducible however
// it is later sampled in parallel.
func NewLattice(cells int, s *Source) *Lattice {
	l := &Lattice{Cells: cells, Data: make([]float32, cells*cells)}
	for i := range l.Data {
		l.Data[i] = float32(s.Float())
	}
	return l
}

// Sample interpolates the periodic lattice at (u, v) in cell units.
func (l *Lattice) Sample(u, v float64) float32 {
	c := l.Cells
	i0 := int(math.Floor(u))
	j0 := int(math.Floor(v))
	tu := smoothstep32(float32(u - math.Floor(u)))
	tv := smoothstep32(float32(v - math.Floor(v)))
	i0 = ((i0 % c) + c) % c
	j0 = ((j0 % c) + c) % c
	i1 := (i0 + 1) % c
	j1 := (j0 + 1) % c
	top := l.Data[j0*c+i0]*(1-tu) + l.Data[j0*c+i1]*tu
	bot := l.Data[j1*c+i0]*(1-tu) + l.Data[j1*c+i1]*tu
	return top*(1-tv) + bot*tv
}

// Params are one fBm stack. Gain is the constant that matters most: eight octaves at 0.5 makes every octave
// as steep as the last and puts a third of the land above 50 degrees. 0.42 to 0.45.
type Params struct {
	BaseCells int
	Octaves   int
	Gain      float64
	Ridged    bool
}

// FBM fills a size x size field with fractional Brownian motion in [0, 1] on the regular grid.
func FBM(size int, s *Source, p Params) *field.Field {
	uv := identity(size)
	return FBMAt(uv.u, uv.v, s, p)
}

// FBMAt samples the same stack at map coordinates, where 0..1 spans the map once. Feed it warped or
// anisotropic coordinates and the noise bends and stretches with them, which is how the ranges come out as
// long chains rather than blobs.
func FBMAt(u, v *field.Field, s *Source, p Params) *field.Field {
	out := field.NewLike(u)
	amplitude, cells, norm := 1.0, p.BaseCells, 0.0
	// Lattices are built up front, in octave order, before any sampling: the draw order must not depend on
	// the parallel loop below.
	lattices := make([]*Lattice, p.Octaves)
	for o := 0; o < p.Octaves; o++ {
		lattices[o] = NewLattice(cells, s)
		cells *= 2
	}
	cells = p.BaseCells
	for o := 0; o < p.Octaves; o++ {
		l := lattices[o]
		amp := float32(amplitude)
		c := float64(cells)
		field.Rows(u.H, func(y0, y1 int) {
			for y := y0; y < y1; y++ {
				for x := 0; x < u.W; x++ {
					i := y*u.W + x
					n := l.Sample(float64(u.Data[i])*c, float64(v.Data[i])*c)
					if p.Ridged {
						n = 1 - float32(math.Abs(float64(n)*2-1))
						n = n * n
					}
					out.Data[i] += n * amp
				}
			}
		})
		norm += amplitude
		amplitude *= p.Gain
		cells *= 2
	}
	inv := float32(1 / norm)
	for i := range out.Data {
		out.Data[i] *= inv
	}
	return out
}

// CellularEdges is Worley F2-F1 through periodic jittered feature points, mapped so the borders between cells
// read 1 and the interiors 0: a network of thin, branching crest lines.
//
// Kept light on purpose. At 30 % of the mountain height this turned the ranges into a honeycomb of polygon
// walls with flat floors (2026-09-17); 12 % through a stronger warp is the setting that survived.
func CellularEdges(u, v *field.Field, s *Source, cells int, jitter float64) *field.Field {
	pts := make([]float32, cells*cells*2)
	for i := range pts {
		pts[i] = float32(s.Float()*jitter + (1-jitter)*0.5)
	}
	out := field.NewLike(u)
	field.Rows(u.H, func(y0, y1 int) {
		for y := y0; y < y1; y++ {
			for x := 0; x < u.W; x++ {
				i := y*u.W + x
				su := float64(u.Data[i]) * float64(cells)
				sv := float64(v.Data[i]) * float64(cells)
				i0 := int(math.Floor(su))
				j0 := int(math.Floor(sv))
				fu := su - math.Floor(su)
				fv := sv - math.Floor(sv)
				f1, f2 := math.Inf(1), math.Inf(1)
				for dj := -1; dj <= 1; dj++ {
					for di := -1; di <= 1; di++ {
						ci := ((i0+di)%cells + cells) % cells
						cj := ((j0+dj)%cells + cells) % cells
						px := float64(pts[(cj*cells+ci)*2]) + float64(di) - fu
						py := float64(pts[(cj*cells+ci)*2+1]) + float64(dj) - fv
						d := math.Hypot(px, py)
						if d < f1 {
							f2, f1 = f1, d
						} else if d < f2 {
							f2 = d
						}
					}
				}
				e := 1 - (f2-f1)/0.6
				if e < 0 {
					e = 0
				} else if e > 1 {
					e = 1
				}
				out.Data[i] = float32(e * e)
			}
		}
	})
	return out
}

type uvPair struct{ u, v *field.Field }

// identity is the unwarped coordinate pair: 0..1 across the map, matching numpy's mgrid / (size - 1).
func identity(size int) uvPair {
	u := field.New(size, size, 1)
	v := field.New(size, size, 1)
	inv := float32(1) / float32(size-1)
	for y := 0; y < size; y++ {
		for x := 0; x < size; x++ {
			u.Data[y*size+x] = float32(x) * inv
			v.Data[y*size+x] = float32(y) * inv
		}
	}
	return uvPair{u, v}
}

// Identity exposes the coordinate pair so passes can warp it.
func Identity(size int) (u, v *field.Field) {
	p := identity(size)
	return p.u, p.v
}

// Warp offsets a coordinate pair by a low-frequency field scaled by amount. This is what bends ridges so
// ranges curve instead of running straight.
func Warp(u, v, dx, dy *field.Field, amount float64) (*field.Field, *field.Field) {
	wu := u.Clone()
	wv := v.Clone()
	a := float32(amount)
	for i := range wu.Data {
		wu.Data[i] += (dx.Data[i] - 0.5) * 2 * a
		wv.Data[i] += (dy.Data[i] - 0.5) * 2 * a
	}
	return wu, wv
}

// WorldUV is the coordinate pair for noise that must not move when the map does: u and v count periods of
// *world* space from the world origin rather than fractions of the map.
//
// Everything else in this package takes coordinates where 0..1 spans the map once, which is right for a field
// whose whole job is to be shaped like the continent — the range grain, the continent outline, the lithology.
// It is wrong for anything that will one day be generated a tile at a time, because two tiles asking about the
// same physical place would get different values and every seam would show. That is rule 1 of the tiling plan
// in Docs/Terrain-Next.md, and this function is what obeying it looks like.
//
// periodM is how far the lattice runs before it repeats, so it must be comfortably larger than any world that
// will ever be generated; the octave count then sets the finest wavelength, periodM / (BaseCells * 2^octaves).
func WorldUV(w, h int, cellM, originXM, originYM, periodM float64) (u, v *field.Field) {
	u = field.New(w, h, cellM)
	v = field.New(w, h, cellM)
	for y := 0; y < h; y++ {
		vy := float32((originYM + float64(y)*cellM) / periodM)
		for x := 0; x < w; x++ {
			i := y*w + x
			u.Data[i] = float32((originXM + float64(x)*cellM) / periodM)
			v.Data[i] = vy
		}
	}
	return u, v
}