path: root/src/fluid_simulator.adb

with

    Ada.Numerics.Generic_Complex_Types,
    Ada.Containers.Vectors,
    Ada.Characters.Latin_1,
    Ada.Strings.Fixed,
    Ada.Text_IO;

procedure Fluid_Simulator is

    package Latin renames Ada.Characters.Latin_1;
    package IO renames Ada.Text_IO;



    procedure Clear_Screen is
    begin
        --  ANSI control sequence Erase in Display
        --  Variant to clear entire screen
        IO.Put (Latin.ESC & "[2J");
    end Clear_Screen;

    procedure Reset_Cursor is
    begin
        --  ANSI control sequence Cursor Position
        --  Parameters to move cursor to top left corner
        IO.Put (Latin.ESC & "[1;1H");
    end Reset_Cursor;



    type Quantity is digits 18;

    package Fixed is new Ada.Numerics.Generic_Complex_Types (Real => Quantity);
    use type Fixed.Complex;

    type Particle is record
        Place        : Fixed.Complex;
        Solid        : Boolean;
        Density      : Quantity := 0.0;
        Acceleration : Fixed.Complex := Fixed.Compose_From_Cartesian (0.0, 0.0);
        Velocity     : Fixed.Complex := Fixed.Compose_From_Cartesian (0.0, 0.0);
    end record;

    function Create
           (X, Y  : Quantity;
            Solid : in Boolean)
        return Particle is
    begin
        return (Place => Fixed.Compose_From_Cartesian (X, Y), Solid => Solid, others => <>);
    end Create;



    --  Constant properties of particles
    Particle_Radius : constant Quantity := 2.0;
    Particle_Mass   : constant Quantity := 1.0;

    --  Constant used in calculating fluid interaction forces
    P0 : constant Quantity := 1.5;

    --  Other constant force factors
    Gravity_Factor : constant Fixed.Complex := Fixed.Compose_From_Cartesian (0.0, 1.0);
    Pressure_Factor : constant Quantity := 4.0;
    Viscosity_Factor : constant Quantity := 8.0;



    package Particle_Vectors is new Ada.Containers.Vectors
       (Index_Type   => Positive,
        Element_Type => Particle);

    Particles : Particle_Vectors.Vector := Particle_Vectors.Empty_Vector;



    procedure Read_Input
           (Store : out Particle_Vectors.Vector)
    is
        Input : Character;
        X, Y : Quantity := 1.0;
    begin
        while not IO.End_Of_File (IO.Standard_Input) loop
            IO.Get_Immediate (Input);
            if Input = Latin.LF then
                X := 1.0;
                Y := Y + 2.0;
            else
                if Input > Latin.Space and Input < Latin.DEL then
                    Store.Append (Create (X, Y, (Input = '#')));
                    Store.Append (Create (X, Y + 1.0, (Input = '#')));
                end if;
                X := X + 1.0;
            end if;
        end loop;
    end Read_Input;



    --  Liquid_Chars : constant String (1 .. 16) := " .,_`/[/']\\-/\#";
    Liquid_Chars : constant String (1 .. 16) := " ,.-`[//'\]\-\/#";

    type Liquidex is mod 2**4;

    type March_Cell is record
        Index : Liquidex := 0;
        Density : Quantity := 0.0;
    end record;

    type March_Cell_Grid is array (Integer range <>, Integer range <>) of March_Cell;



    function BG_Color_Code
           (Value : in Natural)
        return String
    is
        use Ada.Strings;
        use Ada.Strings.Fixed;
    begin
        --  Total length is always 11 characters
        return Latin.ESC & "[48;5;" & Tail (Trim (Integer'Image (Value), Left), 3, '0') & "m";
    end BG_Color_Code;



    function Lookup
           (Input : in March_Cell_Grid;
            X, Y  : in Integer)
        return String
    is
        Average_Density : Natural := Integer (Quantity'Ceiling (Input (X, Y).Density / 4.0));
        Bit_Index : Positive := Integer (Input (X, Y).Index) + 1;
        Choice : Natural;
    begin
        case Average_Density is
            when 1 .. 2 => Choice := 19; --  dark blue
            when 3 .. 4 => Choice := 20; --  slightly less dark blue
            when 5 .. 6 => Choice := 21; --  slightly dark blue
            when 7 .. 8 => Choice := 12; --  blue
            when 9 .. 10 => Choice := 14; --  cyan
            when 11 .. 12 => Choice := 10; --  green
            when 13 .. 14 => Choice := 11; --  yellow
            when 15 .. 16 => Choice := 3; --  dark yellow
            when 17 .. 18 => Choice := 9; --  red
            when 19 .. 20 => Choice := 1; --  dark red
            when others => Choice := 0; --  black
        end case;
        --  Total length should always be 12 characters
        return BG_Color_Code (Choice) & Liquid_Chars (Bit_Index);
    end Lookup;



    function Marching_Squares
           (Input : in Particle_Vectors.Vector)
        return String
    is
        --  Having the grid be one bigger around the edges simplifies calculations
        Grid : March_Cell_Grid (0 .. 81, 0 .. 26);

        --  80 cols * 25 rows * 12 chars/cell + 24 linefeeds + 4 char color reset = 24028
        --  Oh yeah, baby, big strings
        Output : String (1 .. 24028);

        X, Y, S : Integer;
    begin
        for P of Input loop
            X := Integer (Fixed.Re (P.Place) - 0.5);
            Y := Integer (Fixed.Im (P.Place) / 2.0 - 0.5);
            if X >= 0 and X <= 80 and Y >= 0 and Y <= 25 then
                for J in Integer range 0 .. 1 loop
                    for I in Integer range 0 .. 1 loop
                        Grid (X + I, Y + J).Index :=
                            Grid (X + I, Y + J).Index or (2 ** (I + 2 * J));
                        Grid (X + I, Y + J).Density :=
                            Grid (X + I, Y + J).Density + P.Density;
                    end loop;
                end loop;
            end if;
        end loop;
        for J in Integer range 1 .. 25 loop
            for I in Integer range 1 .. 80 loop
                S := (J - 1) * 961 + (I - 1) * 12 + 1;
                Output (S .. S + 11) := Lookup (Grid, I, J);
            end loop;
            Output (J * 961) := Latin.LF;
        end loop;
        Output (24025 .. 24028) := Latin.ESC & "[0m";
        return Output;
    end Marching_Squares;



    procedure Calculate_Density
           (Store : in out Particle_Vectors.Vector)
    is
        Rij, W : Quantity;
    begin
        for P of Store loop
            P.Density := (if P.Solid then 9.0 else 0.0);
            for Q of Store loop
                Rij := Fixed.Modulus (P.Place - Q.Place);
                W := (Rij / Particle_Radius - 1.0) ** 2;
                if Rij < Particle_Radius then
                    P.Density := P.Density + Particle_Mass * W;
                end if;
            end loop;
        end loop;
    end Calculate_Density;



    procedure Calculate_Interaction
           (Store : in out Particle_Vectors.Vector)
    is
        Displacement, Pressure, Viscosity : Fixed.Complex;
        Rij : Quantity;
    begin
        for P of Store loop
            P.Acceleration := Gravity_Factor;
            for Q of Store loop
                Displacement := P.Place - Q.Place;
                Rij := Fixed.Modulus (Displacement);
                if Rij < Particle_Radius then
                    Pressure := (P.Density + Q.Density - 2.0 * P0) * Pressure_Factor * Displacement;
                    Viscosity := (P.Velocity - Q.Velocity) * Viscosity_Factor;
                    P.Acceleration := P.Acceleration +
                        Fixed.Compose_From_Cartesian (1.0 - Rij / Particle_Radius) /
                        P.Density * (Pressure - Viscosity);
                end if;
            end loop;
        end loop;
    end Calculate_Interaction;



    procedure Update_Position
           (Store : in out Particle_Vectors.Vector) is
    begin
        for P of Store loop
            if not P.Solid then
                P.Velocity := P.Velocity + P.Acceleration / 10.0;
                P.Place := P.Place + P.Velocity;
            end if;
        end loop;
    end Update_Position;



    procedure Cull_Outside_Bounds
           (Store     : in out Particle_Vectors.Vector;
            Threshold : in Quantity) is
    begin
        for C in reverse Store.First_Index .. Store.Last_Index loop
            if Fixed.Re (Store (C).Place) < 1.0 - Threshold or
               Fixed.Re (Store (C).Place) > 80.0 + Threshold or
               Fixed.Im (Store (C).Place) < 1.0 - Threshold or
               Fixed.Im (Store (C).Place) > 50.0 + Threshold
            then
                Store.Delete (C);
            end if;
        end loop;
    end Cull_Outside_Bounds;

begin

    Read_Input (Particles);
    loop
        Clear_Screen;
        Reset_Cursor;
        IO.Put (Marching_Squares (Particles));
        Calculate_Density (Particles);
        Calculate_Interaction (Particles);
        Update_Position (Particles);
        Cull_Outside_Bounds (Particles, 50.0);
        delay 0.012321;
    end loop;

end Fluid_Simulator;