with

    Ada.Numerics.Generic_Complex_Types,
    Ada.Containers.Vectors,
    Ada.Characters.Latin_1,
    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;



    function Marching_Squares
           (Input : in Particle_Vectors.Vector)
        return String
    is
        Liquid_Chars : String (1 .. 16) := " .,_`/[/']\\-/\#";
        Grid : array (0 .. 80, 0 .. 25) of Boolean := (others => (others => False));
        Output : String (1 .. 2025);
        X, Y : Integer;
    begin
        for P of Input loop
            X := Integer (Fixed.Re (P.Place));
            Y := Integer (Fixed.Im (P.Place));
            if X >= 0 and X <= 80 and Y >= 0 and Y <= 25 then
                Grid (X, Y) := True;
            end if;
        end loop;
        for J in Integer range 1 .. 25 loop
            for I in Integer range 1 .. 80 loop
                Output ((J - 1) * 81 + I) := Liquid_Chars
                   (Boolean'Pos (Grid (I - 1, J)) + 1 +
                    Boolean'Pos (Grid (I, J)) * 2 +
                    Boolean'Pos (Grid (I, J - 1)) * 4 +
                    Boolean'Pos (Grid (I - 1, J - 1)) * 8);
            end loop;
            Output (J * 81) := Latin.LF;
        end loop;
        return Output;
    end Marching_Squares;

begin

    declare
        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 + 1.0;
            else
                if Input > Latin.Space and Input < Latin.DEL then
                    Particles.Append (Create (X, Y, (Input = '#')));
                end if;
                X := X + 1.0;
            end if;
        end loop;
    end;

    loop
        --  Calculate density
        for P of Particles loop
            if P.Solid then
                P.Density := 9.0;
            else
                P.Density := 0.0;
            end if;
            for Q of Particles loop
                declare
                    Rij : Quantity := Fixed.Modulus (P.Place - Q.Place);
                    W : Quantity := (Rij / Particle_Radius - 1.0) ** 2;
                begin
                    if Rij < Particle_Radius then
                        P.Density := P.Density + Particle_Mass * W;
                    end if;
                end;
            end loop;
        end loop;

        --  Calculate interaction forces
        for P of Particles loop
            P.Acceleration := Gravity_Factor;
            for Q of Particles loop
                declare
                    Displacement : Fixed.Complex := P.Place - Q.Place;
                    Rij : Quantity := Fixed.Modulus (Displacement);
                begin
                    if Rij < Particle_Radius then
                        declare
                            Pressure : Fixed.Complex :=
                                (P.Density + Q.Density - 2.0 * P0) * Pressure_Factor * Displacement;
                            Viscosity : Fixed.Complex :=
                                (P.Velocity - Q.Velocity) * Viscosity_Factor;
                        begin
                            P.Acceleration := P.Acceleration +
                                Fixed.Compose_From_Cartesian (1.0 - Rij / Particle_Radius) /
                                P.Density * (Pressure - Viscosity);
                        end;
                    end if;
                end;
            end loop;
        end loop;

        --  Render using marching squares
        Clear_Screen;
        Reset_Cursor;
        IO.Put (Marching_Squares (Particles));

        --  Update acceleration, velocity, and position
        for P of Particles loop
            if not P.Solid then
                P.Velocity := P.Velocity + P.Acceleration / 10.0;
                P.Place := P.Place + P.Velocity;
            end if;
        end loop;

        --  Cull particles going out of bounds before any overflows happen
        for C in reverse Particles.First_Index .. Particles.Last_Index loop
            if Fixed.Re (Particles (C).Place) < -50.0 or
               Fixed.Re (Particles (C).Place) > 130.0 or
               Fixed.Im (Particles (C).Place) < -50.0 or
               Fixed.Im (Particles (C).Place) > 75.0
            then
                Particles.Delete (C);
            end if;
        end loop;

        --  Sleep for a bit
        delay 0.012321;
    end loop;

end Fluid_Simulator;