Attempt at tectonic plate based terrain generation

Program.cs

using System;
using System.Collections.Generic;
using System.Numerics;
using System.Drawing;
using System.Drawing.Imaging;
using System.Linq;
using System.Text;
using System.Threading;
using System.Threading.Tasks;

namespace FlowTest
{
    class Program
    {
        // Currently Z axis is in kilometers, X and Y are in pixels (arbitrary)

        const int WIDTH = 200;
        const int HEIGHT = 200;
        const int NUMBER_OF_PLATES = 32;
        const double MANTLE_DENSITY = 3.3; // kg/m^3
        const double OCEANIC_CRUST_START_THICKNESS = 6; // km
        const double OCEANIC_CRUST_START_DENSITY = 3; // kg/m^3
        const double OCEANIC_CRUST_START_TOP = 1; // km
        const double TOTAL_WATER_VOLUME = WIDTH * HEIGHT / 2;

        const int EROSION_WALKS_PER_ITERATION = WIDTH * HEIGHT;
        const double EROSION_AMOUNT_ROCK = 0.0049846;
        const double EROSION_AMOUNT_SEDIMENT = 0.0103156;
        const double EROSION_DEPOSIT_SEDIMENT = 0.00332105;
        const double EROSION_DEPOSIT_SEDIMENT_UNDERWATER = 0.014216;


        static Random rand = new Random();

        static List<Vector2> hotspots = new List<Vector2>();
        static Vector2[] plateForces = new Vector2[NUMBER_OF_PLATES];
        static Vector2[] oldPlateForces = new Vector2[NUMBER_OF_PLATES];
        static double[] plateMass = new double[NUMBER_OF_PLATES];
        static int frame = 0;
        static double waterLevel = 3;

        private static int CountNorth = 0;
        private static int CountEast = 0;
        private static int CountSouth = 0;
        private static int CountWest = 0;

        static void Main(string[] args)
        {
            double[,] mantleHeat = new double[HEIGHT, WIDTH];
            double[,] mantleHeatVelocity = new double[HEIGHT, WIDTH];

            Vector2[,] tectonicForces = new Vector2[HEIGHT, WIDTH];

            Crust[,] crust = new Crust[HEIGHT, WIDTH];
            double[,] crustShifts = new double[HEIGHT, WIDTH];

            Bitmap b = new Bitmap(WIDTH, HEIGHT, PixelFormat.Format24bppRgb);

            mantleHeat.Do((m, y, x) => m[y, x] = Math.Sin(x / 15.234) + Math.Cos(y / 13.645));
            mantleHeatVelocity.Do((m, y, x) => m[y, x] = 0);

            // Mantle bootstrap
            Console.WriteLine("Bootstrapping the mantle...");
            for (int i = 0; i < 200; i++)
            {
                UpdateMantleHeat(ref mantleHeat, ref mantleHeatVelocity);
                DrawMantleHeat(b, mantleHeat);

                frame++;
            }

            // Form a crust
            Console.WriteLine("Forming a crust...");
            crust.Do((c, y, x) =>
            {
                var averageHeat = GetAverageOfCross(mantleHeat, x, y);
                c[y, x].Thickness = OCEANIC_CRUST_START_THICKNESS + 3 * averageHeat;
                c[y, x].Density = OCEANIC_CRUST_START_DENSITY - averageHeat / 10;
            });

            SplitCrustIntoPlates(crust);

            for (int i = 0; i < 150; i++)
            {
                UpdateMantleHeat(ref mantleHeat, ref mantleHeatVelocity);
                DrawMantleHeat(b, mantleHeat);

                if (frame < 320)
                    UpdateTectonicForces(tectonicForces, mantleHeat);
                else
                {
                    tectonicForces.Do((t, y, x) => { t[y, x] = Vector2.Zero; });
                }

                UpdateCrust(crust, tectonicForces, mantleHeat, crustShifts);
                UpdateWeather(crust);
                DrawCrust(b, crust);
                DrawCrustForce(b, crust);

                WriteStats(crust, tectonicForces, mantleHeat);

                frame++;
            }

            DrawHeightmap(b, crust);
        }

        private static void WriteStats(Crust[,] crust, Vector2[,] tectonicForces, double[,] mantleHeat)
        {
            var count = 0;

            var minThickness = double.PositiveInfinity;
            var sumThickness = 0.0;
            var maxThickness = double.NegativeInfinity;

            var minDensity = double.PositiveInfinity;
            var sumDensity = 0.0;
            var maxDensity = double.NegativeInfinity;

            var minTop = double.PositiveInfinity;
            var sumTop = 0.0;
            var maxTop = double.NegativeInfinity;

            crust.Do((c, y, x) =>
            {
                var cr = c[y, x];

                count++;

                sumThickness += cr.Thickness;
                sumDensity += cr.Density;
                sumTop += cr.Top;

                if (cr.Thickness > maxThickness) maxThickness = cr.Thickness;
                if (cr.Thickness < minThickness) minThickness = cr.Thickness;

                if (cr.Density > maxDensity) maxDensity = cr.Density;
                if (cr.Density < minDensity) minDensity = cr.Density;

                if (cr.Top > maxTop) maxTop = cr.Top;
                if (cr.Top < minTop) minTop = cr.Top;
            });

            Console.WriteLine($"Thickness: ({minThickness} | {sumThickness / count} | {maxThickness})");
            Console.WriteLine($"Density: ({minDensity} | {sumDensity / count} | {maxDensity})");
            Console.WriteLine($"Top: ({minTop} | {sumTop / count} | {maxTop})");
            Console.WriteLine($"Water level: {waterLevel}");
            Console.WriteLine($"{CountNorth} {CountEast} {CountSouth} {CountWest}");
        }

        private static void SplitCrustIntoPlates(Crust[,] crust)
        {
            Console.WriteLine("Splitting into plates...");
            List<Point> plateSeeds = new List<Point>();
            for (int i = 0; i < NUMBER_OF_PLATES; i++)
            {
                plateSeeds.Add(new Point(rand.Next(HEIGHT), rand.Next(WIDTH)));
            }

            crust.Do((c, y, x) =>
            {
                var minDist = double.PositiveInfinity;
                var id = 0;
                for (int i = 0; i < NUMBER_OF_PLATES; i++)
                {
                    var dist = Math.Sqrt(Math.Pow(x - plateSeeds[i].X, 2) + Math.Pow(y - plateSeeds[i].Y, 2));
                    if (dist < minDist)
                    {
                        minDist = dist;
                        id = i;
                    }
                }

                crust[y, x].PlateId = id;
            });
        }

        private static void UpdateWeather(Crust[,] crust)
        {
            // Recalculate water. Try to balance it out to a set volume by adjusting water level.
            var totalVolume = 0.0;
            crust.Do((c, y, x) =>
            {
                totalVolume += Math.Max(0, waterLevel - c[y, x].Top);
            });
            if (totalVolume > TOTAL_WATER_VOLUME)
            {
                waterLevel -= 0.0625;
            }
            if (totalVolume < TOTAL_WATER_VOLUME)
            {
                waterLevel += 0.0625;
            }

            // TODO: Wind calculation

            // Erosion and sediment redistribution
            for (int i = 0; i < EROSION_WALKS_PER_ITERATION; i++)
            {
                var x = rand.Next(WIDTH);
                var y = rand.Next(HEIGHT);

                ErosionWalk(crust, y, x);
            }
        }

        private static void ErosionWalk(Crust[,] crust, int y, int x)
        {
            var currentSediment = 0.0;

            var maxWalk = 50;
            while (maxWalk > 0)
            {
                var aboveWater = crust[y, x].TotalTop() > waterLevel - 0.125;

                // Only pick up sediment above water
                if (aboveWater)
                {
                    // Pick up some sediment
                    var sedimentAmount = Math.Min(crust[y, x].Sediment, rand.NextDouble() * EROSION_AMOUNT_SEDIMENT);
                    crust[y, x].Sediment -= sedimentAmount;
                    currentSediment += sedimentAmount;

                    // Deposit some sediment
                    var depositAmount = Math.Min(currentSediment, rand.NextDouble() * EROSION_DEPOSIT_SEDIMENT);
                    currentSediment -= depositAmount;

                    // Erode some rock if it is on the surface
                    if (crust[y, x].Sediment < EROSION_AMOUNT_SEDIMENT && crust[y, x].Thickness > 1)
                    {
                        var erodeAmount = rand.NextDouble() * EROSION_AMOUNT_ROCK;
                        crust[y, x].Thickness -= erodeAmount;
                        currentSediment += erodeAmount;
                    }

                    // Deposit after rock erosion. Give rock erosion a chance.
                    crust[y, x].Sediment += depositAmount;
                }
                else
                {
                    // Just deposit a bunch of sediment
                    var depositAmount = Math.Min(currentSediment, EROSION_DEPOSIT_SEDIMENT_UNDERWATER);
                    crust[y, x].Sediment += depositAmount;
                    currentSediment -= depositAmount;

                    // If we run out then we can stop
                    if (currentSediment <= 0) return;
                }


                var north = y > 0 && crust[y - 1, x].TotalTop() < crust[y, x].TotalTop();
                var east = x < WIDTH - 1 && crust[y, x + 1].TotalTop() < crust[y, x].TotalTop();
                var south = y < HEIGHT - 1 && crust[y + 1, x].TotalTop() < crust[y, x].TotalTop();
                var west = x > 0 && crust[y, x - 1].TotalTop() < crust[y, x].TotalTop();

                if (!(north || east || south || west))
                {
                    crust[y, x].Sediment += currentSediment;
                    return;
                }

                var found = false;
                while (!found)
                {
                    switch (rand.Next(4))
                    {
                        case 0:
                            if (north) { y -= 1; found = true; CountNorth++; }
                            break;
                        case 1:
                            if (east) { x += 1; found = true; CountEast++; }
                            break;
                        case 2:
                            if (south) { y += 1; found = true; CountSouth++; }
                            break;
                        case 3:
                            if (west) { x -= 1; found = true; CountWest++; }
                            break;
                    }
                }

                maxWalk--;
            }
        }

        private static void UpdateCrust(Crust[,] crust, Vector2[,] tectonicForces, double[,] mantleHeat, double[,] crustShifts)
        {
            // Undensify (totally a word) by fractional differentiation of heavier elements
            // Also standardise the thickness
            crust.Do((c, y, x) =>
            {
                var densityChange = (2.4 - c[y, x].Density) * (mantleHeat[y, x] / 50);
                ChangeDensity(c, y, x, densityChange);
                //c[y, x].Thickness -= (c[y, x].Thickness - 15 + rand.NextDouble() * 15) / 5000;

                if (c[y, x].Top < 0.2) c[y, x].Top = 0.2;
            });

            // Move crust based on isostasy
            crust.Do((c, y, x) =>
            {
                var idealHeight = (1 - crust[y, x].Density / 3.3) * crust[y, x].Thickness;
                crust[y, x].Force.Z += (float)(idealHeight - crust[y, x].Top) / 4;
            });

            // Apply tectonic forces from convection (super hack)
            crust.Do((c, y, x) =>
            {
                var v = tectonicForces[y, x];
                c[y, x].Force.X += v.X;
                c[y, x].Force.Y += v.Y;
            });

            // Distribute forces
            for (int i = 0; i < 5; i++)
            {
                crust.Do((c, y, x) =>
                {
                    DistributeForce(c, y, x);
                });
            }

            var temp = oldPlateForces;
            oldPlateForces = plateForces;
            plateForces = temp;
            for (int i = 0; i < NUMBER_OF_PLATES; i++)
            {
                plateMass[i] = 0;
                plateForces[i] = default(Vector2);
            }

            // Resolve forces
            crust.Do((c, y, x) =>
            {
                var upForce = c[y, x].Force.Z;
                c[y, x].Top += upForce;
                c[y, x].Force.Z -= upForce;

                Vector2 lateralForce = new Vector2(c[y, x].Force.X, c[y, x].Force.Y);
                plateForces[c[y, x].PlateId] += lateralForce;
                plateMass[c[y, x].PlateId] += 1;
            });

            for (int i = 0; i < NUMBER_OF_PLATES; i++)
                plateForces[i] = plateForces[i] / (float)plateMass[i] + oldPlateForces[i];

            bool doing = true;
            while (doing)
            {
                doing = false;
                for (int i = 0; i < NUMBER_OF_PLATES; i++)
                {
                    if (plateForces[i].X >= 1)
                    {
                        ShiftPlate(crust, i, 1, 0);
                        plateForces[i].X -= 1;
                        doing = true;
                    }
                    if (plateForces[i].X <= -1)
                    {
                        ShiftPlate(crust, i, -1, 0);
                        plateForces[i].X += 1;
                        doing = true;
                    }
                    if (plateForces[i].Y >= 1)
                    {
                        ShiftPlate(crust, i, 0, 1);
                        plateForces[i].Y -= 1;
                        doing = true;
                    }
                    if (plateForces[i].Y <= -1)
                    {
                        ShiftPlate(crust, i, 0, -1);
                        plateForces[i].Y += 1;
                        doing = true;
                    }
                }
            }


            // TODO: Calculate plate splits

            // Smooth out the tops
            crust.Do((c, y, x) =>
            {
                crustShifts[y, x] = crust[y, x].Top;
            });

            crust.Do((c, y, x) =>
            {
                crust[y, x].Top += (GetAverageOfCross(crustShifts, x, y) - crust[y, x].Top) / 100;
            });
        }

        private static void ChangeDensity(Crust[,] c, int y, int x, double densityChange)
        {
            c[y, x].Thickness = (c[y, x].Thickness * c[y, x].Density) / (c[y, x].Density + densityChange);
            c[y, x].Density += densityChange;
        }

        private static void ShiftPlate(Crust[,] crust, int plateId, int xOff, int yOff)
        {
            // Stops continents going "off screen"
            var fuckoff = false;
            crust.Do((c, y, x) =>
            {
                if (c[y, x].PlateId == plateId && (
                    y + yOff < 0 || y + yOff >= WIDTH ||
                    x + xOff < 0 || x + xOff >= HEIGHT)
                )
                {
                    fuckoff = true;
                }
            });
            if (fuckoff) return;

            crust.Do((c, y, x) =>
            {
                if (c[y, x].PlateId == plateId)
                {
                    var donorX = x - xOff;
                    var donorY = y - yOff;

                    if (donorX >= 0 && donorX < WIDTH && donorY >= 0 && donorY < HEIGHT && crust[donorY, donorX].PlateId == crust[y, x].PlateId)
                    {
                        var toX = x + xOff;
                        var toY = y + yOff;

                        // Subduct the lower crust more often
                        var p = Math.Pow((Math.Tanh(c[toY, toX].Top - c[y, x].Top) + 1) / 2, 5);

                        if (rand.NextDouble() < p)
                        {
                            Subduct(c, toY, toX, c[y, x]);
                        }
                        else
                        {
                            var temp = c[toY, toX];
                            c[toY, toX] = c[y, x];

                            Subduct(c, toY, toX, temp);
                        }

                        crust[y, x] = crust[donorY, donorX];
                    }
                    else
                    {
                        // Spawn new oceanic crust
                        crust[y, x] = new Crust()
                        {
                            Density = OCEANIC_CRUST_START_DENSITY,
                            Thickness = OCEANIC_CRUST_START_THICKNESS,
                            Top = OCEANIC_CRUST_START_TOP,
                            PlateId = crust[y, x].PlateId
                        };
                    }
                }
            }, xOff > 0 || yOff > 0); // Go backwards if we are shifting x or y forwards
        }

        private static void Subduct(Crust[,] c, int y, int x, Crust crust)
        {
            // Give the crust being subducted under a bit of a bump
            c[y, x].Thickness += crust.Top / 100;

            // And give it a bit of lower density magma
            var densityChange = (2.4 - c[y, x].Density) / 50;
            ChangeDensity(c, y, x, densityChange);
        }

        private static void DistributeForce(Crust[,] crust, int y, int x)
        {
            var plateId = crust[y, x].PlateId;

            var north = y > 0;// && crust[y - 1, x].PlateId == plateId;
            var east = x < WIDTH - 1;// && crust[y, x + 1].PlateId == plateId;
            var south = y < HEIGHT - 1;// && crust[y + 1, x].PlateId == plateId;
            var west = x > 0;// && crust[y, x - 1].PlateId == plateId;

            if (!(north || east || south || west))
                return;

            // Let's distribute half the force because why not.
            var excess = crust[y, x].Force / 2;
            crust[y, x].Force -= excess;

            var splitAmount = (excess / 4);

            if (north) crust[y - 1, x].Force += splitAmount;
            if (east) crust[y, x + 1].Force += splitAmount;
            if (south) crust[y + 1, x].Force += splitAmount;
            if (west) crust[y, x - 1].Force += splitAmount;
        }


        private static double GetAverageOfCross(double[,] m, int x, int y)
        {
            return ((y > 0 ? m[y - 1, x] : m[y, x]) +
                    (y < HEIGHT - 1 ? m[y + 1, x] : m[y, x]) +
                    (x > 0 ? m[y, x - 1] : m[y, x]) +
                    (x < WIDTH - 1 ? m[y, x + 1] : m[y, x]) +
                    m[y, x]) / 5;
        }

        private static void UpdateTectonicForces(Vector2[,] tectonicForces, double[,] mantleHeat)
        {
            mantleHeat.Map(tectonicForces, (m, y, x) =>
            {
                // Hot flows towards cold.
                var xForce = (m[y, x] - (x > 0 ? m[y, x - 1] : m[y, x])) +
                             ((x < WIDTH - 1 ? m[y, x + 1] : m[y, x]) - m[y, x]) * 100;

                var yForce = (m[y, x] - (y > 0 ? m[y - 1, x] : m[y, x])) +
                             ((y < HEIGHT - 1 ? m[y + 1, x] : m[y, x]) - m[y, x]) * 100;

                return new Vector2((float)xForce, (float)yForce);
            });
        }

        static void UpdateMantleHeat(ref double[,] mantleHeat, ref double[,] velocity)
        {
            var vel = velocity;
            var mhe = mantleHeat;

            // Randomly add and remove hotspots. These act kind of like mantle plumes. Which totally might not be a real thing.
            if (rand.NextDouble() > 0.85 + 0.01 * hotspots.Count)
                hotspots.Add(new Vector2(rand.Next(WIDTH), rand.Next(HEIGHT)));
            if (rand.NextDouble() > 0.85 && hotspots.Count > 3)
                hotspots.RemoveAt(0);

            // What's that? Heat doesn't travel in waves? Screw you, if I say it travels in waves then it travels in waves!
            foreach (Vector2 hotspot in hotspots)
            {
                vel[(int)hotspot.Y, (int)hotspot.X] += 1;
            }

            mantleHeat.Do((m, y, x) =>
            {
                vel[y, x] += 
                    (((y > 0 ? m[y - 1, x] : m[y, x]) +
                    (y < HEIGHT - 1 ? m[y + 1, x] : m[y, x]) +
                    (x > 0 ? m[y, x - 1] : m[y, x]) +
                    (x < WIDTH - 1 ? m[y, x + 1] : m[y, x])) / 4
                    - m[y, x]);

                vel[y, x] *= 0.9999;
            });

            velocity.Do((v, y, x) =>
            {
                mhe[y, x] += v[y, x];
            });

            //Normalise the energy in the system to prevent it growing or shrinking
            var totalEnergy = 0.0;
            mhe.Do((m, y, x) => { totalEnergy += Math.Abs(m[y, x]); });

            var adjust = (WIDTH * HEIGHT / 4) / totalEnergy;
            mhe.Do((m, y, x) => { m[y, x] = m[y, x] * adjust; });
        }

        private static void DrawMantleHeat(Bitmap bitmap, double[,] mantleHeat)
        {
            var data = bitmap.LockBits(Rectify(bitmap.Size), ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);

            unsafe
            {
                byte* start = (byte*)data.Scan0;

                mantleHeat.Do((m, y, x) =>
                {
                    var val = m[y, x];

                    byte r = 0;
                    byte g = 0;
                    byte b = 0;

                    if (val < 0)
                    {
                        r = 0;
                        b = (byte)(-val * 100);
                    }
                    else
                    {
                        b = 0;
                        r = (byte)(val * 100);
                    }

                    if (val > 3) g = 255;

                    *(start) = b;
                    *(start + 1) = g;
                    *(start + 2) = r;

                    start += 3;
                });
            }

            bitmap.UnlockBits(data);

            bitmap.Save("mantle" + frame + ".png");
        }


        private static void DrawHeightmap(Bitmap bitmap, Crust[,] crust)
        {
            var data = bitmap.LockBits(Rectify(bitmap.Size), ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);

            var maxTop = 3.0;
            crust.Do((c, y, x) =>
            {
                if (c[y,x].TotalTop() > maxTop)
                    maxTop = c[y,x].TotalTop();
            });

            unsafe
            {
                byte* start = (byte*)data.Scan0;

                crust.Do((m, y, x) =>
                {
                    var c = m[y, x];

                    byte r = 0;
                    byte g = 0;
                    byte b = 0;

                    r = b = g = (byte)((c.TotalTop() / maxTop) * 255);
                    
                    *(start) = b;
                    *(start + 1) = g;
                    *(start + 2) = r;

                    start += 3;
                });
            }

            bitmap.UnlockBits(data);

            bitmap.Save("crust" + frame + ".png");
        }

        private static void DrawCrust(Bitmap bitmap, Crust[,] crust)
        {
            var data = bitmap.LockBits(Rectify(bitmap.Size), ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);

            var maxTop = 3.0;
            var maxDepth = 3.0;
            var maxSediment = 1.0;
            crust.Do((m, y, x) =>
            {
                var c = m[y, x];
                if (waterLevel - c.TotalTop() > maxDepth)
                    maxDepth = waterLevel - c.TotalTop();
                if (c.Top > maxTop)
                    maxTop = c.Top;
                if (c.Sediment > maxSediment)
                    maxSediment = c.Sediment;
            });

            unsafe
            {
                byte* start = (byte*)data.Scan0;

                crust.Do((m, y, x) =>
                {
                    var c = m[y, x];

                    byte r = 0;
                    byte g = 0;
                    byte b = 0;

                    //r = b = g = (byte)(50 + (c.TotalTop() / maxTop) * 120);

                    r = c.TotalTop() >= waterLevel ? (byte)(50 + (c.Top / maxTop) * 120) : (byte)0;

                    var water = waterLevel - c.TotalTop();
                    b = water > 0 ? (byte)(128 - (water / maxDepth) * 120) : (byte)0;

                    g = (byte)((c.Sediment / maxSediment) * 250);

                    *(start) = b;
                    *(start + 1) = g;
                    *(start + 2) = r;

                    start += 3;
                });
            }

            bitmap.UnlockBits(data);

            bitmap.Save("crust" + frame + ".png");
        }

        private static void DrawCrustForce(Bitmap bitmap, Crust[,] crust)
        {
            var data = bitmap.LockBits(Rectify(bitmap.Size), ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);

            var maxForce = 0.1f;
            crust.Do((m, y, x) =>
            {
                var d = crust[y, x].Force.Length();
                if (d > maxForce)
                    maxForce = d;
            });
            maxForce = Math.Max(maxForce, 1);


            unsafe
            {
                byte* start = (byte*)data.Scan0;

                crust.Do((m, y, x) =>
                {
                    var c = m[y, x];

                    var force = c.Force.Length();

                    *(start) = (byte)((force / maxForce) * 250);
                    *(start + 1) = (byte)((force / maxForce) * 250);
                    *(start + 2) = (byte)((force / maxForce) * 250);

                    start += 3;
                });
            }

            bitmap.UnlockBits(data);

            bitmap.Save("crustForces" + frame + ".png");
        }

        private static Rectangle Rectify(Size size)
        {
            return new Rectangle(0, 0, size.Width, size.Height);
        }

        private static void Swap(ref double[,] mantleHeat, ref double[,] temp)
        {
            var swap = temp;
            temp = mantleHeat;
            mantleHeat = swap;
        }


        static void Display(double[,] array)
        {
            array.Do((m, y, x) =>
            {
                if (y < 20 && x < 20)
                {
                    Console.SetCursorPosition(x * 4, y);
                    if (m[y, x] > 0) Console.ForegroundColor = ConsoleColor.Red;
                    if (m[y, x] == 0) Console.ForegroundColor = ConsoleColor.White;
                    if (m[y, x] < 0) Console.ForegroundColor = ConsoleColor.Blue;
                    Console.Write((int)m[y, x] + (x == WIDTH - 1 ? "\n" : ""));
                }
            });
        }
    }

    internal class CrustPointComparer : IComparer<Point>
    {
        private Crust[,] crust;

        public CrustPointComparer(Crust[,] crust)
        {
            this.crust = crust;
        }

        public int Compare(Point a, Point b)
        {
            return crust[a.Y, a.X].CompareTo(crust[b.Y, b.X]);
        }
    }

    internal struct Crust : IComparable<Crust>
    {
        // Most of these are in kilometer units

        public double Top; // Height of the top of solid matter of this block (without sediment) above mantle equilibrium point
        public double Density; // Density of the solid matter in this
        public double Thickness; // Height of solid matter in this block
        public double Sediment; // The amount of sediment on top of this block
        public Vector3 Force; // Direction and magnitude of force being applied to this section of crust
        public int PlateId;

        public override string ToString()
        {
            return $"Density (g/cm^3): {Density}  Thickness (km): {Thickness}  PlateId: {PlateId}";
        }

        public int CompareTo(Crust other)
        {
            return (Top).CompareTo(other.Top);
        }

        public double TotalTop()
        {
            return Top + Sediment;
        }
    }

    static class Extensions
    {
        public static void Map<T, TOut>(this T[,] array, TOut[,] outArray, Func<T[,], int, int, TOut> action)
        {
            for (int y = 0; y < array.GetLength(0); y++)
            {
                for (int x = 0; x < array.GetLength(1); x++)
                {
                    outArray[y, x] = action(array, y, x);
                }
            }
        }

        public static void Do<T>(this T[,] array, Action<T[,], int, int> action)
        {
            Do(array, action, false);
        }

        public static void Do<T>(this T[,] array, Action<T[,], int, int> action, bool backwards)
        {
            if (backwards)
            {
                for (int y = array.GetLength(0) - 1; y >= 0; y--)
                {
                    for (int x = array.GetLength(1) - 1; x >= 0; x--)
                    {
                        action(array, y, x);
                    }
                }
            }
            else
            {
                for (int y = 0; y < array.GetLength(0); y++)
                {
                    for (int x = 0; x < array.GetLength(1); x++)
                    {
                        action(array, y, x);
                    }
                }
            }
        }

        public static double XYSum(this Vector3 me)
        {
            return me.X + me.Y;
        }
    }
}