with

    Datatypes,
    ANSI_Terminal,
    Ada.Characters.Latin_1,
    Ada.Text_IO;

use

    Datatypes;

use type

    Datatypes.Plane.Complex;

procedure Fluid_Simulator is

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



    Particles : Particle_Vector := Particle_Vectors.Empty_Vector;

    --  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 Plane.Complex := Plane.Compose_From_Cartesian (0.0, 1.0);
    Pressure_Factor  : constant Quantity := 4.0;
    Viscosity_Factor : constant Quantity := 8.0;



    procedure Read_Input
           (Store : out Particle_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;



    procedure Calculate_Density
           (Store : in out Particle_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 := Plane.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_Vector)
    is
        Displacement, Pressure, Viscosity : Plane.Complex;
        Rij : Quantity;
    begin
        for P of Store loop
            P.Acceleration := Gravity_Factor;
            for Q of Store loop
                Displacement := P.Place - Q.Place;
                Rij := Plane.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 +
                        Plane.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_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_Vector;
            Threshold : in     Quantity) is
    begin
        for C in reverse Store.First_Index .. Store.Last_Index loop
            if Plane.Re (Store (C).Place) < 1.0 - Threshold or
               Plane.Re (Store (C).Place) > 80.0 + Threshold or
               Plane.Im (Store (C).Place) < 1.0 - Threshold or
               Plane.Im (Store (C).Place) > 50.0 + Threshold
            then
                Store.Delete (C);
            end if;
        end loop;
    end Cull_Outside_Bounds;


begin

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

end Fluid_Simulator;