subtree

using System;
using System.Diagnostics;
using System.Linq;
using TMPro;
using Unity.Collections;
using Unity.Jobs;
using Unity.VisualScripting;
using UnityEngine;
using UnityEngine.Profiling;
using Unity.Mathematics;
using UnityEngine.Serialization;
using Debug = UnityEngine.Debug;
using Random = UnityEngine.Random;

public class Simulation : MonoBehaviour
{
    public struct Particle
    {
        public float2 position, velocity;
        public float mass;
        public Color color;
    }
    public int particleCount;
    public float radius;
    public float mass;
    public float timeScale;
    [Range(1, 2)]
    public float centerConcentration;
    public float galaxySize;
    public Gradient velocityGradient;

    [Range(0, 1)] public float approximationThreshold;
    [Range(0, 1)] public float approximationSmoothing;
    public int treeSize;

    [Header("")]
    public bool centerAttraction;
    public float centerMass;

    public float particleForce;

    public Material material;

    NativeArray<Particle> particles;
    NativeArray<Particle> renderParticles;
    QuadTree quadTree;
    new ParticleRenderer renderer;
    Particle center;
    TextMeshProUGUI latency;
    System.Random random;
    NativeArray<Color> sampledGradient;
    QuadTree[,] trees;

    JobHandle handle;

    public static int HalfScreenSize = 7;

    void Start()
    {
        center = new Particle
        {
            position = Vector2.zero,
            velocity = Vector2.zero,
            mass = centerMass
        };

        random = new System.Random(0);
        particles = new NativeArray<Particle>(particleCount + 1, Allocator.Persistent);
        renderParticles = new NativeArray<Particle>(particleCount + 1, Allocator.Persistent);
        particles[0] = center;
        trees = new QuadTree[treeSize, treeSize];
        sampledGradient = SampleGradient(velocityGradient, 256);

        for (int i = 1; i < particleCount; i++)
        {
            particles[i] = ParticlePositionAndVelocity();
        }

        //latency = GameObject.Find("Latency").GetComponent<TextMeshProUGUI>();
        renderer = new(particles, material, radius);
    }

    Particle ParticlePositionAndVelocity()
    {
        Particle particle;
        while (true)
        {
            // Generate radius using a quadratic distribution for higher central density
            float radius = Mathf.Pow(Random.value, centerConcentration) * galaxySize; // Adjust power for density control
            float angle = Random.Range(0f, Mathf.PI * 2f); // Uniform angle

            // Convert polar coordinates to Cartesian
            Vector2 pos = new Vector2(radius * Mathf.Cos(angle), radius * Mathf.Sin(angle));
            float centerDist = pos.magnitude;

            // Calculate circular orbit velocity
            float force = Mathf.Sqrt(centerMass / centerDist) * particleForce;
            Vector2 tangentDirection = new Vector2(-pos.y, pos.x).normalized; // Tangential direction
            Vector2 velocity = tangentDirection * force;

            // Skip if position is too close to the center to avoid instabilities
            float noise = Mathf.PerlinNoise((pos.x + 1000) * 4, (pos.y + 1000) * 4);
            if (centerDist > 0.2f && (noise > .8f || random.NextDouble() < .1f))
            {
                particle = new Particle
                {
                    position = pos,
                    velocity = velocity,
                    mass = mass,
                    color = velocityGradient.Evaluate(noise + 1 - centerDist / galaxySize)
                };
                break;
            }
        }

        return particle;
    }

    void FixedUpdate()
    {
        RebuildQuadTree();
        float dt = Time.fixedDeltaTime * timeScale;
        UpdateVelocities(dt);
        UpdatePositions(dt);
    }

    void Update()
    {
        renderer.Render(particles);
        foreach (var tree in trees)
        {
            tree.DrawTree();
        }
    }

    void OnDestroy()
    {
        handle.Complete();
        renderer.ReleaseBuffers();
        particles.Dispose();
        renderParticles.Dispose();
        sampledGradient.Dispose();
    }

    void RebuildQuadTree()
    {
        for (int y = 0; y < treeSize; y++)
        {
            for (int x = 0; x < treeSize; x++)
            {
                float size = (float)HalfScreenSize * 2 / treeSize;
                trees[x, y] = new(size, new Vector2(-HalfScreenSize + x * size, -HalfScreenSize + y * size));
            }
        }

        foreach (var particle in particles)
        {
            foreach (var tree in trees)
            {
                tree.Insert(particle);
            }
        }

        foreach (var tree in trees)
        {
            tree.Propagate();
        }
    }

    void UpdateVelocities(float deltaTime)
    {
        for (int i = 0; i < particles.Length; i++)
        {
            Particle p = particles[i];
            p.velocity += CalculateAcceleration(particles[i], deltaTime);
            particles[i] = p;
        }
    }

    float2 CalculateAcceleration(Particle particle, float deltaTime)
    {
        Vector2 acceleration = new();
        foreach (var tree in trees)
        {
            acceleration +=
                tree.CalculateAcceleration(particle, approximationThreshold, approximationSmoothing, deltaTime);
        }

        return acceleration;
    }

    void UpdatePositions(float deltaTime)
    {
        float minPosX = -HalfScreenSize + radius;
        float maxPosX = HalfScreenSize - radius;
        float minPosY = -HalfScreenSize + radius;
        float maxPosY = HalfScreenSize - radius;

        UpdateParticlesJob job = new()
        {
            particles = particles,
            gradient = sampledGradient,
            gradientResolution = sampledGradient.Length,
            deltaTime = deltaTime,
            minPosX = minPosX,
            maxPosX = maxPosX,
            minPosY = minPosY,
            maxPosY = maxPosY,
            rng = new Unity.Mathematics.Random(1)
        };
        job.Schedule(particles.Length, 64).Complete();
    }

    NativeArray<Color> SampleGradient(Gradient gradient, int resolution)
    {
        NativeArray<Color> gradientArray = new(resolution, Allocator.Persistent);

        for (int i = 0; i < resolution; i++)
        {
            gradientArray[i] = gradient.Evaluate((float)i / resolution);
        }

        return gradientArray;
    }
}
using System.Collections.Generic;
using Unity.VisualScripting;
using UnityEngine;
using Unity.Mathematics;
using Unity.Collections;
using Unity.Jobs;
using Unity.Burst;
using Particle = Simulation.Particle;

public struct Node
{
    public Quad boundary;
    public float mass;
    public float2 centerOfMass;
    public int childIndex;
    public int next;
}
public struct Quad
{
    //Bottom left of rectangle
    public Vector2 position;
    public float size;

    public int FindNode(Vector2 pos)
    {
        bool wX = pos.x - position.x < size / 2;
        bool sY = pos.y - position.y < size / 2;

        if (wX && !sY) return 0;
        if (!wX && !sY) return 1;
        if (wX && sY) return 2;
        if(!wX && sY) return 3;

        return 0;
    }

    public bool Intersects(Vector2 pos)
    {
        return pos.x > position.x && pos.x < position.x + size &&
               pos.y > position.y && pos.y < position.y + size;
    }
}

public class QuadTree
{
    public NativeList<Node> nodes;

    public QuadTree(float size, Vector2 position)
    {
        nodes = new(Allocator.Persistent);
        Quad boundary = new() { position = position, size = size };
        nodes.Add(new Node
        {
            boundary = boundary,
            next = 0
        });
    }

    ~QuadTree()
    {
        nodes.Dispose();
    }

    public void DrawTree()
    {
        foreach (var node in nodes)
        {
            DrawSquare(node.boundary.position, node.boundary.size);
        }
    }
    static void DrawSquare(Vector2 position, float size)
    {
        Debug.DrawLine(position, new Vector2(position.x, position.y + size), Color.red);
        Debug.DrawLine(new Vector2(position.x, position.y + size), new Vector2(position.x + size, position.y + size), Color.red);
        Debug.DrawLine(new Vector2(position.x + size, position.y + size), new Vector2(position.x + size, position.y), Color.red);
        Debug.DrawLine(new Vector2(position.x + size, position.y), new Vector2(position.x, position.y), Color.red);
    }

    public Vector2 CalculateAcceleration(Particle particle, float theta, float epsilon, float deltaTime)
    {
        Vector2 acceleration = new();
        float thetaSqr = theta * theta;
        float epsilonSqr = epsilon * epsilon;

        int nodeIndex = 0;
        while (true)
        {
            Node node = nodes[nodeIndex];
            Vector2 diff = node.centerOfMass - particle.position;
            float sqrDst = Mathf.Max(diff.sqrMagnitude, 0.05f);

            if (node.childIndex == 0 || node.boundary.size * node.boundary.size < sqrDst * thetaSqr)
            {
                float denom = (sqrDst + epsilonSqr) * Mathf.Sqrt(sqrDst);
                acceleration += diff * node.mass / denom * deltaTime;

                if (node.next == 0) return acceleration;
                nodeIndex = node.next;
            }
            else nodeIndex = node.childIndex;
        }
    }

    public void Propagate()
    {
        int count = nodes.Length - 1;
        for (int i = count; i > 0; i--)
        {
            int child = nodes[i].childIndex;
            if(child == 0) continue;

            Node node = nodes[i];
            Node nw = nodes[child];
            Node ne = nodes[child + 1];
            Node sw = nodes[child + 2];
            Node se = nodes[child + 3];
            node.centerOfMass = nw.centerOfMass * nw.mass
                                + ne.centerOfMass * ne.mass
                                + sw.centerOfMass * sw.mass
                                + se.centerOfMass * se.mass;
            node.mass = nw.mass
                        + ne.mass
                        + sw.mass
                        + se.mass;
            node.centerOfMass /= node.mass;
            nodes[i] = node;

        }
    }

    public void Insert(Particle particle)
    {
        if (!nodes[0].boundary.Intersects(particle.position)) return;

        int nodeIndex = 0;
        Node node = nodes[0];
        while (node.childIndex != 0)
        {
            int quad = node.boundary.FindNode(particle.position);
            nodeIndex = quad + node.childIndex;
            node = nodes[nodeIndex];
        }

        if (node.mass == 0)
        {
            node.centerOfMass = particle.position;
            node.mass = particle.mass;
            nodes[nodeIndex] = node;
            return;
        }

        if (Equals(node.centerOfMass, particle.position))
        {
            node.mass += particle.mass;
            nodes[nodeIndex] = node;
            return;
        }

        Vector2 nodePos = node.centerOfMass;
        float nodeMass = node.mass;
        while (true)
        {
            Subdivide(nodeIndex);
            node = nodes[nodeIndex];
            int childIndex = node.childIndex;
            int q1 = node.boundary.FindNode(nodePos);
            int q2 = node.boundary.FindNode(particle.position);

            if (q1 == q2)
            {
                nodeIndex = q1 + childIndex;
            }
            else
            {
                int n1 = q1 + childIndex;
                int n2 = q2 + childIndex;

                nodes[n1] = new Node
                {
                    centerOfMass = nodePos,
                    mass = nodeMass,
                    boundary = nodes[n1].boundary,
                    next = nodes[n1].next
                };
                nodes[n2] = new Node
                {
                    centerOfMass = particle.position,
                    mass = particle.mass,
                    boundary = nodes[n2].boundary,
                    next = nodes[n2].next
                };
                return;
            }
        }
    }

    void Subdivide(int parentIndex)
    {
        Node parent = nodes[parentIndex];
        parent.childIndex = nodes.Length;
        nodes[parentIndex] = parent;
        Vector2 parentPos = parent.boundary.position;
        Vector2 nwPosition = new(parentPos.x, parentPos.y + parent.boundary.size / 2);
        Vector2 nePosition = new(parentPos.x + parent.boundary.size / 2, parentPos.y + parent.boundary.size / 2);
        Vector2 swPosition = new(parentPos.x, parentPos.y);
        Vector2 sePosition = new(parentPos.x + parent.boundary.size / 2, parentPos.y);

        Node nw = new()
        {
            boundary = new Quad { position = nwPosition, size = parent.boundary.size / 2 },
            next = parent.childIndex + 1
        };
        Node ne = new()
        {
            boundary = new Quad { position = nePosition, size = parent.boundary.size / 2 },
            next = parent.childIndex + 2
        };
        Node sw = new()
        {
            boundary = new Quad { position = swPosition, size = parent.boundary.size / 2 },
            next = parent.childIndex + 3
        };
        Node se = new()
        {
            boundary = new Quad { position = sePosition, size = parent.boundary.size / 2 },
            next = parent.next
        };

        nodes.Add(nw);
        nodes.Add(ne);
        nodes.Add(sw);
        nodes.Add(se);
    }
}