XFEEDToolBox

using System.Numerics;
using XFEEDToolBox.Core.Models;

namespace XFEEDToolBox.Core.Rendering;

/// <summary>A source MilkyWayLayerJob 28-byte diffuse-cloud record, never a system/star.</summary>
public readonly record struct GalaxyOriginalCloud(Vector3 AbsolutePosition,float Alpha,float Size,float Rotation,int AtlasIndex);
public readonly record struct GalaxyOriginalCloudWindow(int MinX,int MinY,int MinZ,int MaxX,int MaxY,int MaxZ,int Step,int Ordinal,bool Dust);
public sealed record GalaxyOriginalCloudLayer(GalaxyOriginalCloudWindow Window,IReadOnlyList<GalaxyOriginalCloud> Clouds,int SampledCells);

/// <summary>
/// Pure local port of original 3A64C60/3A469D0 and constructor 3A08F00 defaults.
/// The source maps retain exact integers. Original runtime debug overrides,
/// shared output multipliers and the render/material stage are independent.
/// </summary>
public static class GalaxyOriginalClouds
{
    public const int Step=1066,MilkyWayHeight=16000,DustHeight=1280;
    public const float MinimumDensity=.001f,MinimumSize=.75f,MaximumSize=1.5f;

    public static GalaxyOriginalCloudWindow Window(Coordinates focus,int ordinal=0,bool dust=false)
    {
        if(ordinal<0||ordinal>(dust?2:0))throw new ArgumentOutOfRangeException(nameof(ordinal));
        var source=GalaxyOriginalDensity.Reference;var absoluteX=Math.Clamp(focus.X+source.Origin.X,0,100000);
        var absoluteZ=Math.Clamp(focus.Z+source.Origin.Z,0,100000);var step=Step>>ordinal;
        var half=ordinal==0?int.MaxValue:(100000>>(ordinal+1))>>1;
        var lowX=ordinal==0?0:Math.Max(0,(int)absoluteX-half)/step*step;
        var lowZ=ordinal==0?0:Math.Max(0,(int)absoluteZ-half)/step*step;
        var height=dust?DustHeight:MilkyWayHeight;
        return new(lowX,40955-(height>>1),lowZ,ordinal==0?100000:Math.Min(100000,(int)absoluteX+half),
            40955+(height>>1),ordinal==0?100000:Math.Min(100000,(int)absoluteZ+half),step,ordinal,dust);
    }

    public static GalaxyOriginalCloudLayer Generate(GalaxyOriginalCloudWindow window,CancellationToken cancellation=default)
    {
        if(window.Step!=(Step>>window.Ordinal)||window.Ordinal<0||window.Ordinal>(window.Dust?2:0))throw new ArgumentOutOfRangeException(nameof(window));
        var source=GalaxyOriginalDensity.Reference;var result=new List<GalaxyOriginalCloud>(window.Dust?500:10000);var sampled=0;
        var quantization=~((1<<window.Ordinal)-1);
        for(var y=window.MinY;y<window.MaxY;y+=window.Step)
        {
            cancellation.ThrowIfCancellationRequested();
            if(window.Dust&&(y<(81919-window.Step)/2||y>=(81919+window.Step)/2))continue;
            for(var x=window.MinX;x<window.MaxX;x+=window.Step)
            for(var z=window.MinZ;z<window.MaxZ;z+=window.Step)
            {
                var cell=GalaxyOriginalDensityMap.CellAddress(x&quantization,y&quantization,z&quantization);sampled++;
                var seed=Hash(cell);var random=seed;float density;uint atlasRandom,rotationRandom,sizeRandom;
                if(!window.Dust)
                {
                    density=source.StellarDensity(cell);
                    if(density<=MinimumDensity||source.MilkyWaySurface(cell)<=MinimumDensity)continue;
                    atlasRandom=Next(ref random);rotationRandom=Next(ref random);sizeRandom=Next(ref random);
                }
                else
                {
                    var midplane=(cell&0xfffc001fff8UL)|0x20000000;
                    if(source.HeightIntegral(cell)<source.HeightIntegral(midplane)*(1f-95f*.01f))continue;
                    density=source.DustDensity(cell);
                    if(Next(ref random)*(1f/2147483648)>=density*.1f)continue;
                    atlasRandom=Next(ref random);rotationRandom=Next(ref random);sizeRandom=Next(ref random);
                }
                var cx=(uint)(cell>>30)&0x3fff;var cy=(uint)(cell>>17)&0x1fff;var cz=(uint)(cell>>3)&0x3fff;
                var jitter=seed;
                // Original call order is Z, Y, X. Size/angle use a separate
                // sequence from the same cell seed, not the post-Gaussian state.
                var pz=Gaussian(ref jitter,cz*10+5,window.Step);
                var py=Gaussian(ref jitter,cy*10+5,window.Step);
                var px=Gaussian(ref jitter,cx*10+5,window.Step);
                result.Add(new(new(px,py,pz),window.Dust?density*.5f:density,
                    sizeRandom*(1f/2147483648)*(MaximumSize-MinimumSize)+MinimumSize,
                    rotationRandom*(1f/2147483648)*6.2831854820251465f,
                    (window.Dust?3:window.Ordinal*3)+(int)(atlasRandom%3)));
            }
        }
        return new(window,result.AsReadOnly(),sampled);
    }

    public static uint Hash(ulong cell)
    {
        var hash=unchecked(cell*0x40000+~cell);hash=unchecked((hash^(hash>>31))*21);hash=unchecked((hash^(hash>>11))*65);
        var seed=(uint)(hash^(hash>>22));return seed==0?1:seed;
    }
    private static uint Next(ref uint state)
    {
        var quotient=state/44488;var product=(state%44488)*48271;var subtract=quotient*3399;
        state=unchecked(product-subtract+(subtract<product?0U:2147483647U));return state;
    }
    private static float Gaussian(ref uint random,float mean,float sigma)
    {
        float first,squared;
        do
        {
            var u=Next(ref random)*(1f/2147483648);first=u+u-1;
            var v=Next(ref random)*(1f/2147483648);var second=v+v-1;squared=first*first+second*second;
        }while(squared>1||squared==0);
        return MathF.Sqrt(MathF.Log(squared)*-2/squared)*first*sigma+mean;
    }
}