Initial attempt

1 files changed, 197 insertions, 0 deletions

diff --git a/src/fluid_simulator.adb b/src/fluid_simulator.adb
new file mode 100644
index 0000000..01d3772
--- /dev/null
+++ b/src/fluid_simulator.adb
@@ -0,0 +1,197 @@
+
+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;
+