Minor cleanup.

Added particle toggle.
Added boundary shape adjustment.

edt/edt.go

@@ -0,0 +1,293 @@
+package edt
+
+import "math"
+
+type Bounds struct {
+	X int
+	Y int
+}
+
+type Coordinate struct {
+	X int
+	Y int
+}
+
+/*
+Represents a Euclidean Distance Transform of dimensions Dimensions with the transform
+stored in D and a lookup in L. The lookup represents the nearest occupied coordinate (x, y)
+relative to the current coordinate.
+*/
+type EDT struct {
+	Dimensions Bounds       //spatial dimensions
+	D          []float64    //distance transform
+	L          []Coordinate //transform lookup
+	gx         []float64    //buffer for our vector function, per row
+	gy         []float64    //buffer for our vector function, per column
+}
+
+func NewEDT(dimensions Bounds) *EDT {
+	edt := &EDT{
+		Dimensions: dimensions,
+	}
+
+	edt.InitializeD()
+	return edt
+}
+
+func (edt *EDT) InitializeD() {
+
+	n := edt.Dimensions.X * edt.Dimensions.Y
+	if n <= 0 {
+		return
+	}
+
+	edt.D = make([]float64, n)
+	edt.L = make([]Coordinate, n)
+	edt.gx = make([]float64, edt.Dimensions.X)
+	edt.gy = make([]float64, edt.Dimensions.Y)
+
+	for i := range n {
+		edt.D[i] = math.Inf(+1)
+	}
+}
+
+/*
+Reinitializes EDT using occupancy values. Occupancy dimensions must match EDT.
+*/
+func (edt *EDT) AssignOccupancy(occupancy []bool) {
+
+	n := edt.Dimensions.X * edt.Dimensions.Y
+	if len(occupancy) != n {
+		return
+	}
+
+	//assign occupancy values
+	for i := range n {
+		if occupancy[i] {
+			edt.D[i] = 0
+		} else {
+			edt.D[i] = math.Inf(+1)
+		}
+	}
+}
+
+/*
+EDT must already have been initialized with occupancy assigned
+*/
+func (edt *EDT) ComputeDistanceTransform() {
+	if len(edt.D) == 0 {
+		return
+	}
+	//compute edt for each row
+	for j := range edt.Dimensions.Y {
+		//prepare our row vector
+		edt.ConstructGx(j)
+
+		//compute distance transform for this row
+		//rowD, rowL := edt.DistanceTransform1D(edt.gx)
+		rowD, rowL := edt1D_with_labels(edt.gx)
+
+		//update corresponding row
+		edt.UpdateRow(j, rowD, rowL)
+	}
+
+	//compute edt for each column
+	for i := range edt.Dimensions.X {
+		//prepare our column vector
+		edt.ConstructGy(i)
+
+		//compute distance transform for this column
+		//colD, colL := edt.DistanceTransform1D(edt.gy)
+		colD, colL := edt1D_with_labels(edt.gy)
+
+		//update corresponding column
+		edt.UpdateCol(i, colD, colL)
+	}
+}
+
+/*Computes EDT and assigns lookup coordinates for given vector g*/
+func (edt *EDT) DistanceTransform1D(g []float64) (d []float64, l []int) {
+	//perform 1d distance transform
+	n := len(g)
+	d = make([]float64, n)
+	l = make([]int, n)
+
+	v := make([]int, n)
+	// envelope sites (seed indices)
+	z := make([]float64, n+1) // changeover abscissae
+	k := 0
+
+	// top
+	v[0] = 0
+	z[0] = math.Inf(-1)
+	z[1] = math.Inf(+1)
+
+	// initialize first valid seed
+	k = -1
+	for q := 0; q < n; q++ {
+		if g[q] == math.Inf(+1) {
+			continue
+		}
+
+		// push q: pop while it starts before current segment
+		var s float64
+		for {
+			if k < 0 {
+				k = 0
+				v[0] = q
+				z[0] = math.Inf(-1)
+				z[1] = math.Inf(+1)
+				break
+			}
+			i := v[k]
+			// intersection s(i,q)
+			s = ((float64(q*q) + g[q]) - (float64(i*i) + g[i])) / float64(2*(q-i))
+			if s <= z[k] {
+				k--
+				if k < 0 {
+					continue
+				} // force a push with -inf
+			} else {
+				k++
+				v[k] = q
+				z[k] = s
+				z[k+1] = math.Inf(+1)
+				break
+			}
+		}
+	}
+
+	// evaluate
+	k = 0
+	for x := 0; x < n; x++ {
+		for z[k+1] <= float64(x) {
+			k++
+		}
+		i := v[k]
+		l[x] = i
+		dx := float64(x - i)
+		d[x] = dx*dx + g[i]
+	}
+
+	return d, l
+}
+
+func edt1D_with_labels(g []float64) ([]float64, []int) {
+	n := len(g)
+	d := make([]float64, n)
+	l := make([]int, n)
+
+	// --- 1. Gather valid seeds ---
+	type site struct {
+		i  int
+		gi float64
+	}
+	seeds := make([]site, 0, n)
+	for i := 0; i < n; i++ {
+		if !math.IsInf(g[i], +1) { // it's a seed
+			seeds = append(seeds, site{i, g[i]})
+		}
+	}
+
+	// --- 2. Handle empty case directly ---
+	if len(seeds) == 0 {
+		for i := 0; i < n; i++ {
+			d[i] = math.Inf(+1)
+			l[i] = -1
+		}
+		return d, l
+	}
+
+	// --- 3. Lower-envelope construction (Felzenszwalb–Huttenlocher) ---
+	v := make([]int, len(seeds))
+	z := make([]float64, len(seeds)+1)
+	k := 0
+	v[0] = seeds[0].i
+	z[0] = math.Inf(-1)
+	z[1] = math.Inf(+1)
+
+	for q := 1; q < len(seeds); q++ {
+		i := seeds[q].i
+		gi := seeds[q].gi
+		for {
+			j := v[k]
+			gj := g[j]
+			s := ((float64(i*i) + gi) - (float64(j*j) + gj)) / float64(2*(i-j))
+			if s <= z[k] {
+				k--
+				if k < 0 {
+					k = 0
+					break
+				}
+				continue
+			}
+			k++
+			v[k] = i
+			z[k] = s
+			z[k+1] = math.Inf(+1)
+			break
+		}
+	}
+
+	// --- 4. Evaluate ---
+	k = 0
+	for x := 0; x < n; x++ {
+		for z[k+1] <= float64(x) {
+			k++
+		}
+		i := v[k]
+		l[x] = i
+		dx := float64(x - i)
+		d[x] = dx*dx + g[i]
+	}
+
+	return d, l
+}
+
+/*
+Prepares row vector (gx) from specified column in D
+*/
+func (edt *EDT) ConstructGx(col int) {
+	if col < edt.Dimensions.Y {
+		//we can be very efficient due to slice referencing
+		startIdx := col * edt.Dimensions.Y
+		edt.gx = edt.D[startIdx : startIdx+edt.Dimensions.X]
+	}
+}
+
+/*
+Prepares column vector (gy) from specified row in D
+*/
+func (edt *EDT) ConstructGy(row int) {
+	if row < edt.Dimensions.X {
+		for j := range edt.Dimensions.Y {
+			edt.gy[j] = edt.D[j*edt.Dimensions.X+row]
+		}
+	}
+}
+
+/*
+Writes data and lookup values into EDT row specified by rowidx
+*/
+func (edt *EDT) UpdateRow(rowidx int, data []float64, lookup []int) {
+	if rowidx < edt.Dimensions.Y {
+		startidx := rowidx * edt.Dimensions.X
+		for i := range edt.Dimensions.X {
+			edt.D[startidx+i] = data[i]
+			edt.L[startidx+i].X = lookup[i]
+		}
+	}
+}
+
+/*
+Writes data and lookup values into EDT column specified by colidx
+*/
+func (edt *EDT) UpdateCol(colidx int, data []float64, lookup []int) {
+	if colidx < edt.Dimensions.X {
+		for j := range edt.Dimensions.Y {
+			edt.D[j*edt.Dimensions.X+colidx] = data[j]
+			edt.L[j*edt.Dimensions.X+colidx].X = edt.L[lookup[j]*edt.Dimensions.X+colidx].X
+			edt.L[j*edt.Dimensions.X+colidx].Y = lookup[j]
+		}
+	}
+}

elements/alert.go

@@ -0,0 +1,49 @@
+package elements
+
+import (
+	"fluids/fonts"
+
+	"github.com/hajimehoshi/ebiten/v2"
+	"github.com/hajimehoshi/ebiten/v2/text/v2"
+)
+
+const (
+	AlertWidth  = 50
+	AlertHeight = 150
+)
+
+type Alert struct {
+	Sprite         *ebiten.Image
+	text           string
+	textfacesource *text.GoTextFaceSource
+}
+
+func NewAlert() *Alert {
+	a := &Alert{
+		Sprite:         ebiten.NewImage(AlertWidth, AlertHeight),
+		textfacesource: fonts.PixelFont,
+	}
+	return a
+}
+
+func (a *Alert) Draw() {
+	a.Sprite.Clear()
+
+	fnt := &text.GoTextFace{
+		Source: a.textfacesource,
+		Size:   20,
+	}
+
+	_, h := text.Measure(a.text, fnt, 0)
+	top := &text.DrawOptions{}
+	top.GeoM.Translate(0, h)
+	text.Draw(a.Sprite, a.text, fnt, top)
+}
+
+func (a *Alert) Update() {
+
+}
+
+func (a *Alert) SetText(t string) {
+	a.text = t
+}

elements/flaskRoundedBottom.go

@@ -0,0 +1,265 @@
+package elements
+
+import (
+	"fluids/fluid"
+	"fluids/gamedata"
+	"fluids/resources"
+	"image/color"
+
+	"github.com/hajimehoshi/ebiten/v2"
+)
+
+const (
+	RoundedBottomFlaskWidth           = 32  //pixels
+	RoundedBottomFlaskHeight          = 32  //pixels
+	RBFlaskMaskWidth                  = 46  //pixels
+	RBFlaskMaskHeight                 = 46  //pixels
+	RoundedBottomFlaskFluidRadius     = 9   //pixels: represents spherical portion of the flask where fluid will be contained
+	RoundedBottomFlaskFluidOriginX    = 16  //pixels: x origin of fluid area
+	RoundedBottomFlaskFluidOriginY    = 20  //pixels: y origin of fluid area
+	RoundedBottomFlaskFluidWidth      = 0.5 //meters
+	RoundedBottomFlaskFluidHeight     = 0.5 //meters
+	RoundedBottomFlaskFluidResolution = 20
+)
+
+type RoundedBottomFlask struct {
+	MappedEntityBase
+	fluid             *fluid.Fluid    //our physical representation of the fluid
+	fluidbuffS        *ebiten.Image   //predraw for the fluid, static
+	fluidbuffD        *ebiten.Image   //predraw for the fluid, dynamic
+	fluidcellbuff     *ebiten.Image   //persistent fluid sprite, we redraw this everywhere we want to represent fluid
+	flaskbase         *ebiten.Image   //flask background (container)
+	flaskhighlight    *ebiten.Image   //flask foreground (glassware highlights)
+	flaskboundarymask *ebiten.Image   //flask boundary mask
+	fieldscaleStatic  gamedata.Vector //used for transforming from fluid-space to sprite-space: static
+	fieldscaleDyanmic gamedata.Vector //used for transforming from fluid-space to sprite-space: dynamic
+	angle             float64
+
+	//fluid color business
+	fluidcolor  color.RGBA //premultiplied fluid color, as set by external sources
+	fluidcolorF []float32  //for caching of individual color values, we compute rarely
+
+	//boundary
+	boundaryinitialized bool
+	container           *fluid.Boundary
+}
+
+func NewRoundedBottomFlask() *RoundedBottomFlask {
+	flask := &RoundedBottomFlask{
+		fluidbuffS:        ebiten.NewImage(RoundedBottomFlaskFluidRadius*2, RoundedBottomFlaskFluidRadius*2),
+		fluidbuffD:        ebiten.NewImage(RBFlaskMaskWidth, RBFlaskMaskHeight),
+		fluidcellbuff:     ebiten.NewImage(1, 1),
+		flaskbase:         ebiten.NewImageFromImage(resources.ImageBank[resources.RoundedBottomFlaskBase]),
+		flaskhighlight:    ebiten.NewImageFromImage(resources.ImageBank[resources.RoundedBottomFlaskHighlights]),
+		flaskboundarymask: ebiten.NewImageFromImage(resources.ImageBank[resources.RoundedBottomFlaskBoundaryMap46]),
+		//flaskboundarymask:   ebiten.NewImageFromImage(resources.ImageBank[resources.RoundedBottomFlaskBoundaryMap]),
+		angle:               0,
+		fluidcolorF:         make([]float32, 4), //one field each for R,G,B,A
+		boundaryinitialized: false,
+	}
+	flask.Initialize()
+	return flask
+}
+
+func (flask *RoundedBottomFlask) Initialize() {
+	//prepare our internal data
+	flask.dimensions = gamedata.Vector{
+		X: RoundedBottomFlaskWidth,
+		Y: RoundedBottomFlaskHeight,
+	}
+	flask.Sprite = ebiten.NewImage(RoundedBottomFlaskWidth, RoundedBottomFlaskHeight)
+	flask.fluidcellbuff.Fill(color.White)
+
+	//prepare and initialize the fluid
+	fluiddimensions := fluid.FieldVector{
+		X: RoundedBottomFlaskFluidWidth,
+		Y: RoundedBottomFlaskFluidHeight,
+	}
+	flask.fluid = fluid.NewFluid(fluiddimensions, RoundedBottomFlaskFluidHeight/RoundedBottomFlaskFluidResolution)
+	flask.fluid.Initialize()
+	flask.fluid.Block = false //rounded flask, not a rect volume
+
+	//compute fieldscales using newly created fluid
+	flask.fieldscaleStatic = gamedata.Vector{
+		X: RoundedBottomFlaskFluidRadius * 2 / float64(flask.fluid.Field.Nx-1),
+		Y: RoundedBottomFlaskFluidRadius * 2 / float64(flask.fluid.Field.Ny-1),
+	}
+
+	flask.fieldscaleDyanmic = gamedata.Vector{
+		X: RBFlaskMaskWidth / float64(flask.fluid.Field.Nx-2),
+		Y: RBFlaskMaskHeight / float64(flask.fluid.Field.Ny-2),
+	}
+
+	//setup default fluid color
+	flask.SetFluidColor(color.RGBA{R: 0x0, G: 0x0, B: 0xff, A: 0xff})
+}
+
+func (flask *RoundedBottomFlask) Update() {
+
+	if flask.paused {
+		return
+	}
+
+	flask.fluid.Step()
+}
+
+func (flask *RoundedBottomFlask) Draw() {
+	flask.Sprite.Clear()
+
+	//render flask background
+	flask.Sprite.DrawImage(flask.flaskbase, nil)
+
+	//render fluid
+	if flask.boundaryinitialized {
+		flask.RenderFluidDynamic()
+	} else {
+		flask.RenderFluidStatic()
+	}
+
+	//render flask foreground
+	flask.Sprite.DrawImage(flask.flaskhighlight, nil)
+
+}
+
+func (flask *RoundedBottomFlask) SetAngle(angle float64) {
+	flask.angle = angle
+	flask.fluid.SetAngle(float32(flask.angle))
+}
+
+func (flask *RoundedBottomFlask) GetAngle() float64 {
+	return flask.angle
+}
+
+func (flask *RoundedBottomFlask) RenderFluidDynamic() {
+	flask.fluidbuffD.Clear()
+	//flask.fluidbuffD.Fill(color.White)
+	//vector.StrokeRect(flask.fluidbuffD, 0, 0, 46, 46, 1, color.White, true)
+
+	//construct fluid buffer from fluid simulation
+	for i := range flask.fluid.Field.Nx {
+		for j := range flask.fluid.Field.Ny {
+
+			idx := i*flask.fluid.Field.Ny + j
+
+			if flask.fluid.Field.CellType[idx] != fluid.CellTypeFluid {
+				continue
+			}
+
+			celldensity := flask.fluid.ParticleDensity[idx] / flask.fluid.ParticleRestDensity
+			/*if celldensity > 0.8 {
+				celldensity = 1
+			}*/
+
+			ox := float64(i) * flask.fieldscaleDyanmic.X
+			oy := float64(j) * flask.fieldscaleDyanmic.Y
+			op := &ebiten.DrawImageOptions{}
+			op.GeoM.Translate(-.5, -.5)
+			op.GeoM.Scale(flask.fieldscaleDyanmic.X, flask.fieldscaleDyanmic.Y)
+			op.GeoM.Translate(ox, oy)
+			op.ColorScale.ScaleAlpha(celldensity)
+			op.ColorScale.Scale(flask.fluidcolorF[0], flask.fluidcolorF[1], flask.fluidcolorF[2], flask.fluidcolorF[3])
+			flask.fluidbuffD.DrawImage(flask.fluidcellbuff, op)
+		}
+	}
+
+	//transform buffer for our flask space
+
+	s := float64(RoundedBottomFlaskWidth) / RBFlaskMaskWidth * 67. / 104
+	op := &ebiten.DrawImageOptions{}
+	op.GeoM.Translate(-RBFlaskMaskWidth/2, -RBFlaskMaskHeight/2)
+	op.GeoM.Scale(s, -s*1.15)
+	op.GeoM.Translate(RoundedBottomFlaskWidth/2, RoundedBottomFlaskHeight/2+2)
+	//op.GeoM.Translate(RoundedBottomFlaskFluidOriginX, RoundedBottomFlaskFluidOriginY)
+	flask.Sprite.DrawImage(flask.fluidbuffD, op)
+}
+
+func (flask *RoundedBottomFlask) RenderFluidStatic() {
+	flask.fluidbuffS.Clear()
+
+	//construct fluid buffer from fluid simulation
+	for i := range flask.fluid.Field.Nx {
+		for j := range flask.fluid.Field.Ny {
+
+			idx := i*flask.fluid.Field.Ny + j
+
+			if flask.fluid.Field.CellType[idx] != fluid.CellTypeFluid {
+				continue
+			}
+
+			celldensity := flask.fluid.ParticleDensity[idx] / flask.fluid.ParticleRestDensity
+			/*if celldensity > 0.8 {
+				celldensity = 1
+			}*/
+
+			ox := float64(i) * flask.fieldscaleStatic.X
+			oy := float64(j) * flask.fieldscaleStatic.Y
+			op := &ebiten.DrawImageOptions{}
+			op.GeoM.Translate(-.5, -.5)
+			op.GeoM.Scale(flask.fieldscaleStatic.X, flask.fieldscaleStatic.Y)
+			op.GeoM.Translate(ox, oy)
+			op.ColorScale.ScaleAlpha(celldensity)
+			op.ColorScale.Scale(flask.fluidcolorF[0], flask.fluidcolorF[1], flask.fluidcolorF[2], flask.fluidcolorF[3])
+			flask.fluidbuffS.DrawImage(flask.fluidcellbuff, op)
+		}
+	}
+
+	//transform buffer for our flask space
+	op := &ebiten.DrawImageOptions{}
+
+	op.GeoM.Translate(-RoundedBottomFlaskFluidRadius, -RoundedBottomFlaskFluidRadius)
+	op.GeoM.Scale(1, -1)
+	op.GeoM.Translate(RoundedBottomFlaskFluidOriginX, RoundedBottomFlaskFluidOriginY)
+	flask.Sprite.DrawImage(flask.fluidbuffS, op)
+}
+
+func (flask *RoundedBottomFlask) SetFluidColor(c color.RGBA) {
+	flask.fluidcolor = c
+	flask.fluidcolorF[0] = float32(flask.fluidcolor.R) / float32(flask.fluidcolor.A)
+	flask.fluidcolorF[1] = float32(flask.fluidcolor.G) / float32(flask.fluidcolor.A)
+	flask.fluidcolorF[2] = float32(flask.fluidcolor.B) / float32(flask.fluidcolor.A)
+	flask.fluidcolorF[3] = float32(flask.fluidcolor.A) / 0xff
+}
+
+func (flask *RoundedBottomFlask) InitializeBoundary() {
+
+	//prepare the dimensions of our boundary map
+	dimensions := fluid.BoundaryDimensions{
+		X: RBFlaskMaskWidth,
+		Y: RBFlaskMaskHeight,
+	}
+
+	//instantiate the boundary structure
+	flask.container = fluid.NewBoundary(dimensions)
+
+	//load pixel data from boundary
+	var pixels []byte = make([]byte, dimensions.X*dimensions.Y*4)
+	flask.flaskboundarymask.ReadPixels(pixels)
+
+	//populate our boundary map based on the pixel data
+	for i := 0; i < len(pixels); i += 4 {
+		cellidx := i / 4
+		boundary := pixels[i] == 0
+		flask.container.Cells[cellidx] = !boundary
+	}
+
+	//apply to fluid simulation
+	flask.fluid.SetBoundary(flask.container)
+	flask.boundaryinitialized = true
+}
+
+// set new boundary mask and reinitialize associated data, including computing and
+// passing that into the associated fluid simulation
+func (flask *RoundedBottomFlask) SetBoundaryMap(img *ebiten.Image) {
+	flask.flaskboundarymask = img
+	flask.InitializeBoundary()
+}
+
+func (flask *RoundedBottomFlask) ToggleBoundaryMask() {
+	if flask.boundaryinitialized {
+		flask.fluid.SetBoundary(nil)
+		flask.boundaryinitialized = false
+	} else {
+		flask.InitializeBoundary()
+		flask.boundaryinitialized = true
+	}
+}

elements/fluidsim10.go

@@ -0,0 +1,146 @@
+package elements
+
+import (
+	"fluids/fluid"
+	"fluids/gamedata"
+	"image/color"
+
+	"github.com/hajimehoshi/ebiten/v2"
+	"github.com/hajimehoshi/ebiten/v2/vector"
+)
+
+const (
+	FS10PixelWidth  = 100
+	FS10PixelHeight = 100
+	FS10FluidWidth  = 1.0 //meters
+	FS10FluidHeight = 1.0 //meters
+	FS10Resolution  = 20  //need to workshop this
+)
+
+type FluidSim10 struct {
+	MappedEntityBase
+	fluid           *fluid.Fluid
+	angle           float64
+	particlebuff    *ebiten.Image
+	renderparticles bool
+	renderfield     bool
+	fieldscale      gamedata.Vector
+}
+
+func NewFluidSim10() *FluidSim10 {
+	fsim := &FluidSim10{
+		fluid:           fluid.NewFluid(fluid.FieldVector{X: FS10FluidWidth, Y: FS10FluidHeight}, FS10FluidHeight/FS10Resolution),
+		particlebuff:    ebiten.NewImage(1, 1),
+		renderfield:     true,  //false,
+		renderparticles: false, //true,
+	}
+	fsim.Initialize()
+	return fsim
+}
+
+func (f *FluidSim10) Initialize() {
+	f.Sprite = ebiten.NewImage(FS10PixelWidth, FS10PixelHeight)
+	f.particlebuff.Fill(color.White)
+	f.fluid.Initialize()
+
+	//fieldscale for rendering purposes
+	f.fieldscale = gamedata.Vector{
+		X: FS10PixelWidth / float64(f.fluid.Field.Nx-1),
+		Y: FS10PixelHeight / float64(f.fluid.Field.Ny-1),
+	}
+}
+
+func (f *FluidSim10) SetAngle(angle float64) {
+	f.angle = angle
+	f.fluid.SetAngle(float32(angle))
+}
+
+func (f *FluidSim10) GetAngle() float64 {
+	return f.angle
+}
+
+func (f *FluidSim10) Draw() {
+	f.Sprite.Clear()
+
+	if f.renderfield {
+		/*
+			alternatively, single loop (not nested) with x/y cell coordinates given by
+			xi := i % f.fluid.Field.Ny
+			yi := i / f.fluid.Field.Ny
+		*/
+		for i := range f.fluid.Field.Nx {
+			for j := range f.fluid.Field.Ny {
+
+				idx := i*f.fluid.Field.Ny + j
+
+				if f.fluid.Field.CellType[idx] != fluid.CellTypeFluid {
+					continue
+				}
+
+				celldensity := f.fluid.ParticleDensity[idx] / f.fluid.ParticleRestDensity
+				if celldensity > 0.8 {
+					celldensity = 1
+				}
+
+				ox := float64(i) * f.fieldscale.X
+				oy := float64(j) * f.fieldscale.Y
+				op := &ebiten.DrawImageOptions{}
+				op.GeoM.Translate(-.5, -.5)
+				op.GeoM.Scale(f.fieldscale.X, f.fieldscale.Y)
+				op.GeoM.Translate(ox, oy)
+				op.ColorScale.ScaleAlpha(celldensity)
+				op.ColorM.Scale(0, 0, 1, 1)
+				f.Sprite.DrawImage(f.particlebuff, op)
+			}
+		}
+	}
+
+	if f.renderparticles {
+		for i := range f.fluid.Particles {
+			//for each particle, compute its relative position based on its
+			//position within the fluid field
+			p := &f.fluid.Particles[i]
+			percentx := (p.Position.X - f.fluid.Field.H/2) / f.fluid.Field.Dimensions.X
+			percenty := (p.Position.Y - f.fluid.Field.H/2) / f.fluid.Field.Dimensions.Y
+			ox := float64(percentx * FS10PixelWidth)
+			oy := float64(percenty * FS10PixelHeight)
+
+			op := &ebiten.DrawImageOptions{}
+			op.GeoM.Translate(ox, oy)
+			f.Sprite.DrawImage(f.particlebuff, op)
+		}
+	}
+
+	if f.fluid.Block {
+		vector.StrokeRect(f.Sprite, 0, 0, FS10PixelWidth, FS10PixelHeight, 2, color.White, true)
+	} else {
+		radius := float32(f.fluid.Field.Nx) * f.fluid.Field.H / 2
+		pixelradius := radius/f.fluid.Field.Dimensions.X*FS10PixelWidth - 4
+		vector.StrokeCircle(f.Sprite, FS10PixelWidth/2, FS10PixelHeight/2, pixelradius, 2, color.White, true)
+	}
+}
+
+func (f *FluidSim10) Update() {
+
+	if f.paused {
+		return
+	}
+
+	f.fluid.Step()
+}
+
+func (f *FluidSim10) EnableParticleRender(v bool) {
+	f.renderparticles = v
+}
+
+func (f *FluidSim10) EnableFieldRender(v bool) {
+	f.renderfield = v
+}
+
+func (f *FluidSim10) ToggleShape() {
+	f.fluid.ToggleShape()
+}
+
+func (f *FluidSim10) ToggleParticles() {
+	f.renderparticles = !f.renderparticles
+}

elements/fluidsimd.go

@@ -0,0 +1,443 @@
+package elements
+
+import (
+	"fluids/gamedata"
+	"fluids/quadtree"
+	"fmt"
+	"image/color"
+	"math"
+	"math/rand/v2"
+
+	"github.com/hajimehoshi/ebiten/v2"
+	"github.com/hajimehoshi/ebiten/v2/vector"
+)
+
+const (
+	FSDWidth           = 640
+	FSDHeight          = 360
+	FSDParticleCount   = 2000
+	FSDGravity         = 2
+	FSDParticleRadius  = 2.5
+	FSDDamping         = .7
+	FSDDeltaTimeStep   = 0.5
+	FSDInfluenceRadius = 30
+)
+
+type FluidSimD struct {
+	MappedEntityBase
+	particles         []*Particle
+	cycle             int
+	particlebox       *gamedata.Vector
+	particlebuff      *ebiten.Image
+	quadtree          *quadtree.Quadtree
+	collisionquad     quadtree.Quadrant
+	paused            bool
+	renderquads       bool
+	resolvecollisions bool
+	resolvers         []func(particle *Particle)
+	resolveridx       int
+	angle             float64
+}
+
+func NewFluidSimD() *FluidSimD {
+	fsd := &FluidSimD{
+		particlebuff:      ebiten.NewImage(FSDParticleRadius*2, FSDParticleRadius*2),
+		paused:            false,
+		renderquads:       false,
+		resolvecollisions: false,
+		resolveridx:       0,
+		angle:             0,
+	}
+	fsd.dimensions = gamedata.Vector{X: FSDWidth, Y: FSDHeight}
+
+	//prepare root quadtree and subsequent collision search quadrant
+	quad := quadtree.Quadrant{
+		Position:   gamedata.Vector{X: FSDWidth / 2, Y: FSDHeight / 2},
+		Dimensions: gamedata.Vector{X: FSDWidth, Y: FSDHeight},
+	}
+	fsd.quadtree = quadtree.New(quad, 0)
+
+	fsd.collisionquad = quadtree.Quadrant{
+		Position: gamedata.Vector{
+			X: 0,
+			Y: 0,
+		},
+		Dimensions: gamedata.Vector{
+			X: FSDParticleRadius * 2,
+			Y: FSDParticleRadius * 2,
+		},
+	}
+
+	//add all resolvers to strategy list
+	fsd.resolvers = append(fsd.resolvers, fsd.ResolveCollisionsA)
+	fsd.resolvers = append(fsd.resolvers, fsd.ResolveCollisionsB)
+	fsd.resolvers = append(fsd.resolvers, fsd.ResolveCollisionsC)
+
+	//initialize particles, set bounding box, prepare image buffer
+	fsd.Sprite = ebiten.NewImage(FSDWidth, FSDHeight)
+	fsd.particlebox = &gamedata.Vector{
+		X: FSDWidth - 50,
+		Y: FSDHeight - 50,
+	}
+	fsd.InitializeParticles()
+
+	vector.FillCircle(fsd.particlebuff, FSDParticleRadius, FSDParticleRadius, FSDParticleRadius, color.White, true)
+
+	return fsd
+
+}
+
+func (f *FluidSimD) Draw() {
+	f.Sprite.Clear()
+
+	f.RenderParticles()
+	f.RenderBox()
+
+	if f.renderquads {
+		f.RenderQuadrants()
+	}
+}
+
+func (f *FluidSimD) Update() {
+
+	if f.paused {
+		return
+	}
+
+	f.RebuildQuadtree()
+	f.UpdateParticles()
+
+	f.cycle++
+}
+
+func (f *FluidSimD) UpdateParticles() {
+
+	dt := FSDDeltaTimeStep
+
+	mx, my := ebiten.CursorPosition()
+	mpos := gamedata.Vector{X: float64(mx), Y: float64(my)}
+	maxdeflect := 40 * dt * FSDGravity
+
+	gravity := gamedata.Vector{
+		X: FSDGravity * dt * math.Sin(f.angle),
+		Y: FSDGravity * dt * math.Cos(f.angle),
+	}
+
+	for _, particle := range f.particles {
+		//particle.Velocity.Y += FSDGravity * dt
+
+		particle.Velocity = particle.Velocity.Add(gravity)
+
+		//if inpututil.IsMouseButtonJustPressed(ebiten.MouseButtonLeft) {
+
+		if ebiten.IsMouseButtonPressed(ebiten.MouseButtonLeft) {
+			delta := gamedata.Vector{
+				X: mpos.X - particle.Position.X,
+				Y: mpos.Y - particle.Position.Y,
+			}
+
+			dist := math.Sqrt(delta.X*delta.X + delta.Y*delta.Y)
+			theta := math.Atan2(delta.Y, delta.X)
+
+			if dist < FSDInfluenceRadius {
+
+				dx := dist * math.Cos(theta)
+				dy := dist * math.Sin(theta)
+
+				if dx != 0 {
+					gainx := (-1./FSDInfluenceRadius)*math.Abs(dx) + 1.
+					particle.Velocity.X += 20 * gainx * -1 * math.Copysign(1, dx)
+				}
+
+				if dy != 0 {
+					gainy := (-1./FSDInfluenceRadius)*math.Abs(dy) + 1.
+					particle.Velocity.Y += maxdeflect * gainy * -1 * math.Copysign(1, dy)
+				}
+			}
+
+		}
+
+		particle.Position.X += particle.Velocity.X * dt
+		particle.Position.Y += particle.Velocity.Y * dt
+
+		if f.resolvecollisions {
+			f.resolvers[f.resolveridx](particle)
+		}
+	}
+
+	for _, p := range f.particles {
+		f.BoundParticle(p)
+	}
+
+}
+
+func (f *FluidSimD) InitializeParticles() {
+
+	f.particles = f.particles[:0]
+
+	xmin := (FSDWidth-f.particlebox.X)/2 + FSDParticleRadius
+	xmax := f.particlebox.X - FSDParticleRadius*2
+	ymin := (FSDHeight-f.particlebox.Y)/2 + FSDParticleRadius
+	ymax := f.particlebox.Y - FSDParticleRadius*2
+
+	for i := 0; i < FSDParticleCount; i++ {
+
+		p := &Particle{
+			Position: gamedata.Vector{
+				X: xmin + rand.Float64()*xmax,
+				Y: ymin + rand.Float64()*ymax,
+			},
+			Velocity: gamedata.Vector{
+				X: 0,
+				Y: 0,
+			},
+			Radius: FSDParticleRadius,
+		}
+
+		f.particles = append(f.particles, p)
+	}
+
+}
+
+func (f *FluidSimD) BoundParticle(p *Particle) {
+
+	xmin := (FSDWidth-f.particlebox.X)/2 + p.Radius
+	xmax := xmin + f.particlebox.X - p.Radius*2
+
+	if p.Position.X > xmax {
+		p.Velocity.X *= -1 * FSDDamping
+		p.Position.X = xmax
+	}
+
+	if p.Position.X < xmin {
+		p.Velocity.X *= -1 * FSDDamping
+		p.Position.X = xmin
+	}
+
+	ymin := (FSDHeight-f.particlebox.Y)/2 + p.Radius
+	ymax := ymin + f.particlebox.Y - p.Radius*2
+
+	if p.Position.Y > ymax {
+		p.Velocity.Y *= -1 * FSDDamping
+		p.Position.Y = ymax
+	}
+
+	if p.Position.Y < ymin {
+		p.Velocity.Y *= -1 * FSDDamping
+		p.Position.Y = ymin
+	}
+
+}
+
+func (f *FluidSimD) RenderQuadrants() {
+	clr := color.RGBA{R: 0xff, G: 0x00, B: 0x00, A: 0xff}
+
+	quadrants := f.quadtree.GetQuadrants()
+	for _, quad := range quadrants {
+		ox := float32(quad.Position.X - quad.Dimensions.X/2)
+		oy := float32(quad.Position.Y - quad.Dimensions.Y/2)
+		vector.StrokeRect(f.Sprite, ox, oy, float32(quad.Dimensions.X), float32(quad.Dimensions.Y), 1, clr, true)
+	}
+
+}
+
+func (f *FluidSimD) RenderParticles() {
+	for _, particle := range f.particles {
+		x0 := particle.Position.X - particle.Radius
+		y0 := particle.Position.Y - particle.Radius
+		op := &ebiten.DrawImageOptions{}
+		op.GeoM.Translate(x0, y0)
+		f.Sprite.DrawImage(f.particlebuff, op)
+	}
+}
+
+func (f *FluidSimD) RenderBox() {
+	x0 := (FSDWidth - f.particlebox.X) / 2
+	y0 := (FSDHeight - f.particlebox.Y) / 2
+	vector.StrokeRect(f.Sprite, float32(x0), float32(y0), float32(f.particlebox.X), float32(f.particlebox.Y), 2, color.White, true)
+}
+
+func (f *FluidSimD) RebuildQuadtree() {
+	f.quadtree.Clear()
+
+	for _, p := range f.particles {
+		if !f.quadtree.Insert(p) {
+			fmt.Println("quadtree insertion failed")
+		}
+	}
+}
+
+func (f *FluidSimD) ResolveCollisionsA(particle *Particle) {
+	//construct search quadrant from current particle
+	quadrant := quadtree.Quadrant{
+		Position:   particle.Position,
+		Dimensions: particle.GetDimensions(),
+	}
+
+	//find list of possible maybe collisions, we inspect those in more detail
+	maybes := f.quadtree.FindAll(quadrant)
+
+	sqdist := float64(particle.Radius*particle.Radius) * 4
+
+	for _, p := range maybes {
+		if p == particle {
+			continue
+		}
+
+		pos := p.GetPosition()
+		delta := gamedata.Vector{
+			X: pos.X - particle.Position.X,
+			Y: pos.Y - particle.Position.Y,
+		}
+		dist2 := delta.X*delta.X + delta.Y*delta.Y
+
+		if dist2 == 0 {
+			// Same position: pick a fallback direction to avoid NaN
+			delta.X = 1
+			delta.Y = 0
+			dist2 = 1
+		}
+
+		if dist2 < sqdist {
+			d := math.Sqrt(dist2)
+			nx, ny := delta.X/d, delta.Y/d
+			overlap := particle.Radius*2 - d
+
+			pos.X += nx * overlap
+			pos.Y += ny * overlap
+			p.SetPosition(pos)
+			/*
+				newpos := gamedata.Vector{
+					X: pos.X + delta.X,
+					Y: pos.Y + delta.Y,
+				}
+				p.SetPosition(newpos)
+			*/
+		}
+	}
+}
+
+func (f *FluidSimD) ResolveCollisionsB(particle *Particle) {
+	//construct search quadrant from current particle
+	quadrant := quadtree.Quadrant{
+		Position:   particle.Position,
+		Dimensions: particle.GetDimensions(),
+	}
+
+	//find list of possible maybe collisions, we inspect those in more detail
+	maybes := f.quadtree.FindAll(quadrant)
+
+	sqdist := float64(particle.Radius*particle.Radius) * 4
+
+	for _, p := range maybes {
+		if p == particle {
+			continue
+		}
+
+		pos := p.GetPosition()
+		delta := gamedata.Vector{
+			X: pos.X - particle.Position.X,
+			Y: pos.Y - particle.Position.Y,
+		}
+
+		dist2 := delta.X*delta.X + delta.Y*delta.Y
+
+		if dist2 < sqdist {
+			d := math.Sqrt(dist2)
+			overlap := particle.Radius*2 - d
+			theta := math.Atan2(delta.Y, delta.X)
+			pos.X += overlap * math.Cos(theta)
+			pos.Y += overlap * math.Sin(theta)
+			p.SetPosition(pos)
+
+			m := p.(*Particle)
+			m.Velocity.X *= -1 * FSDDamping
+			m.Velocity.Y *= -1 * FSDDamping
+
+		}
+	}
+}
+
+func (f *FluidSimD) ResolveCollisionsC(particle *Particle) {
+	//construct search quadrant from current particle
+	quadrant := quadtree.Quadrant{
+		Position:   particle.Position,
+		Dimensions: particle.GetDimensions(),
+	}
+
+	//find list of possible maybe collisions, we inspect those in more detail
+	maybes := f.quadtree.FindAll(quadrant)
+
+	sqdist := float64(particle.Radius*particle.Radius) * 4
+
+	for _, p := range maybes {
+		if p == particle {
+			continue
+		}
+
+		pos := p.GetPosition()
+		delta := gamedata.Vector{
+			X: pos.X - particle.Position.X,
+			Y: pos.Y - particle.Position.Y,
+		}
+
+		dist2 := delta.X*delta.X + delta.Y*delta.Y
+
+		if dist2 < sqdist {
+
+			/*
+				//compute impact to this particle
+				deltav1 := particle.Velocity.Subtract(p.(*Particle).Velocity)
+				deltax1 := particle.Position.Subtract(p.(*Particle).Position)
+				dot1 := deltav1.DotProduct(deltax1)
+				mag1 := deltax1.Magnitude() * deltax1.Magnitude()
+				particle.Velocity = deltax1.Scale(dot1 / mag1)
+
+				//compute impact to other particle
+				deltav2 := p.(*Particle).Velocity.Subtract(particle.Velocity)
+				deltax2 := p.(*Particle).Position.Subtract(particle.Position)
+				dot2 := deltav2.DotProduct(deltax2)
+				mag2 := deltax2.Magnitude() * deltax2.Magnitude()
+				p.(*Particle).Velocity = deltax2.Scale(dot2 / mag2)
+			*/
+
+			dist := math.Sqrt(dist2)
+			s := 0.5 * (particle.Radius*2 - dist) / dist
+
+			dnorm := delta.Scale(s)
+
+			particle.Position = particle.Position.Subtract(dnorm)
+			p.(*Particle).Position = p.(*Particle).Position.Add(dnorm)
+
+		}
+	}
+
+}
+
+func (f *FluidSimD) SetRenderQuads(v bool) {
+	f.renderquads = v
+}
+
+func (f *FluidSimD) RenderQuads() bool {
+	return f.renderquads
+}
+
+func (f *FluidSimD) SetResolveCollisions(v bool) {
+	f.resolvecollisions = v
+}
+
+func (f *FluidSimD) ResolveCollisions() bool {
+	return f.resolvecollisions
+}
+
+func (f *FluidSimD) NextSolver() {
+	f.resolveridx = (f.resolveridx + 1) % len(f.resolvers)
+
+}
+
+func (f *FluidSimD) PreviousSolver() {
+	f.resolveridx = f.resolveridx - 1
+	if f.resolveridx < 0 {
+		f.resolveridx = len(f.resolvers) - 1
+	}
+}

elements/mappedentity.go

@@ -0,0 +1,18 @@
+package elements
+
+import (
+	"fluids/gamedata"
+
+	"github.com/hajimehoshi/ebiten/v2"
+)
+
+type MappedEntity interface {
+	Draw()
+	Update()
+	GetSprite() *ebiten.Image
+	GetDimensions() gamedata.Vector
+	GetPosition() gamedata.Vector
+	SetPosition(gamedata.Vector)
+	SetPaused(bool)
+	Paused() bool
+}

elements/mappedentitybase.go

@@ -0,0 +1,51 @@
+package elements
+
+import (
+	"fluids/gamedata"
+
+	"github.com/hajimehoshi/ebiten/v2"
+)
+
+type MappedEntityBase struct {
+	Sprite     *ebiten.Image
+	dimensions gamedata.Vector
+	position   gamedata.Vector
+	paused     bool
+}
+
+func NewMappedEntityBase(dimensions gamedata.Vector) *MappedEntityBase {
+	meb := &MappedEntityBase{
+		dimensions: dimensions,
+	}
+	return meb
+}
+
+func (m *MappedEntityBase) Draw() {
+}
+
+func (m *MappedEntityBase) Update() {
+}
+
+func (m *MappedEntityBase) GetSprite() *ebiten.Image {
+	return m.Sprite
+}
+
+func (m *MappedEntityBase) GetDimensions() gamedata.Vector {
+	return m.dimensions
+}
+
+func (m *MappedEntityBase) GetPosition() gamedata.Vector {
+	return m.position
+}
+
+func (m *MappedEntityBase) SetPosition(pos gamedata.Vector) {
+	m.position = pos
+}
+
+func (m *MappedEntityBase) SetPaused(p bool) {
+	m.paused = p
+}
+
+func (m *MappedEntityBase) Paused() bool {
+	return m.paused
+}

fluid/boundary.go

@@ -0,0 +1,19 @@
+package fluid
+
+type BoundaryDimensions struct {
+	X int
+	Y int
+}
+
+type Boundary struct {
+	Dimensions BoundaryDimensions
+	Cells      []bool
+}
+
+func NewBoundary(dimensions BoundaryDimensions) *Boundary {
+	b := &Boundary{
+		Dimensions: dimensions,
+		Cells:      make([]bool, dimensions.X*dimensions.Y),
+	}
+	return b
+}

fluid/field.go

@@ -0,0 +1,74 @@
+package fluid
+
+import "math"
+
+type CellType int
+
+const (
+	CellTypeFluid = iota
+	CellTypeAir
+	CellTypeSolid
+	CellTypeMax
+)
+
+type VelocityField struct {
+	Dimensions   FieldVector //in meters
+	Numcells     int
+	Nx, Ny       int     //number of cells in x, y
+	H            float32 //field spacing, in meters
+	InvH         float32
+	U, V         []float32 //field components < u(x,y), v(x,y) >
+	PrevU, PrevV []float32 //previous field components < u(x,y), v(x,y) >
+	DU, DV       []float32 //weighted sums for components < u(x,y), v(x,y) >
+	S            []float32 //fluid cell scales
+	CellType     []float32 //type of cell
+}
+
+func NewVelocityField(dimensions FieldVector, spacing float32) *VelocityField {
+	vf := &VelocityField{
+		Dimensions: dimensions,
+		H:          spacing,
+		InvH:       1 / spacing,
+	}
+
+	vf.Initialize()
+	return vf
+}
+
+func (vf *VelocityField) Initialize() {
+
+	vf.Nx = int(math.Floor(float64(vf.Dimensions.X/vf.H))) + 1
+	vf.Ny = int(math.Floor(float64(vf.Dimensions.Y/vf.H))) + 1
+
+	cellcount := vf.Nx * vf.Ny
+
+	vf.Numcells = cellcount
+
+	vf.U = make([]float32, cellcount)
+	vf.V = make([]float32, cellcount)
+	vf.PrevU = make([]float32, cellcount)
+	vf.PrevV = make([]float32, cellcount)
+	vf.DU = make([]float32, cellcount)
+	vf.DV = make([]float32, cellcount)
+	vf.S = make([]float32, cellcount)
+	vf.CellType = make([]float32, cellcount)
+
+	for i := 0; i < vf.Nx; i++ {
+		for j := 0; j < vf.Ny; j++ {
+			var stype float32 = CellTypeAir //port, code corrected
+			//var stype float32 = 1 //10mp value, claims fluid but isn't
+			//var stype float32 = CellTypeSolid //port, with 10mp intent
+			if i == 0 || i == vf.Nx-1 || j == 0 {
+				stype = CellTypeFluid //port, code corrected
+				//stype = 0 //10mp value, claims solid but isn't
+				//stype = CellTypeFluid //port, with 10mp intent
+			}
+			vf.S[i*vf.Ny+j] = stype
+		}
+	}
+
+}
+
+func (vf *VelocityField) GetIndex(i, j int) int {
+	return i*vf.Ny + j
+}

fluid/fieldvector.go

@@ -0,0 +1,36 @@
+package fluid
+
+import "math"
+
+type FieldVector struct {
+	X float32
+	Y float32
+}
+
+func (v FieldVector) DotProduct(p FieldVector) float32 {
+	return v.X*p.X + v.Y*p.Y
+}
+
+func (v FieldVector) Magnitude() float32 {
+	return float32(math.Sqrt(float64(v.X*v.X + v.Y*v.Y)))
+}
+
+func (v FieldVector) Add(p FieldVector) FieldVector {
+	return FieldVector{X: v.X + p.X, Y: v.Y + p.Y}
+}
+
+func (v FieldVector) Subtract(p FieldVector) FieldVector {
+	return FieldVector{X: v.X - p.X, Y: v.Y - p.Y}
+}
+
+func (v FieldVector) Scale(s float32) FieldVector {
+	return FieldVector{X: v.X * s, Y: v.Y * s}
+}
+
+func (v FieldVector) Normalize() FieldVector {
+	mag := v.Magnitude()
+	return FieldVector{
+		X: v.X / mag,
+		Y: v.Y / mag,
+	}
+}

fluid/fluid.go

@@ -0,0 +1,713 @@
+// Based on code by Matthias Müller (Ten Minute Physics), MIT License
+// See licenses/TenMinutePhysics for full license text
+
+package fluid
+
+import (
+	"fluids/edt"
+	"fluids/utils"
+	"math"
+)
+
+const (
+	FluidDefaultDensity               = 1000 //kilogram per cubic meter
+	FluidDefaultGravity               = -9.8
+	FluidDefaultTimeStep              = 1. / 120
+	FluidDefaultParticleRadiusPercent = 0.3
+	FluidDefaultAngle                 = 0
+	FluidIncompressibilityIterations  = 100
+	FluidDefaultRestDensity           = 0
+	FluidDefaultFlipPicRatio          = 0.9
+	FluidDefaultOverRelaxation        = 1.9
+	FluidDefaultSubSteps              = 2
+)
+
+type Fluid struct {
+	Particles           []Particle
+	Field               *VelocityField
+	particleRadius      float32
+	numCellParticles    []int     //SPATIAL HASHING: count of particles per cell
+	firstCellParticle   []int     //SPATIAL HASHING: index of first particle for given cell
+	cellParticleIds     []int     //SPATIAL HASHING: id of all particles
+	pNumX, pNumY        int       //SPATIAL HASHING: number of particle cells x,y
+	pNumCells           int       //SPATIAL HASHING: number of particle cells in total
+	pInvSpacing         float32   //SPATIAL HASHING: one over the particle spacing
+	ParticleDensity     []float32 //density of particles within each field cell
+	ParticleRestDensity float32
+	fluidDensity        float32
+	angle               float32
+	dt                  float32     //delta time step
+	gravitynormal       FieldVector //unit vector for gravity
+	stepacceleration    FieldVector //step acceleration due to gravity
+	flipPicRatio        float32
+	numSubSteps         int  //number of simulation substeps
+	Block               bool //rectangular or circular field container
+
+	//for arbitrary boundaries
+	edt *edt.EDT
+}
+
+func NewFluid(dimensions FieldVector, spacing float32) *Fluid {
+	f := &Fluid{
+		Field:               NewVelocityField(dimensions, spacing),
+		dt:                  FluidDefaultTimeStep,
+		fluidDensity:        FluidDefaultDensity,
+		ParticleRestDensity: FluidDefaultRestDensity,
+		flipPicRatio:        FluidDefaultFlipPicRatio,
+		numSubSteps:         FluidDefaultSubSteps,
+		Block:               true,
+		edt:                 nil,
+	}
+
+	f.Initialize()
+	return f
+}
+
+func (f *Fluid) Initialize() {
+	f.InitializeParticles()
+
+	f.SetAngle(FluidDefaultAngle)
+	f.ParticleDensity = make([]float32, f.Field.Nx*f.Field.Ny)
+}
+
+func (f *Fluid) InitializeParticles() {
+	//compute particle radius based on our field spacing and radius scaling
+	f.particleRadius = FluidDefaultParticleRadiusPercent * f.Field.H
+
+	//we want to prepare a hexagonally packed grid of particles within our fluid field
+	dx := float32(2.0 * f.particleRadius)
+	dy := float32(math.Sqrt(3.) / 2. * float64(dx))
+
+	//number of particles that fit within our specified dimensions
+	numPx := int(math.Floor(float64(f.Field.Dimensions.X-2*f.Field.H-2*f.particleRadius) / float64(dx)))
+	numPy := int(math.Floor(float64(f.Field.Dimensions.Y/2-2*f.Field.H-2*f.particleRadius) / float64(dy)))
+
+	particlecount := numPx * numPy
+	f.Particles = make([]Particle, particlecount)
+
+	//now prepare our hashing structures, purely for collision adjustment
+	//first let's compute some boundaries
+	f.pInvSpacing = 1 / (2.2 * f.particleRadius)
+	f.pNumX = int(math.Floor(float64(f.Field.Dimensions.X*f.pInvSpacing))) + 1
+	f.pNumY = int(math.Floor(float64(f.Field.Dimensions.Y*f.pInvSpacing))) + 1
+	f.pNumCells = f.pNumX * f.pNumY
+
+	//now use boundaries to construct our spatial hashing containers
+	f.numCellParticles = make([]int, f.pNumCells)
+	f.firstCellParticle = make([]int, f.pNumCells+1)
+	f.cellParticleIds = make([]int, int(particlecount))
+
+	//create particles
+	var p int = 0
+	for i := 0; i < numPx; i++ {
+		for j := 0; j < numPy; j++ {
+
+			var rowAdjust float32 = 0
+			if j%2 == 1 {
+				rowAdjust = f.particleRadius
+			}
+
+			f.Particles[p].Position.X = f.Field.H + f.particleRadius + dx*float32(i) + rowAdjust
+			f.Particles[p].Position.Y = f.Field.H + f.particleRadius + dy*float32(j)
+
+			p++
+		}
+	}
+}
+
+func (f *Fluid) Step() {
+
+	for range f.numSubSteps {
+		f.IntegrateParticles()
+		f.SeparateParticles()
+		f.HandleParticleCollisions()
+		f.TransferVelocities(true)
+		f.UpdateParticleDensity()
+		f.SolveIncompressibility()
+		f.TransferVelocities(false)
+	}
+	//f.UpdateParticleColors()
+}
+
+func (f *Fluid) RecalculateStepAcceleration() {
+	scale := FluidDefaultGravity * f.dt
+	f.stepacceleration = FieldVector{
+		X: scale * f.gravitynormal.X,
+		Y: scale * f.gravitynormal.Y,
+	}
+}
+
+// change the time step
+func (f *Fluid) SetTimeStep(dt float32) {
+	f.dt = dt
+	f.RecalculateStepAcceleration()
+}
+
+// set new fluid angle, recompute gravity normal, initiate step accel recomputation
+func (f *Fluid) SetAngle(angle float32) {
+	f.angle = angle
+	sin, cos := math.Sincos(float64(angle))
+	f.gravitynormal = FieldVector{
+		X: -float32(sin),
+		Y: float32(cos),
+	}
+	f.RecalculateStepAcceleration()
+}
+
+func (f *Fluid) IntegrateParticles() {
+	ax, ay := f.stepacceleration.X, f.stepacceleration.Y
+	dt := f.dt
+	for i := range f.Particles {
+		p := &f.Particles[i]
+		p.Velocity.X += ax
+		p.Velocity.Y += ay
+		p.Position.X += p.Velocity.X * dt
+		p.Position.Y += p.Velocity.Y * dt
+	}
+}
+
+func (f *Fluid) PerformParticleSpatialHash() {
+
+	clear(f.numCellParticles)
+
+	//find spatcial hash cell for this particle's position, increment count there
+	for i := range f.Particles {
+		p := &f.Particles[i]
+		xi := utils.Clamp(int(p.Position.X*f.pInvSpacing), 0, f.pNumX-1)
+		yi := utils.Clamp(int(p.Position.Y*f.pInvSpacing), 0, f.pNumY-1)
+		cellnumber := xi*f.pNumY + yi
+		f.numCellParticles[cellnumber]++
+	}
+
+	//build up firstCellParticle bucket
+	var sum int = 0
+	for i := range f.pNumCells {
+		sum += f.numCellParticles[i]
+		f.firstCellParticle[i] = sum
+	}
+	f.firstCellParticle[f.pNumCells] = sum
+
+	//next work backwards on it to assign particle ids, adjusting the
+	//firstCellParticle bucket as we go
+	for i := range f.Particles {
+		p := &f.Particles[i]
+		xi := utils.Clamp(int(p.Position.X*f.pInvSpacing), 0, f.pNumX-1)
+		yi := utils.Clamp(int(p.Position.Y*f.pInvSpacing), 0, f.pNumY-1)
+		cellnumber := xi*f.pNumY + yi
+		f.firstCellParticle[cellnumber]--
+		f.cellParticleIds[f.firstCellParticle[cellnumber]] = i
+	}
+}
+
+func (f *Fluid) SeparateParticles() {
+
+	//setup
+	minDist := 2 * f.particleRadius
+	minDist2 := minDist * minDist
+
+	//build our spatial hashing table
+	f.PerformParticleSpatialHash()
+
+	//perform broad phase search
+	for i := range f.Particles {
+
+		p := &f.Particles[i]
+
+		pxi := int(p.Position.X * f.pInvSpacing)
+		pyi := int(p.Position.Y * f.pInvSpacing)
+		x0 := max(pxi-1, 0)
+		y0 := max(pyi-1, 0)
+		x1 := min(pxi+1, f.pNumX)
+		y1 := min(pyi+1, f.pNumY)
+
+		for xi := x0; xi < x1; xi++ {
+			for yi := y0; yi < y1; yi++ {
+				//find current cell
+				cell := xi*f.pNumY + yi
+
+				//set start and ending particles to check between this and next cell
+				first := f.firstCellParticle[cell]
+				last := f.firstCellParticle[cell+1]
+
+				//iterate over these candidates
+				for j := first; j < last; j++ {
+					//skip self
+					id := f.cellParticleIds[j]
+					if id == i {
+						continue
+					}
+
+					//compute delta
+					dx := f.Particles[id].Position.X - p.Position.X
+					dy := f.Particles[id].Position.Y - p.Position.Y
+					dist2 := dx*dx + dy*dy
+
+					//rule out non-collisions, ignore overlaps (we screwed)
+					if (dist2 == 0) || (dist2 > minDist2) {
+						continue
+					}
+
+					//need to use more expensive sqrt now to find our penetration
+					//then compute a scalar for vector normalization, applied
+					//equally in half to both particles
+					dist := float32(math.Sqrt(float64(dist2)))
+					penetration := (minDist - dist)
+					scale := penetration / dist
+					dx *= scale / 2
+					dy *= scale / 2
+
+					p.Position.X -= dx
+					p.Position.Y -= dy
+					f.Particles[id].Position.X += dx
+					f.Particles[id].Position.Y += dy
+				}
+			}
+		}
+	}
+}
+
+func (f *Fluid) HandleParticleCollisions() {
+	//setup
+	/*
+		minDist := 2 * f.particleRadius
+		minDist2 := minDist * minDist
+	*/
+
+	if f.edt != nil {
+		f.HandleBoundaryCollisions()
+		return
+	}
+
+	if f.Block {
+		minX := f.Field.H + f.particleRadius
+		maxX := float32(f.Field.Nx-1)*f.Field.H - f.particleRadius
+		minY := f.Field.H + f.particleRadius
+		maxY := float32(f.Field.Ny-1)*f.Field.H - f.particleRadius
+
+		for i := range f.Particles {
+			p := &f.Particles[i]
+
+			if p.Position.X < minX {
+				p.Position.X = minX
+				p.Velocity.X = 0
+			}
+
+			if p.Position.X > maxX {
+				p.Position.X = maxX
+				p.Velocity.X = 0
+			}
+
+			if p.Position.Y < minY {
+				p.Position.Y = minY
+				p.Velocity.Y = 0
+			}
+
+			if p.Position.Y > maxY {
+				p.Position.Y = maxY
+				p.Velocity.Y = 0
+			}
+
+		}
+	} else {
+		//find fluid cell in which the particle resides
+		//compute distance of this cell from 'center' cell of our circle
+		//if distance exceeds radius, clamp particle velocity and position
+		//to edge cell
+
+		//find origin of circle, conveniently these dimensions are also our radius
+		originX := float32(f.Field.Nx) * f.Field.H / 2
+		originY := float32(f.Field.Ny) * f.Field.H / 2
+		radius := originX - f.Field.H - f.particleRadius
+
+		//cxi := utils.Clamp(int(math.Floor(float64(originX*f.Field.InvH))), 0, f.Field.Nx-1)
+		//cyi := utils.Clamp(int(math.Floor(float64(originY*f.Field.InvH))), 0, f.Field.Ny-1)
+		//originCellIdx := cxi*f.Field.Ny + cyi
+
+		//find fluid cell for particle
+		for i := range f.Particles {
+			p := &f.Particles[i]
+
+			//xi := utils.Clamp(int(math.Floor(float64(p.Position.X*f.Field.InvH))), 0, f.Field.Nx-1)
+			//yi := utils.Clamp(int(math.Floor(float64(p.Position.Y*f.Field.InvH))), 0, f.Field.Ny-1)
+			//
+			//dx := originX - px
+			//dy := originY - py
+
+			dx := p.Position.X - originX
+			dy := p.Position.Y - originY
+			dist2 := dx*dx + dy*dy
+			if dist2 < radius*radius || dist2 == 0 {
+				continue
+			}
+
+			dist := float32(math.Sqrt(float64(dist2)))
+
+			newx := originX + dx*radius/dist
+			newy := originY + dy*radius/dist
+
+			p.Position.X = newx
+			p.Position.Y = newy
+			p.Velocity.X = 0
+			p.Velocity.Y = 0
+
+		}
+
+	}
+}
+
+func (f *Fluid) TransferVelocities(togrid bool) {
+
+	//setup
+	n := f.Field.Ny
+	h1 := f.Field.InvH
+	h2 := 0.5 * f.Field.H
+
+	if togrid {
+		copy(f.Field.PrevU, f.Field.U)
+		copy(f.Field.PrevV, f.Field.V)
+		clear(f.Field.DU)
+		clear(f.Field.DV)
+		clear(f.Field.U)
+		clear(f.Field.V)
+
+		for i := 0; i < f.Field.Numcells; i++ {
+			if f.Field.S[i] == 0 {
+				f.Field.CellType[i] = CellTypeSolid
+			} else {
+				f.Field.CellType[i] = CellTypeAir
+			}
+		}
+
+		for i := range f.Particles {
+			x := f.Particles[i].Position.X
+			y := f.Particles[i].Position.Y
+			xi := utils.Clamp(int(math.Floor(float64(x*h1))), 0, f.Field.Nx-1)
+			yi := utils.Clamp(int(math.Floor(float64(y*h1))), 0, f.Field.Ny-1)
+			cellnumber := xi*n + yi
+			if f.Field.CellType[cellnumber] == CellTypeAir {
+				f.Field.CellType[cellnumber] = CellTypeFluid
+			}
+		}
+	}
+
+	for component := 0; component < 2; component++ {
+		var dx, dy float32
+		var fl, prevF, d *[]float32
+
+		if component == 0 {
+			dx = 0
+			dy = h2
+			fl = &f.Field.U
+			prevF = &f.Field.PrevU
+			d = &f.Field.DU
+		} else {
+			dx = h2
+			dy = 0
+			fl = &f.Field.V
+			prevF = &f.Field.PrevV
+			d = &f.Field.DV
+		}
+
+		for i := range f.Particles {
+			p := &f.Particles[i]
+			x := p.Position.X
+			y := p.Position.Y
+
+			x = utils.Clamp32(x, f.Field.H, float32((f.Field.Nx-1))*f.Field.H)
+			y = utils.Clamp32(y, f.Field.H, float32((f.Field.Ny-1))*f.Field.H)
+
+			//find cell components neighbouring x,y
+			x0 := min(int(math.Floor(float64((x-dx)*h1))), f.Field.Nx-2)
+			x1 := min(x0+1, f.Field.Nx-2)
+			y0 := min(int(math.Floor(float64((y-dy)*h1))), f.Field.Ny-2)
+			y1 := min(y0+1, f.Field.Ny-2)
+
+			//compute scales
+			tx := ((x - dx) - float32(x0)*f.Field.H) * h1
+			ty := ((y - dy) - float32(y0)*f.Field.H) * h1
+			sx := 1.0 - tx
+			sy := 1.0 - ty
+
+			//compute weights
+			d0 := sx * sy
+			d1 := tx * sy
+			d2 := tx * ty
+			d3 := sx * ty
+
+			//cell indexes
+			nr0 := x0*n + y0
+			nr1 := x1*n + y0
+			nr2 := x1*n + y1
+			nr3 := x0*n + y1
+
+			if togrid {
+				//extract velocity for this component
+				var pv float32 // = this.particleVel[2 * i + component];
+				if component == 0 {
+					pv = p.Velocity.X
+				} else {
+					pv = p.Velocity.Y
+				}
+
+				//to our velocity component we will add the contribution of
+				//the particle velocity scaled according to its weight
+				(*fl)[nr0] += pv * d0
+				(*fl)[nr1] += pv * d1
+				(*fl)[nr2] += pv * d2
+				(*fl)[nr3] += pv * d3
+
+				//we also keep a running total of all the weights (d is r in the notes)
+				(*d)[nr0] += d0
+				(*d)[nr1] += d1
+				(*d)[nr2] += d2
+				(*d)[nr3] += d3
+			} else {
+
+				var offset int
+				var v float32
+				if component == 0 {
+					offset = n
+					v = p.Velocity.X
+				} else {
+					offset = 1
+					v = p.Velocity.Y
+				}
+
+				var valid0, valid1, valid2, valid3 float32 = 0., 0., 0., 0.
+
+				if f.Field.CellType[nr0] != CellTypeAir || f.Field.CellType[nr0-offset] != CellTypeAir {
+					valid0 = 1
+				}
+
+				if f.Field.CellType[nr1] != CellTypeAir || f.Field.CellType[nr1-offset] != CellTypeAir {
+					valid1 = 1
+				}
+
+				if f.Field.CellType[nr2] != CellTypeAir || f.Field.CellType[nr2-offset] != CellTypeAir {
+					valid2 = 1
+				}
+
+				if f.Field.CellType[nr3] != CellTypeAir || f.Field.CellType[nr3-offset] != CellTypeAir {
+					valid3 = 1
+				}
+
+				//weights
+				wsum := valid0*d0 + valid1*d1 + valid2*d2 + valid3*d3
+
+				if wsum > 0.0 {
+					picV := (valid0*d0*(*fl)[nr0] + valid1*d1*(*fl)[nr1] + valid2*d2*(*fl)[nr2] + valid3*d3*(*fl)[nr3]) / wsum
+					corr := (valid0*d0*((*fl)[nr0]-(*prevF)[nr0]) + valid1*d1*((*fl)[nr1]-(*prevF)[nr1]) +
+						valid2*d2*((*fl)[nr2]-(*prevF)[nr2]) + valid3*d3*((*fl)[nr3]-(*prevF)[nr3])) / wsum
+					flipV := v + corr
+					vNew := (1-f.flipPicRatio)*picV + f.flipPicRatio*flipV
+
+					if component == 0 {
+						p.Velocity.X = vNew
+					} else {
+						p.Velocity.Y = vNew
+					}
+				}
+			}
+		}
+
+		if togrid {
+			//finally, for all q's, we divide by the weighted sum of that cell
+			for i := 0; i < len(*fl); i++ {
+				if (*d)[i] > 0 {
+					(*fl)[i] = (*fl)[i] / (*d)[i]
+				}
+			}
+
+			//next we restore cells are solid to their previous values
+			for i := 0; i < f.Field.Nx; i++ {
+				for j := 0; j < f.Field.Ny; j++ {
+					idx := i*n + j
+					issolid := f.Field.CellType[idx] == CellTypeSolid
+
+					if issolid || (i > 0 && f.Field.CellType[(i-1)*n+j] == CellTypeSolid) {
+						f.Field.U[idx] = f.Field.PrevU[idx]
+					}
+
+					if issolid || (j > 0 && f.Field.CellType[i*n+j-1] == CellTypeSolid) {
+						f.Field.V[idx] = f.Field.PrevV[idx]
+					}
+
+				}
+			}
+		}
+
+	}
+
+}
+
+func (f *Fluid) UpdateParticleDensity() {
+	n := f.Field.Ny
+	h := f.Field.H
+	h1 := f.Field.InvH
+	h2 := 0.5 * h
+
+	clear(f.ParticleDensity)
+
+	for i := range f.Particles {
+		p := &f.Particles[i]
+		x := p.Position.X
+		y := p.Position.Y
+
+		x = utils.Clamp32(x, h, float32(f.Field.Nx-1)*h)
+		y = utils.Clamp32(y, h, float32(f.Field.Ny-1)*h)
+
+		x0 := int(math.Floor(float64((x - h2) * h1)))
+		x1 := min(x0+1, f.Field.Nx-2)
+
+		y0 := int(math.Floor(float64((y - h2) * h1)))
+		y1 := min(y0+1, f.Field.Ny-2)
+
+		tx := ((x - h2) - float32(x0)*h) * h1
+		ty := ((y - h2) - float32(y0)*h) * h1
+
+		sx := 1.0 - tx
+		sy := 1.0 - ty
+
+		if x0 < f.Field.Nx && y0 < f.Field.Ny {
+			f.ParticleDensity[x0*n+y0] += sx * sy
+		}
+		if x1 < f.Field.Nx && y0 < f.Field.Ny {
+			f.ParticleDensity[x1*n+y0] += tx * sy
+		}
+		if x1 < f.Field.Nx && y1 < f.Field.Ny {
+			f.ParticleDensity[x1*n+y1] += tx * ty
+		}
+		if x0 < f.Field.Nx && y1 < f.Field.Ny {
+			f.ParticleDensity[x0*n+y1] += sx * ty
+		}
+	}
+
+	if f.ParticleRestDensity == 0.0 {
+		var sum float32 = 0.0
+		numFluidCells := 0
+
+		for i := 0; i < f.Field.Nx*f.Field.Ny; i++ {
+			if f.Field.CellType[i] == CellTypeFluid {
+				sum += f.ParticleDensity[i]
+				numFluidCells++
+			}
+		}
+
+		if numFluidCells > 0 {
+			f.ParticleRestDensity = sum / float32(numFluidCells)
+		}
+	}
+}
+
+func (f *Fluid) SolveIncompressibility() {
+
+	//setup
+	copy(f.Field.PrevU, f.Field.U)
+	copy(f.Field.PrevV, f.Field.V)
+
+	n := f.Field.Ny
+	//cp := f.fluidDensity * f.Field.H / FluidDefaultTimeStep
+
+	for iter := 0; iter < FluidIncompressibilityIterations; iter++ {
+		for i := 1; i < f.Field.Nx-1; i++ {
+			for j := 1; j < f.Field.Ny-1; j++ {
+
+				if f.Field.CellType[i*n+j] != CellTypeFluid {
+					continue
+				}
+
+				//compute indexes of surrounding cells
+				center := i*n + j
+				left := (i-1)*n + j
+				right := (i+1)*n + j
+				bottom := i*n + j - 1
+				top := i*n + j + 1
+
+				sx0 := f.Field.S[left]
+				sx1 := f.Field.S[right]
+				sy0 := f.Field.S[bottom]
+				sy1 := f.Field.S[top]
+				s := sx0 + sx1 + sy0 + sy1
+				if s == 0 {
+					continue
+				}
+
+				//divergence -> we want this to be zero
+				div := f.Field.U[right] - f.Field.U[center] + f.Field.V[top] - f.Field.V[center]
+
+				if f.ParticleRestDensity > 0.0 {
+					var k float32 = 1.0
+					var compression float32 = f.ParticleDensity[i*n+j] - f.ParticleRestDensity
+					if compression > 0.0 {
+						div = div - k*compression
+					}
+				}
+
+				p := -div / s
+				p *= FluidDefaultOverRelaxation //overRelaxation
+				//f.Field.Pressure[center] += cp * p
+
+				f.Field.U[center] -= sx0 * p
+				f.Field.U[right] += sx1 * p
+				f.Field.V[center] -= sy0 * p
+				f.Field.V[top] += sy1 * p
+
+			}
+		}
+	}
+
+}
+
+func (f *Fluid) ToggleShape() {
+	f.Block = !f.Block
+}
+
+func (f *Fluid) SetBoundary(b *Boundary) {
+
+	if b != nil {
+		dim := edt.Bounds{
+			X: b.Dimensions.X,
+			Y: b.Dimensions.Y,
+		}
+
+		f.edt = edt.NewEDT(dim)
+		f.edt.AssignOccupancy(b.Cells)
+		f.edt.ComputeDistanceTransform()
+	} else {
+		f.edt = nil
+	}
+}
+
+// we perform our boundary resolution on particles according to our defined fluid boundary
+func (f *Fluid) HandleBoundaryCollisions() {
+
+	//grid spacing within our distance field
+	dx := f.Field.Dimensions.X / float32(f.edt.Dimensions.X)
+	dy := f.Field.Dimensions.Y / float32(f.edt.Dimensions.Y)
+
+	//find cell of distance field for which particles belong, check if it's in bounds
+	for i := range f.Particles {
+		p := &f.Particles[i]
+
+		xi := utils.Clamp(int(math.Floor(float64(p.Position.X/dx))), 0, f.edt.Dimensions.X-1)
+		yi := utils.Clamp(int(math.Floor(float64(p.Position.Y/dy))), 0, f.edt.Dimensions.Y-1)
+
+		cellidx := xi + yi*f.edt.Dimensions.X
+
+		//if our current cell isn't a boundary, then we skip
+		if f.edt.D[cellidx] == 0 {
+			continue
+		}
+
+		//find where the particle actually belongs
+		pos := f.edt.L[cellidx]
+		newx := (float32(pos.X-xi) + f.particleRadius) * dx
+		newy := (float32(pos.Y-yi) + f.particleRadius) * dy
+		p.Position.X += newx
+		p.Position.Y += newy
+		p.Velocity.X = 0
+		p.Velocity.Y = 0
+
+	}
+
+}

fluid/particle.go

@@ -0,0 +1,6 @@
+package fluid
+
+type Particle struct {
+	Position FieldVector
+	Velocity FieldVector
+}

fonts/fonts.go

@@ -0,0 +1,25 @@
+package fonts
+
+import (
+	"bytes"
+	_ "embed"
+	"log"
+
+	"github.com/hajimehoshi/ebiten/v2/text/v2"
+)
+
+const ()
+
+var (
+	PixelFont *text.GoTextFaceSource
+	//go:embed Rockboxcond12.ttf
+	pixel_fnt []byte
+)
+
+func init() {
+	s, err := text.NewGoTextFaceSource(bytes.NewReader(pixel_fnt))
+	if err != nil {
+		log.Fatal(err)
+	}
+	PixelFont = s
+}

game/game.go

@@ -3,92 +3,78 @@ package game
 import (
 	"fluids/elements"
 	"fluids/gamedata"
-	"fluids/quadtree"
 	"fmt"
 	"image/color"
 	"math"
-	"math/rand"
 
 	"github.com/hajimehoshi/ebiten/v2"
+	"github.com/hajimehoshi/ebiten/v2/ebitenutil"
 	"github.com/hajimehoshi/ebiten/v2/inpututil"
-	"github.com/hajimehoshi/ebiten/v2/vector"
 )
 
 const (
 	GameWidth  = 640
 	GameHeight = 360
-	GameParticleCount  = 2000
-	GameGravity        = 2
-	GameParticleRadius = 5
-	GameDamping        = .7
-	GameDeltaTimeStep  = 0.5
-
-	GameInfluenceRadius = 30
+	GameFSDW   = 200
+	GameFSDH   = 100
+	GameSims   = 3 //number of total sims
 )
 
 type Game struct {
-	particles         []*elements.Particle
+	//control data
 	cycle       int
-	particlebox       *gamedata.Vector
-	particlebuff      *ebiten.Image
-	quadtree          *quadtree.Quadtree
-	collisionquad     quadtree.Quadrant
 	paused      bool
-	renderquads       bool
-	resolvecollisions bool
-	resolvers         []func(particle *elements.Particle)
-	resolveridx       int
+	fluidsimidx int
+
+	//key elements
+	fluidsimd  *elements.FluidSimD
+	fluidsim10 []*elements.FluidSim10
+	alertbox   *elements.Alert
+	flask      *elements.RoundedBottomFlask
+
+	//cache elements
+	fluidsim10width  float64
+	fluidsim10height float64
+
+	//other
+	fluidsim10angle float64 //purely for debugging
+	mdown           bool
+	mdx, mdy        int
+	mdangle         float64 //angle at mousedown
+	mfangle         float64 //last flask angle
+
+	//colors
+	flaskcolor color.RGBA
 }
 
 func NewGame() *Game {
 	g := &Game{
-		particlebuff:      ebiten.NewImage(GameParticleRadius*2, GameParticleRadius*2),
 		paused:      false,
-		renderquads:       false,
-		resolvecollisions: false,
-		resolveridx:       0,
+		fluidsimidx: 0,
+		alertbox:    elements.NewAlert(),
+		fluidsimd:   elements.NewFluidSimD(),
+		flask:       elements.NewRoundedBottomFlask(),
 	}
 
-	g.particlebox = &gamedata.Vector{
-		X: GameWidth - 50,
-		Y: GameHeight - 50,
-	}
-
-	quad := quadtree.Quadrant{
-		Position:   gamedata.Vector{X: GameWidth / 2, Y: GameHeight / 2},
-		Dimensions: gamedata.Vector{X: GameWidth, Y: GameHeight},
-	}
-	g.quadtree = quadtree.New(quad)
-
-	vector.FillCircle(g.particlebuff, GameParticleRadius, GameParticleRadius, GameParticleRadius, color.White, true)
-
-	g.collisionquad = quadtree.Quadrant{
-		Position: gamedata.Vector{
-			X: 0,
-			Y: 0,
-		},
-		Dimensions: gamedata.Vector{
-			X: GameParticleRadius * 2,
-			Y: GameParticleRadius * 2,
-		},
-	}
-
-	g.resolvers = append(g.resolvers, g.ResolveCollisionsA)
-	g.resolvers = append(g.resolvers, g.ResolveCollisionsB)
-
-	//g.InitializeColliders()
-	g.InitializeParticles()
-
+	g.Initialize()
 	return g
 }
 
 func (g *Game) Update() error {
 
 	g.ParseInputs()
-	g.RebuildQuadtree()
 
 	if !g.paused {
-		g.UpdateParticles()
+		switch g.fluidsimidx {
+		case 0:
+			g.UpdateFlask()
+		case 1:
+			g.fluidsimd.Update()
+		case 2:
+			g.UpdateFluidsim10()
+		default:
+			break
+		}
 	}
 
 	g.cycle++
@@ -96,315 +82,282 @@ func (g *Game) Update() error {
 	return nil
 }
 
+func (g *Game) UpdateFluidsim10() {
+	for _, sim := range g.fluidsim10 {
+		sim.Update()
+	}
+}
+
 func (g *Game) Draw(screen *ebiten.Image) {
 	screen.Clear()
-	g.RenderParticles(screen)
-	g.RenderBox(screen)
 
-	if g.renderquads {
-		g.RenderQuadrants(screen)
+	switch g.fluidsimidx {
+	case 0:
+		g.RenderFlask(screen)
+	case 1:
+		g.RenderFluidSimD(screen)
+	case 2:
+		g.RenderFluidSim10(screen)
+
+	default:
+		break
 	}
+
+	if g.paused {
+		op := &ebiten.DrawImageOptions{}
+		op.GeoM.Translate(50, 50)
+		screen.DrawImage(g.alertbox.Sprite, op)
+	}
+}
+
+func (g *Game) RenderFluidSimD(img *ebiten.Image) {
+	g.fluidsimd.Draw()
+	pos := g.fluidsimd.GetPosition()
+	op := &ebiten.DrawImageOptions{}
+	op.GeoM.Translate(pos.X, pos.Y)
+	img.DrawImage(g.fluidsimd.GetSprite(), op)
+}
+
+func (g *Game) RenderFluidSim10(img *ebiten.Image) {
+
+	for _, sim := range g.fluidsim10 {
+		sim.Draw()
+
+		angle := sim.GetAngle()
+		pos := sim.GetPosition()
+		op := &ebiten.DrawImageOptions{}
+		op.GeoM.Translate(-g.fluidsim10width/2, -g.fluidsim10height/2)
+		op.GeoM.Scale(1, -1)
+		op.GeoM.Rotate(angle)
+		//	op.GeoM.Translate(g.fluidsim10width/2, g.fluidsim10height/2)
+		op.GeoM.Translate(pos.X, pos.Y)
+		img.DrawImage(sim.GetSprite(), op)
+	}
+
+	//debug info
+	x, y := ebiten.CursorPosition()
+	deg := g.fluidsim10angle * 180 / math.Pi
+
+	str := fmt.Sprintf("Mouse (x: %d, y: %d) Origin (x: %d, y: %d) Angle (%f)", x, y, GameWidth/2, GameHeight/2, deg)
+	ebitenutil.DebugPrint(img, str)
 }
 
 func (g *Game) Layout(x, y int) (int, int) {
 	return GameWidth, GameHeight
 }
 
-func (g *Game) RenderParticles(img *ebiten.Image) {
-	// mx, my := ebiten.CursorPosition()
-	//clr := color.RGBA{R: 0xff, G: 0x00, B: 0x00, A: 0xff}
-
-	for _, particle := range g.particles {
-
-		x0 := particle.Position.X - GameParticleRadius
-		y0 := particle.Position.Y - GameParticleRadius
-
-		//redness := float32(particle.Position.Y / g.particlebox.Y)
-		//blueness := 1 - redness
-
-		op := &ebiten.DrawImageOptions{}
-		op.GeoM.Translate(x0, y0)
-		//op.ColorScale.Scale(redness, 0, blueness, 1)
-		img.DrawImage(g.particlebuff, op)
-
-		//vector.FillCircle(img, float32(particle.Position.X), float32(particle.Position.Y), float32(particle.Radius), color.White, true)
-		// vector.StrokeCircle(img, float32(particle.Position.X), float32(particle.Position.Y), GameInfluenceRadius, 1, clr, true)
-		// vector.StrokeLine(img, float32(mx), float32(my), float32(particle.Position.X), float32(particle.Position.Y), 2, clr, true)
-	}
-}
-
-func (g *Game) RenderBox(img *ebiten.Image) {
-	x0 := (GameWidth - g.particlebox.X) / 2
-	y0 := (GameHeight - g.particlebox.Y) / 2
-	vector.StrokeRect(img, float32(x0), float32(y0), float32(g.particlebox.X), float32(g.particlebox.Y), 2, color.White, true)
-}
-
-func (g *Game) InitializeColliders() {
-	/*
-		//box container
-		box := &colliders.BaseRectCollider{}
-		box.SetDimensions(gamedata.Vector{
-			X: GameWidth - 50,
-			Y: GameHeight - 50,
-		})
-		box.SetPosition(gamedata.Vector{
-			X: GameWidth / 2,
-			Y: GameHeight / 2,
-		})
-		box.SetContainer(true)
-		g.rectColliders = append(g.rectColliders, box)
-	*/
-}
-
-func (g *Game) InitializeParticles() {
-
-	g.particles = g.particles[:0]
-
-	xmin := (GameWidth-g.particlebox.X)/2 + GameParticleRadius
-	xmax := g.particlebox.X - GameParticleRadius*2
-	ymin := (GameHeight-g.particlebox.Y)/2 + GameParticleRadius
-	ymax := g.particlebox.Y - GameParticleRadius*2
-
-	for i := 0; i < GameParticleCount; i++ {
-
-		p := &elements.Particle{
-			Position: gamedata.Vector{
-				X: xmin + rand.Float64()*xmax,
-				Y: ymin + rand.Float64()*ymax,
-			},
-			Velocity: gamedata.Vector{
-				X: 0,
-				Y: 0,
-			},
-			Radius: GameParticleRadius,
-		}
-
-		g.particles = append(g.particles, p)
-	}
-
-}
-
-func (g *Game) UpdateParticles() {
-
-	dt := GameDeltaTimeStep
-
-	mx, my := ebiten.CursorPosition()
-	mpos := gamedata.Vector{X: float64(mx), Y: float64(my)}
-	maxdeflect := 40 * GameDeltaTimeStep * GameGravity
-
-	for _, particle := range g.particles {
-		particle.Velocity.Y += GameGravity * dt
-
-		//if inpututil.IsMouseButtonJustPressed(ebiten.MouseButtonLeft) {
-
-		if ebiten.IsMouseButtonPressed(ebiten.MouseButtonLeft) {
-			delta := gamedata.Vector{
-				X: mpos.X - particle.Position.X,
-				Y: mpos.Y - particle.Position.Y,
-			}
-
-			dist := math.Sqrt(delta.X*delta.X + delta.Y*delta.Y)
-			theta := math.Atan2(delta.Y, delta.X)
-
-			if dist < GameInfluenceRadius {
-
-				dx := dist * math.Cos(theta)
-				dy := dist * math.Sin(theta)
-
-				if dx != 0 {
-					gainx := (-1./GameInfluenceRadius)*math.Abs(dx) + 1.
-					particle.Velocity.X += 20 * gainx * -1 * math.Copysign(1, dx)
-				}
-
-				if dy != 0 {
-					gainy := (-1./GameInfluenceRadius)*math.Abs(dy) + 1.
-					particle.Velocity.Y += maxdeflect * gainy * -1 * math.Copysign(1, dy)
-				}
-			}
-
-		}
-
-		particle.Position.X += particle.Velocity.X * dt
-		particle.Position.Y += particle.Velocity.Y * dt
-
-		if g.resolvecollisions {
-			g.resolvers[g.resolveridx](particle)
-		}
-
-		g.BoundParticle(particle)
-	}
-
-}
-
-func (g *Game) BoundParticle(p *elements.Particle) {
-
-	xmin := (GameWidth-g.particlebox.X)/2 + p.Radius
-	xmax := xmin + g.particlebox.X - p.Radius*2
-
-	if p.Position.X > xmax {
-		p.Velocity.X *= -1 * GameDamping
-		p.Position.X = xmax
-	}
-
-	if p.Position.X < xmin {
-		p.Velocity.X *= -1 * GameDamping
-		p.Position.X = xmin
-	}
-
-	ymin := (GameHeight-g.particlebox.Y)/2 + p.Radius
-	ymax := ymin + g.particlebox.Y - p.Radius*2
-
-	if p.Position.Y > ymax {
-		p.Velocity.Y *= -1 * GameDamping
-		p.Position.Y = ymax
-	}
-
-	if p.Position.Y < ymin {
-		p.Velocity.Y *= -1 * GameDamping
-		p.Position.Y = ymin
-	}
-
-}
-
-func (g *Game) RenderQuadrants(img *ebiten.Image) {
-	clr := color.RGBA{R: 0xff, G: 0x00, B: 0x00, A: 0xff}
-
-	quadrants := g.quadtree.GetQuadrants()
-	for _, quad := range quadrants {
-		ox := float32(quad.Position.X - quad.Dimensions.X/2)
-		oy := float32(quad.Position.Y - quad.Dimensions.Y/2)
-		vector.StrokeRect(img, ox, oy, float32(quad.Dimensions.X), float32(quad.Dimensions.Y), 1, clr, true)
-	}
-
-}
-
 func (g *Game) ParseInputs() {
-	if inpututil.IsKeyJustPressed(ebiten.KeyR) {
-		g.InitializeParticles()
+
+	//simulation specific updates
+	switch g.fluidsimidx {
+	case 0:
+		g.ManageFlaskInputs()
+	case 1:
+		g.ManageFluidSimDInputs()
+	case 2:
+		g.ManageFluidSim10Inputs()
+
+	default:
+		break
 	}
 
+	//common updates
 	if inpututil.IsKeyJustPressed(ebiten.KeyP) {
 		g.paused = !g.paused
+		g.alertbox.SetText("PAUSED")
+		g.alertbox.Draw()
 	}
 
+	//swap fluid simulations
+	if inpututil.IsKeyJustPressed(ebiten.KeyPageUp) {
+		g.fluidsimidx = (g.fluidsimidx + 1) % GameSims
+	}
+
+	if inpututil.IsKeyJustPressed(ebiten.KeyPageDown) {
+		g.fluidsimidx = g.fluidsimidx - 1
+		if g.fluidsimidx < 0 {
+			g.fluidsimidx = GameSims - 1
+		}
+	}
+
+}
+
+func (g *Game) ManageFluidSimDInputs() {
+	//refresh particles
+	if inpututil.IsKeyJustPressed(ebiten.KeyR) {
+		g.fluidsimd.InitializeParticles()
+	}
+
+	//pause simulation
+	if inpututil.IsKeyJustPressed(ebiten.KeyP) {
+		g.fluidsimd.SetPaused(!g.fluidsimd.Paused())
+	}
+
+	//show quadtree quadrants
 	if inpututil.IsKeyJustPressed(ebiten.KeyQ) {
-		g.renderquads = !g.renderquads
+		g.fluidsimd.SetRenderQuads(!g.fluidsimd.RenderQuads())
 	}
 
+	//enable collision resolution
 	if inpututil.IsKeyJustPressed(ebiten.KeyC) {
-		g.resolvecollisions = !g.resolvecollisions
+		g.fluidsimd.SetResolveCollisions(!g.fluidsimd.ResolveCollisions())
 	}
 
+	//switch between collision resolvers
 	if inpututil.IsKeyJustPressed(ebiten.KeyLeft) {
-		g.resolveridx = g.resolveridx - 1
-		if g.resolveridx < 0 {
-			g.resolveridx = len(g.resolvers) - 1
-		}
+		g.fluidsimd.PreviousSolver()
 	}
 
 	if inpututil.IsKeyJustPressed(ebiten.KeyRight) {
-		g.resolveridx = (g.resolveridx + 1) % len(g.resolvers)
+		g.fluidsimd.NextSolver()
+	}
 }
 
+func (g *Game) ManageFluidSim10Inputs() {
+	//refresh particles
+	if inpututil.IsKeyJustPressed(ebiten.KeyR) {
+		for _, sim := range g.fluidsim10 {
+			sim.Initialize()
+			g.fluidsim10angle = 0
+		}
 	}
 
-func (g *Game) RebuildQuadtree() {
-	g.quadtree.Clear()
+	if inpututil.IsMouseButtonJustPressed(ebiten.MouseButtonLeft) {
+		g.mdown = true
+		g.mdx, g.mdy = ebiten.CursorPosition()
+	}
 
-	for _, p := range g.particles {
-		if !g.quadtree.Insert(p) {
-			fmt.Println("quadtree insertion failed")
+	if ebiten.IsMouseButtonPressed(ebiten.MouseButtonLeft) {
+		if g.mdown {
+			mx, my := ebiten.CursorPosition()
+
+			for _, sim := range g.fluidsim10 {
+				dx := float64(mx) - sim.GetPosition().X
+				dy := float64(my) - sim.GetPosition().Y
+				angle := math.Atan2(dy, dx)
+				g.fluidsim10angle = angle
+				sim.SetAngle(angle)
 			}
 		}
 	}
 
-func (g *Game) ResolveCollisionsA(particle *elements.Particle) {
-	//construct search quadrant from current particle
-	quadrant := quadtree.Quadrant{
-		Position:   particle.Position,
-		Dimensions: particle.GetDimensions(),
+	if inpututil.IsMouseButtonJustReleased(ebiten.MouseButtonLeft) {
+		g.mdown = false
 	}
 
-	//find list of possible maybe collisions, we inspect those in more detail
-	maybes := g.quadtree.FindAll(quadrant)
-
-	sqdist := float64(GameParticleRadius*GameParticleRadius) * 4
-
-	for _, p := range maybes {
-		if p == particle {
-			continue
-		}
-
-		pos := p.GetPosition()
-		delta := gamedata.Vector{
-			X: pos.X - particle.Position.X,
-			Y: pos.Y - particle.Position.Y,
-		}
-		dist2 := delta.X*delta.X + delta.Y*delta.Y
-
-		if dist2 == 0 {
-			// Same position: pick a fallback direction to avoid NaN
-			delta.X = 1
-			delta.Y = 0
-			dist2 = 1
-		}
-
-		if dist2 < sqdist {
-			d := math.Sqrt(dist2)
-			nx, ny := delta.X/d, delta.Y/d
-			overlap := GameParticleRadius*2 - d
-
-			pos.X += nx * overlap
-			pos.Y += ny * overlap
-			p.SetPosition(pos)
-			/*
-				newpos := gamedata.Vector{
-					X: pos.X + delta.X,
-					Y: pos.Y + delta.Y,
-				}
-				p.SetPosition(newpos)
-			*/
-		}
+	if inpututil.IsKeyJustPressed(ebiten.KeyB) {
+		for _, sim := range g.fluidsim10 {
+			sim.ToggleShape()
 		}
 	}
 
-func (g *Game) ResolveCollisionsB(particle *elements.Particle) {
-	//construct search quadrant from current particle
-	quadrant := quadtree.Quadrant{
-		Position:   particle.Position,
-		Dimensions: particle.GetDimensions(),
-	}
-
-	//find list of possible maybe collisions, we inspect those in more detail
-	maybes := g.quadtree.FindAll(quadrant)
-
-	sqdist := float64(GameParticleRadius*GameParticleRadius) * 4
-
-	for _, p := range maybes {
-		if p == particle {
-			continue
-		}
-
-		pos := p.GetPosition()
-		delta := gamedata.Vector{
-			X: pos.X - particle.Position.X,
-			Y: pos.Y - particle.Position.Y,
-		}
-
-		dist2 := delta.X*delta.X + delta.Y*delta.Y
-
-		if dist2 < sqdist {
-			d := math.Sqrt(dist2)
-			overlap := GameParticleRadius*2 - d
-			theta := math.Atan2(delta.Y, delta.X)
-			pos.X += overlap * math.Cos(theta)
-			pos.Y += overlap * math.Sin(theta)
-			p.SetPosition(pos)
-
-			m := p.(*elements.Particle)
-			m.Velocity.X *= -1 * GameDamping
-			m.Velocity.Y *= -1 * GameDamping
-
+	if inpututil.IsKeyJustPressed(ebiten.KeyV) {
+		for _, sim := range g.fluidsim10 {
+			sim.ToggleParticles()
 		}
 	}
+
+}
+
+func (g *Game) Initialize() {
+
+	//10MP Fluid Simulation Initialization
+	g.fluidsim10 = append(g.fluidsim10, elements.NewFluidSim10())
+	g.fluidsim10 = append(g.fluidsim10, elements.NewFluidSim10())
+
+	g.fluidsim10width = float64(g.fluidsim10[0].GetSprite().Bounds().Dx())
+	g.fluidsim10height = float64(g.fluidsim10[0].GetSprite().Bounds().Dy())
+
+	x0 := float64(GameWidth / (len(g.fluidsim10) + 1))
+
+	for i, sim := range g.fluidsim10 {
+
+		pos := gamedata.Vector{
+			X: x0 * float64(i+1),
+			Y: GameHeight / 2.,
+		}
+		sim.SetPosition(pos)
+	}
+
+	//Flask Initialization
+	pos := gamedata.Vector{
+		X: GameWidth / 2,
+		Y: GameHeight / 2,
+	}
+	g.flask.SetPosition(pos)
+	g.flaskcolor = color.RGBA{R: 0xff, G: 0x00, B: 0x00, A: 0xff}
+	g.flask.SetFluidColor(g.flaskcolor)
+
+}
+
+func (g *Game) ManageFlaskInputs() {
+	if inpututil.IsKeyJustPressed(ebiten.KeyR) {
+		g.flask.Initialize()
+		g.mfangle = 0
+		g.flask.SetAngle(g.mfangle)
+		g.flask.SetFluidColor(g.flaskcolor)
+	}
+
+	if inpututil.IsMouseButtonJustPressed(ebiten.MouseButtonLeft) {
+		//angle at mouse down is the zero angle
+		g.mdown = true
+		g.mdx, g.mdy = ebiten.CursorPosition()
+
+		dx := float64(g.mdx) - g.flask.GetPosition().X
+		dy := float64(g.mdy) - g.flask.GetPosition().Y
+		g.mdangle = math.Atan2(dy, dx)
+	}
+
+	if ebiten.IsMouseButtonPressed(ebiten.MouseButtonLeft) {
+
+		if g.mdown {
+			mx, my := ebiten.CursorPosition()
+
+			dx := float64(mx) - g.flask.GetPosition().X
+			dy := float64(my) - g.flask.GetPosition().Y
+			angle := math.Atan2(dy, dx)
+
+			g.flask.SetAngle(g.mfangle + (angle - g.mdangle))
+		}
+	}
+
+	if inpututil.IsMouseButtonJustReleased(ebiten.MouseButtonLeft) {
+		g.mdown = false
+
+		mx, my := ebiten.CursorPosition()
+
+		dx := float64(mx) - g.flask.GetPosition().X
+		dy := float64(my) - g.flask.GetPosition().Y
+		angle := math.Atan2(dy, dx)
+
+		g.mfangle = g.mfangle + (angle - g.mdangle)
+	}
+
+	if inpututil.IsKeyJustPressed(ebiten.KeyM) {
+		g.flask.ToggleBoundaryMask()
+	}
+}
+
+func (g *Game) UpdateFlask() {
+	g.flask.Update()
+}
+
+func (g *Game) RenderFlask(img *ebiten.Image) {
+
+	g.flask.Draw()
+
+	angle := g.flask.GetAngle()
+	dim := g.flask.GetDimensions()
+
+	// //TODO: use flask position for rendering, not screen
+	pos := g.flask.GetPosition()
+
+	op := &ebiten.DrawImageOptions{}
+	op.GeoM.Translate(-dim.X/2, -dim.Y/2)
+	op.GeoM.Rotate(angle)
+	op.GeoM.Scale(2, 2)
+	op.GeoM.Translate(pos.X, pos.Y)
+
+	img.DrawImage(g.flask.GetSprite(), op)
+
 }

gamedata/vector.go

@@ -1,6 +1,28 @@
 package gamedata
 
+import "math"
+
 type Vector struct {
 	X float64
 	Y float64
 }
+
+func (v Vector) DotProduct(p Vector) float64 {
+	return v.X*p.X + v.Y + p.Y
+}
+
+func (v Vector) Magnitude() float64 {
+	return math.Sqrt(v.X*v.X + v.Y*v.Y)
+}
+
+func (v Vector) Add(p Vector) Vector {
+	return Vector{X: v.X + p.X, Y: v.Y + p.Y}
+}
+
+func (v Vector) Subtract(p Vector) Vector {
+	return Vector{X: v.X - p.X, Y: v.Y - p.Y}
+}
+
+func (v Vector) Scale(s float64) Vector {
+	return Vector{X: v.X * s, Y: v.Y * s}
+}

go.mod

@@ -8,7 +8,11 @@ require (
 	github.com/ebitengine/gomobile v0.0.0-20250923094054-ea854a63cce1 // indirect
 	github.com/ebitengine/hideconsole v1.0.0 // indirect
 	github.com/ebitengine/purego v0.9.0 // indirect
+	github.com/go-text/typesetting v0.3.0 // indirect
 	github.com/jezek/xgb v1.1.1 // indirect
+	github.com/rivo/uniseg v0.4.7 // indirect
+	golang.org/x/image v0.31.0 // indirect
 	golang.org/x/sync v0.17.0 // indirect
 	golang.org/x/sys v0.36.0 // indirect
+	golang.org/x/text v0.29.0 // indirect
 )

licenses/TenMinutePhysics

@@ -0,0 +1,11 @@
+Copyright 2022 Matthias Müller - Ten Minute Physics,
+www.youtube.com/c/TenMinutePhysics
+www.matthiasMueller.info/tenMinutePhysics
+
+MIT License
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

main.go

@@ -2,6 +2,7 @@ package main
 
 import (
 	"fluids/game"
+	"fluids/resources"
 	"fmt"
 	"log"
 
@@ -16,6 +17,10 @@ const (
 func main() {
 	fmt.Println("fluid experiments")
 
+	//preload assets
+	resources.LoadImages()
+
+	//initialize new game instance
 	g := game.NewGame()
 
 	ebiten.SetWindowTitle("fluids")

quadtree/quadtree.go

@@ -6,6 +6,7 @@ import (
 )
 
 const (
+	QuadtreeMaxDepth     = 20
 	QuadtreeMaxColliders = 40
 )
 
@@ -18,11 +19,13 @@ type Quadtree struct {
 	quadrant  Quadrant
 	children  []*Quadtree
 	colliders []colliders.Collider
+	depth     int
 }
 
-func New(quadrant Quadrant) *Quadtree {
+func New(quadrant Quadrant, depth int) *Quadtree {
 	qt := &Quadtree{
 		quadrant: quadrant,
+		depth:    depth,
 	}
 	return qt
 }
@@ -101,12 +104,16 @@ func (q *Quadtree) SubdivideAndInsert(obj colliders.Collider) bool {
 
 	var result bool = false
 
-	//initialize up children
+	if q.depth+1 > QuadtreeMaxDepth {
+		return result
+	}
+
+	//initialize children
 	q.children = q.children[:0]
-	q.children = append(q.children, New(Quadrant{Position: gamedata.Vector{X: q.quadrant.Position.X - q.quadrant.Dimensions.X/4, Y: q.quadrant.Position.Y - q.quadrant.Dimensions.Y/4}, Dimensions: gamedata.Vector{X: q.quadrant.Dimensions.X / 2, Y: q.quadrant.Dimensions.Y / 2}}))
-	q.children = append(q.children, New(Quadrant{Position: gamedata.Vector{X: q.quadrant.Position.X + q.quadrant.Dimensions.X/4, Y: q.quadrant.Position.Y - q.quadrant.Dimensions.Y/4}, Dimensions: gamedata.Vector{X: q.quadrant.Dimensions.X / 2, Y: q.quadrant.Dimensions.Y / 2}}))
-	q.children = append(q.children, New(Quadrant{Position: gamedata.Vector{X: q.quadrant.Position.X - q.quadrant.Dimensions.X/4, Y: q.quadrant.Position.Y + q.quadrant.Dimensions.Y/4}, Dimensions: gamedata.Vector{X: q.quadrant.Dimensions.X / 2, Y: q.quadrant.Dimensions.Y / 2}}))
-	q.children = append(q.children, New(Quadrant{Position: gamedata.Vector{X: q.quadrant.Position.X + q.quadrant.Dimensions.X/4, Y: q.quadrant.Position.Y + q.quadrant.Dimensions.Y/4}, Dimensions: gamedata.Vector{X: q.quadrant.Dimensions.X / 2, Y: q.quadrant.Dimensions.Y / 2}}))
+	q.children = append(q.children, New(Quadrant{Position: gamedata.Vector{X: q.quadrant.Position.X - q.quadrant.Dimensions.X/4, Y: q.quadrant.Position.Y - q.quadrant.Dimensions.Y/4}, Dimensions: gamedata.Vector{X: q.quadrant.Dimensions.X / 2, Y: q.quadrant.Dimensions.Y / 2}}, q.depth+1))
+	q.children = append(q.children, New(Quadrant{Position: gamedata.Vector{X: q.quadrant.Position.X + q.quadrant.Dimensions.X/4, Y: q.quadrant.Position.Y - q.quadrant.Dimensions.Y/4}, Dimensions: gamedata.Vector{X: q.quadrant.Dimensions.X / 2, Y: q.quadrant.Dimensions.Y / 2}}, q.depth+1))
+	q.children = append(q.children, New(Quadrant{Position: gamedata.Vector{X: q.quadrant.Position.X - q.quadrant.Dimensions.X/4, Y: q.quadrant.Position.Y + q.quadrant.Dimensions.Y/4}, Dimensions: gamedata.Vector{X: q.quadrant.Dimensions.X / 2, Y: q.quadrant.Dimensions.Y / 2}}, q.depth+1))
+	q.children = append(q.children, New(Quadrant{Position: gamedata.Vector{X: q.quadrant.Position.X + q.quadrant.Dimensions.X/4, Y: q.quadrant.Position.Y + q.quadrant.Dimensions.Y/4}, Dimensions: gamedata.Vector{X: q.quadrant.Dimensions.X / 2, Y: q.quadrant.Dimensions.Y / 2}}, q.depth+1))
 
 	//move colliders into child nodes
 	var failed bool = false

resources/images.go

@@ -0,0 +1,49 @@
+package resources
+
+import (
+	"bytes"
+	"image"
+	"log"
+
+	_ "embed"
+
+	"github.com/hajimehoshi/ebiten/v2"
+)
+
+type ImageName string
+
+const (
+	RoundedBottomFlaskBase          ImageName = "RoundedBottomFlaskBase"
+	RoundedBottomFlaskHighlights    ImageName = "RoundedBottomFlaskHighlights"
+	RoundedBottomFlaskBoundaryMap   ImageName = "RoundedBottomFlaskBoundaryMap"
+	RoundedBottomFlaskBoundaryMap46 ImageName = "RoundedBottomFlaskBoundaryMap46"
+)
+
+var (
+	ImageBank map[ImageName]*ebiten.Image
+
+	//go:embed rounded_bottom_flask_base.png
+	rounded_bottom_flask_base []byte
+	//go:embed rounded_bottom_flask_highlights.png
+	rounded_bottom_flask_highlight []byte
+	//go:embed rounded_bottom_flask_boundary_map.png
+	rounded_bottom_flask_boundary_map []byte
+	//go:embed rounded_bottom_flask_boundary_map_46.png
+	rounded_bottom_flask_boundary_map_46 []byte
+)
+
+func LoadImages() {
+	ImageBank = make(map[ImageName]*ebiten.Image)
+	ImageBank[RoundedBottomFlaskBase] = LoadImagesFatal(rounded_bottom_flask_base)
+	ImageBank[RoundedBottomFlaskHighlights] = LoadImagesFatal(rounded_bottom_flask_highlight)
+	ImageBank[RoundedBottomFlaskBoundaryMap] = LoadImagesFatal(rounded_bottom_flask_boundary_map)
+	ImageBank[RoundedBottomFlaskBoundaryMap46] = LoadImagesFatal(rounded_bottom_flask_boundary_map_46)
+}
+
+func LoadImagesFatal(b []byte) *ebiten.Image {
+	img, _, err := image.Decode(bytes.NewReader(b))
+	if err != nil {
+		log.Fatal(err)
+	}
+	return ebiten.NewImageFromImage(img)
+}

utils/clamp.go

@@ -0,0 +1,27 @@
+package utils
+
+// returns adjusted value between start and end, whichever is closer
+// unless value is between, in which case it returns value
+func Clamp(value, start, end int) int {
+	if value < start {
+		return start
+	}
+
+	if value < end {
+		return value
+	}
+
+	return end
+}
+
+func Clamp32(value, start, end float32) float32 {
+	if value < start {
+		return start
+	}
+
+	if value < end {
+		return value
+	}
+
+	return end
+}