From 0d8f4e413ba7a1f64753d73b757d7cc45e363369 Mon Sep 17 00:00:00 2001 From: Jedidiah Barber Date: Wed, 16 Nov 2022 01:06:13 +1300 Subject: Initial attempt --- src/complex_fixed_points.adb | 529 +++++++++++++++++++++++++++++++++++++++++++ src/complex_fixed_points.ads | 117 ++++++++++ src/fluid_simulator.adb | 197 ++++++++++++++++ 3 files changed, 843 insertions(+) create mode 100644 src/complex_fixed_points.adb create mode 100644 src/complex_fixed_points.ads create mode 100644 src/fluid_simulator.adb (limited to 'src') 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; + -- cgit