Initial attempt

3 files changed, 843 insertions, 0 deletions

diff --git a/src/complex_fixed_points.adb b/src/complex_fixed_points.adb
new file mode 100644
index 0000000..a0c7be1
--- /dev/null
+++ b/src/complex_fixed_points.adb
@@ -0,0 +1,529 @@
+
+package body Complex_Fixed_Points is
+
+    --  I have no idea what visibility conflict is requiring these two to be defined
+
+    function "*"
+           (Left, Right : in Real'Base)
+        return Real'Base is
+    begin
+        return Standard."*" (Left, Right);
+    end "*";
+
+    function "/"
+           (Left, Right : in Real'Base)
+        return Real'Base is
+    begin
+        return Standard."/" (Left, Right);
+    end "/";
+
+
+
+    function Re
+           (X : in Complex)
+        return Real'Base is
+    begin
+        return X.Re;
+    end Re;
+
+    function Im
+           (X : in Complex)
+        return Real'Base is
+    begin
+        return X.Im;
+    end Im;
+
+    function Im
+           (X : in Imaginary)
+        return Real'Base is
+    begin
+        return Real'Base (X);
+    end Im;
+
+
+
+    procedure Set_Re
+           (X  : in out Complex;
+            Re : in     Real'Base) is
+    begin
+        X.Re := Re;
+    end Set_Re;
+
+    procedure Set_Im
+           (X  : in out Complex;
+            Im : in     Real'Base) is
+    begin
+        X.Im := Im;
+    end Set_Im;
+
+    procedure Set_Im
+           (X  : in out Imaginary;
+            Im : in     Real'Base) is
+    begin
+        X := Imaginary (Im);
+    end Set_Im;
+
+
+
+    function Cartesian
+           (Re, Im : in Real'Base)
+        return Complex is
+    begin
+        return (Re => Re, Im => Im);
+    end Cartesian;
+
+    function Cartesian
+           (Re : in Real'Base)
+        return Complex is
+    begin
+        return (Re => Re, Im => 0.0);
+    end Cartesian;
+
+    function Cartesian
+           (Im : in Imaginary)
+        return Complex is
+    begin
+        return (Re => 0.0, Im => Real'Base (Im));
+    end Cartesian;
+
+
+
+    function Sqrt
+           (X : in Real'Base)
+        return Real'Base
+    is
+        Old_Result : Real'Base;
+        Result : Real'Base := 1.0;
+        Current_Square : Real'Base := 1.0;
+        Next_L : Real'Base := 3.0;
+    begin
+        --  If we don't check for this an overflow will happen due to taking things to the limit
+        if X = 0.0 then
+            return 0.0;
+        end if;
+
+        --  Find the smallest integer greater than Sqrt (X)
+        while Current_Square < X loop
+            Current_Square := Current_Square + Next_L;
+            Next_L := Next_L + 2.0;
+            Result := Result + 1.0;
+        end loop;
+
+        --  Lucky guess
+        if Current_Square = X then
+            return Result;
+        end if;
+
+        --  Babylonian aka Newton's Method
+        --  To be honest I'm kinda suspicious this may not necessarily terminate
+        loop
+            Old_Result := Result;
+            Result := (Old_Result + X / Old_Result) / 2.0;
+            exit when Real (Result) = Real (Old_Result);
+        end loop;
+        return Result;
+    end Sqrt;
+
+    function Modulus
+           (X : in Complex)
+        return Real'Base is
+    begin
+        return Sqrt (X.Re * X.Re + X.Im * X.Im);
+    end Modulus;
+
+
+
+    function "+"
+           (Right : in Complex)
+        return Complex is
+    begin
+        return Right;
+    end "+";
+
+    function "-"
+           (Right : in Complex)
+        return Complex is
+    begin
+        return (Re => -Right.Re, Im => -Right.Im);
+    end "-";
+
+    function Conjugate
+           (X : in Complex)
+        return Complex is
+    begin
+        return (Re => X.Re, Im => -X.Im);
+    end Conjugate;
+
+
+
+    function "+"
+           (Left, Right : in Complex)
+        return Complex is
+    begin
+        return (Re => Left.Re + Right.Re,
+                Im => Left.Im + Right.Im);
+    end "+";
+
+    function "-"
+           (Left, Right : in Complex)
+        return Complex is
+    begin
+        return (Re => Left.Re - Right.Re,
+                Im => Left.Im - Right.Im);
+    end "-";
+
+    function "*"
+           (Left, Right : in Complex)
+        return Complex is
+    begin
+        return (Re => Left.Re * Right.Re - Left.Im * Right.Im,
+                Im => Left.Re * Right.Im + Left.Im * Right.Re);
+    end "*";
+
+    function "/"
+           (Left, Right : in Complex)
+        return Complex
+    is
+        Denominator : Real'Base := Right.Re * Right.Re + Right.Im * Right.Im;
+    begin
+        return (Re => (Left.Re * Right.Re + Left.Im * Right.Im) / Denominator,
+                Im => (Left.Im * Right.Re - Left.Re * Right.Im) / Denominator);
+    end "/";
+
+
+
+    function "**"
+           (Left  : in Complex;
+            Right : in Integer)
+        return Complex
+    is
+        Result : Complex;
+        Exponent : Integer := Right;
+    begin
+        if Exponent >= 1 then
+            Result := Left;
+            while Exponent > 1 loop
+                Result := Result * Left;
+                Exponent := Exponent - 1;
+            end loop;
+        else
+            Result := One;
+            while Exponent < 0 loop
+                Result := Result / Left;
+                Exponent := Exponent + 1;
+            end loop;
+        end if;
+        return Result;
+    end "**";
+
+
+
+    function "+"
+           (Right : in Imaginary)
+        return Imaginary is
+    begin
+        return Right;
+    end "+";
+
+    function "-"
+           (Right : in Imaginary)
+        return Imaginary is
+    begin
+        return Imaginary (-Real'Base (Right));
+    end "-";
+
+    function "abs"
+           (Right : in Imaginary)
+        return Real'Base is
+    begin
+        return abs (Real'Base (Right));
+    end "abs";
+
+
+
+    function "+"
+           (Left, Right : in Imaginary)
+        return Imaginary is
+    begin
+        return Imaginary (Real'Base (Left) + Real'Base (Right));
+    end "+";
+
+    function "-"
+           (Left, Right : in Imaginary)
+        return Imaginary is
+    begin
+        return Imaginary (Real'Base (Left) - Real'Base (Right));
+    end "-";
+
+    function "*"
+           (Left, Right : in Imaginary)
+        return Real'Base is
+    begin
+        return -(Real'Base (Left) * Real'Base (Right));
+    end "*";
+
+    function "/"
+           (Left, Right : in Imaginary)
+        return Real'Base is
+    begin
+        return Real'Base (Left) / Real'Base (Right);
+    end "/";
+
+
+
+    function "**"
+           (Left  : in Imaginary;
+            Right : in Integer)
+        return Complex
+    is
+        Result : Complex;
+        Exponent : Integer := Right;
+    begin
+        if Exponent >= 1 then
+            Result := Cartesian (Left);
+            while Exponent > 1 loop
+                Result := Result * Left;
+                Exponent := Exponent - 1;
+            end loop;
+        else
+            Result := One;
+            while Exponent < 0 loop
+                Result := Result / Left;
+                Exponent := Exponent + 1;
+            end loop;
+        end if;
+        return Result;
+    end "**";
+
+
+
+    function "<"
+           (Left, Right : in Imaginary)
+        return Boolean is
+    begin
+        return Real'Base (Left) < Real'Base (Right);
+    end "<";
+
+    function "<="
+           (Left, Right : in Imaginary)
+        return Boolean is
+    begin
+        return Real'Base (Left) <= Real'Base (Right);
+    end "<=";
+
+    function ">"
+           (Left, Right : in Imaginary)
+        return Boolean is
+    begin
+        return Real'Base (Left) > Real'Base (Right);
+    end ">";
+
+    function ">="
+           (Left, Right : in Imaginary)
+        return Boolean is
+    begin
+        return Real'Base (Left) >= Real'Base (Right);
+    end ">=";
+
+
+
+    function "+"
+           (Left  : in Complex;
+            Right : in Real'Base)
+        return Complex is
+    begin
+        return (Re => Left.Re + Right, Im => Left.Im);
+    end "+";
+
+    function "+"
+           (Left  : in Real'Base;
+            Right : in Complex)
+        return Complex is
+    begin
+        return Right + Left;
+    end "+";
+
+    function "-"
+           (Left  : in Complex;
+            Right : in Real'Base)
+        return Complex is
+    begin
+        return (Re => Left.Re - Right, Im => Left.Im);
+    end "-";
+
+    function "-"
+           (Left  : in Real'Base;
+            Right : in Complex)
+        return Complex is
+    begin
+        return (Re => Left - Right.Re, Im => -Right.Im);
+    end "-";
+
+    function "*"
+           (Left  : in Complex;
+            Right : in Real'Base)
+        return Complex is
+    begin
+        return (Re => Left.Re * Right, Im => Left.Im * Right);
+    end "*";
+
+    function "*"
+           (Left  : in Real'Base;
+            Right : in Complex)
+        return Complex is
+    begin
+        return Right * Left;
+    end "*";
+
+    function "/"
+           (Left  : in Complex;
+            Right : in Real'Base)
+        return Complex is
+    begin
+        return (Re => Left.Re / Right, Im => Left.Im / Right);
+    end "/";
+
+    function "/"
+           (Left  : in Real'Base;
+            Right : in Complex)
+        return Complex is
+    begin
+        return Cartesian (Left) / Right;
+    end "/";
+
+
+
+    function "+"
+           (Left  : in Complex;
+            Right : in Imaginary)
+        return Complex is
+    begin
+        return (Re => Left.Re, Im => Left.Im + Real'Base (Right));
+    end "+";
+
+    function "+"
+           (Left  : in Imaginary;
+            Right : in Complex)
+        return Complex is
+    begin
+        return Right + Left;
+    end "+";
+
+    function "-"
+           (Left  : in Complex;
+            Right : in Imaginary)
+        return Complex is
+    begin
+        return (Re => Left.Re, Im => Left.Im - Real'Base (Right));
+    end "-";
+
+    function "-"
+           (Left  : in Imaginary;
+            Right : in Complex)
+        return Complex is
+    begin
+        return (Re => -Right.Re, Im => Real'Base (Left) - Right.Im);
+    end "-";
+
+    function "*"
+           (Left  : in Complex;
+            Right : in Imaginary)
+        return Complex is
+    begin
+        return (Re => -Left.Im * Real'Base (Right),
+                Im => Left.Re * Real'Base (Right));
+    end "*";
+
+    function "*"
+           (Left  : in Imaginary;
+            Right : in Complex)
+        return Complex is
+    begin
+        return Right * Left;
+    end "*";
+
+    function "/"
+           (Left  : in Complex;
+            Right : in Imaginary)
+        return Complex is
+    begin
+        return (Re => Left.Im / Real'Base (Right),
+                Im => -Left.Re / Real'Base (Right));
+    end "/";
+
+    function "/"
+           (Left  : in Imaginary;
+            Right : in Complex)
+        return Complex is
+    begin
+        return Cartesian (Left) / Right;
+    end "/";
+
+
+
+    function "+"
+           (Left  : in Imaginary;
+            Right : in Real'Base)
+        return Complex is
+    begin
+        return (Re => Right, Im => Real'Base (Left));
+    end "+";
+
+    function "+"
+           (Left  : in Real'Base;
+            Right : in Imaginary)
+        return Complex is
+    begin
+        return (Re => Left, Im => Real'Base (Right));
+    end "+";
+
+    function "-"
+           (Left  : in Imaginary;
+            Right : in Real'Base)
+        return Complex is
+    begin
+        return (Re => -Right, Im => Real'Base (Left));
+    end "-";
+
+    function "-"
+           (Left  : in Real'Base;
+            Right : in Imaginary)
+        return Complex is
+    begin
+        return (Re => Left, Im => -Real'Base (Right));
+    end "-";
+
+    function "*"
+           (Left  : in Imaginary;
+            Right : in Real'Base)
+        return Imaginary is
+    begin
+        return Imaginary (Standard."*" (Right, Real'Base (Left)));
+    end "*";
+
+    function "*"
+           (Left  : in Real'Base;
+            Right : in Imaginary)
+        return Imaginary is
+    begin
+        return Imaginary (Standard."*" (Left, Real'Base (Right)));
+    end "*";
+
+    function "/"
+           (Left  : in Imaginary;
+            Right : in Real'Base)
+        return Imaginary is
+    begin
+        return Imaginary (Standard."/" (Real'Base (Left), Right));
+    end "/";
+
+    function "/"
+           (Left  : in Real'Base;
+            Right : in Imaginary)
+        return Imaginary is
+    begin
+        return -Imaginary (Standard."/" (Left, Real'Base (Right)));
+    end "/";
+
+end Complex_Fixed_Points;
+
diff --git a/src/complex_fixed_points.ads b/src/complex_fixed_points.ads
new file mode 100644
index 0000000..97667b1
--- /dev/null
+++ b/src/complex_fixed_points.ads
@@ -0,0 +1,117 @@
+
+generic
+    type Real is delta <> digits <>;
+package Complex_Fixed_Points is
+
+    pragma Pure (Complex_Fixed_Points);
+
+    type Complex is private;
+    type Imaginary is private;
+    pragma Preelaborable_Initialization (Imaginary);
+
+    i : constant Imaginary;
+    j : constant Imaginary;
+
+    Zero : constant Complex;
+    One  : constant Complex;
+
+    function Re (X : in Complex) return Real'Base;
+    function Im (X : in Complex) return Real'Base;
+    function Im (X : in Imaginary) return Real'Base;
+
+    procedure Set_Re
+           (X  : in out Complex;
+            Re : in     Real'Base);
+    procedure Set_Im
+           (X  : in out Complex;
+            Im : in     Real'Base);
+    procedure Set_Im
+           (X  : in out Imaginary;
+            Im : in     Real'Base);
+
+    function Cartesian (Re, Im : in Real'Base) return Complex;
+    function Cartesian (Re     : in Real'Base) return Complex;
+    function Cartesian (Im     : in Imaginary) return Complex;
+
+    function Modulus (X     : in Complex) return Real'Base;
+    function "abs"   (Right : in Complex) return Real'Base renames Modulus;
+
+    --  Like hell am I writing fixed point trig functions right now.
+
+    --  function Argument (X : in Complex) return Real'Base;
+    --  function Argument (X : in Complex; Cycle : in Real'Base) return Real'Base;
+
+    --  function Polar (Modulus, Argument        : in Real'Base) return Complex;
+    --  function Polar (Modulus, Argument, Cycle : in Real'Base) return Complex;
+
+    function "+" (Right : in Complex) return Complex;
+    function "-" (Right : in Complex) return Complex;
+    function Conjugate (X : in Complex) return Complex;
+
+    function "+" (Left, Right : in Complex) return Complex;
+    function "-" (Left, Right : in Complex) return Complex;
+    function "*" (Left, Right : in Complex) return Complex;
+    function "/" (Left, Right : in Complex) return Complex;
+
+    function "**" (Left : in Complex; Right : in Integer) return Complex;
+
+    function "+"       (Right : in Imaginary) return Imaginary;
+    function "-"       (Right : in Imaginary) return Imaginary;
+    function Conjugate (X     : in Imaginary) return Imaginary renames "-";
+    function "abs"     (Right : in Imaginary) return Real'Base;
+
+    function "+" (Left, Right : in Imaginary) return Imaginary;
+    function "-" (Left, Right : in Imaginary) return Imaginary;
+    function "*" (Left, Right : in Imaginary) return Real'Base;
+    function "/" (Left, Right : in Imaginary) return Real'Base;
+
+    function "**" (Left : in Imaginary; Right : in Integer) return Complex;
+
+    function "<"  (Left, Right : in Imaginary) return Boolean;
+    function "<=" (Left, Right : in Imaginary) return Boolean;
+    function ">"  (Left, Right : in Imaginary) return Boolean;
+    function ">=" (Left, Right : in Imaginary) return Boolean;
+
+    function "+" (Left : in Complex;   Right : in Real'Base) return Complex;
+    function "+" (Left : in Real'Base; Right : in Complex)   return Complex;
+    function "-" (Left : in Complex;   Right : in Real'Base) return Complex;
+    function "-" (Left : in Real'Base; Right : in Complex)   return Complex;
+    function "*" (Left : in Complex;   Right : in Real'Base) return Complex;
+    function "*" (Left : in Real'Base; Right : in Complex)   return Complex;
+    function "/" (Left : in Complex;   Right : in Real'Base) return Complex;
+    function "/" (Left : in Real'Base; Right : in Complex)   return Complex;
+
+    function "+" (Left : in Complex;   Right : in Imaginary) return Complex;
+    function "+" (Left : in Imaginary; Right : in Complex)   return Complex;
+    function "-" (Left : in Complex;   Right : in Imaginary) return Complex;
+    function "-" (Left : in Imaginary; Right : in Complex)   return Complex;
+    function "*" (Left : in Complex;   Right : in Imaginary) return Complex;
+    function "*" (Left : in Imaginary; Right : in Complex)   return Complex;
+    function "/" (Left : in Complex;   Right : in Imaginary) return Complex;
+    function "/" (Left : in Imaginary; Right : in Complex)   return Complex;
+
+    function "+" (Left : in Imaginary; Right : in Real'Base) return Complex;
+    function "+" (Left : in Real'Base; Right : in Imaginary) return Complex;
+    function "-" (Left : in Imaginary; Right : in Real'Base) return Complex;
+    function "-" (Left : in Real'Base; Right : in Imaginary) return Complex;
+    function "*" (Left : in Imaginary; Right : in Real'Base) return Imaginary;
+    function "*" (Left : in Real'Base; Right : in Imaginary) return Imaginary;
+    function "/" (Left : in Imaginary; Right : in Real'Base) return Imaginary;
+    function "/" (Left : in Real'Base; Right : in Imaginary) return Imaginary;
+
+private
+
+    type Complex is record
+        Re, Im : Real'Base;
+    end record;
+
+    type Imaginary is new Real'Base;
+
+    i : constant Imaginary := 1.0;
+    j : constant Imaginary := 1.0;
+
+    Zero : constant Complex := (Re => 0.0, Im => 0.0);
+    One  : constant Complex := (Re => 1.0, Im => 0.0);
+
+end Complex_Fixed_Points;
+
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;
+