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;



    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;



    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 <= 50 then
                Grid (X, Y / 2) := 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;



    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;