{$MODE OBJFPC} { -*- delphi -*- }
{$INCLUDE settings.inc}
unit gravity;

interface

uses
   rpc, storable, autostorable, stringstream, actors, spacemanagers, callbacks;

//{$DEFINE GRAVITY_VERBOSE}

type
   TGravitationalSystem = class(TSpaceManager) // @RegisterActorClass
    protected
      type
       PSystemParticipant = ^TSystemParticipant;
       TSystemParticipant = record
         PosX, PosY: TGravitationalSystemValue;
         VelX, VelY: TGravitationalSystemValue;
         Actor: TChildActor;
         Next: PSystemParticipant;
       end;
       TPointType = (ptStart, ptEnd);
       TCollisionDetectionDatum = record
         Position: TGravitationalSystemValue;
         PointType: TPointType;
         ParticipantIndex: Cardinal; // index into SystemData during tick (same order as FParticipants)
       end;
       TCollisionDetectionData = array of TCollisionDetectionDatum;
      var
       FParticipants: PSystemParticipant;
       FCachedParticipantCount: Cardinal;
       FCurrentTickInterval: TGravitationalSystemValue;
       FCurrentRadius: TGravitationalSystemValue;
       FCollisionFilterListX, FCollisionFilterListY: TCollisionDetectionData;
      function GetRendererClass(): UTF8String; override;
      function GetRadius(): TGravitationalSystemValue; override;
      procedure InternalAddChild(const NewRecord: PSystemParticipant);
      {$IFDEF DEBUG} function GetNameByIndex(Index: Cardinal): UTF8String; {$ENDIF}
      {$IFDEF DEBUG} procedure DebugPrintFilterLists(); {$ENDIF}
    public
      constructor Create(AName: UTF8String); overload; deprecated 'For hacky genesis only';
      constructor Create(Dynasty: TAbstractDynasty; Settings: TAutoStorable); override; overload;
      destructor Destroy(); override;
      constructor Read(Stream: TReadStream); override;
      procedure Write(Stream: TWriteStream); override;
      procedure ReportState(Dynasty: TAbstractDynasty; Stream: TStringStreamWriter); override;
      procedure Tick(Interval: TDateTime); override;
      function ParseLocationFromMessage(Stream: TStringStreamReader): TAbstractLocation; override;
      procedure RemoveChild(const OldChild: TAbstractTreeNode); override;
      function AddChild(const NewChild: TAbstractTreeNode; const Location: TAbstractLocation): Boolean; override;
      procedure KillChildren(); override;
      procedure UpdatePhysicalProperties(); override;
      procedure TeleportToOrbit(Actor: TChildActor; Orbit: TGravitationalSystemValue; Clockwise: Boolean; Around: TChildActor);
      procedure Launch(Actor: TChildActor; Around: TChildActor; Angle: TGravitationalSystemValue); deprecated 'This should be obsolete now we have Thrust'; 
  end;

implementation

uses
   math, exceptions;

const
   MinimumTickInterval = 1.0*60;
   DefaultMaximumTickInterval = 1.0*60; // in seconds (this is absurdly small but ensures we start off conservatively)

constructor TGravitationalSystem.Create(AName: UTF8String);
begin
   inherited Create(AName);
   FCurrentTickInterval := DefaultMaximumTickInterval;
end;

constructor TGravitationalSystem.Create(Dynasty: TAbstractDynasty; Settings: TAutoStorable);
begin
   inherited;
   FCurrentTickInterval := DefaultMaximumTickInterval;
end;

destructor TGravitationalSystem.Destroy();
var
   Participant: PSystemParticipant;
begin
   while (Assigned(FParticipants)) do
   begin
      Participant := FParticipants;
      FParticipants := FParticipants^.Next;
      Participant^.Actor.Free();
      Dispose(Participant);
   end;
   FCachedParticipantCount := 0;
   inherited;
end;

constructor TGravitationalSystem.Read(Stream: TReadStream);
var
   Participant: PSystemParticipant;
   Actor: TChildActor;

   function GetNextActor(): Boolean;
   begin
      Actor := Stream.ReadObject() as TChildActor;
      Result := Assigned(Actor);
   end;

begin
   inherited;
   while (GetNextActor()) do
   begin
      New(Participant);
      Participant^.Actor := Actor; {BOGUS Warning: Local variable "Actor" does not seem to be initialized}
      Participant^.PosX := Stream.ReadDouble();
      Participant^.PosY := Stream.ReadDouble();
      Participant^.VelX := Stream.ReadDouble();
      Participant^.VelY := Stream.ReadDouble();
      Participant^.Next := FParticipants;
      FParticipants := Participant;
      Inc(FCachedParticipantCount);
   end;
   FCurrentTickInterval := Stream.ReadDouble();
end;

procedure TGravitationalSystem.Write(Stream: TWriteStream);
var
   Participant: PSystemParticipant;
begin
   inherited;
   Participant := FParticipants;
   while (Assigned(Participant)) do
   begin
      Stream.WriteObject(Participant^.Actor);
      Stream.WriteDouble(Participant^.PosX);
      Stream.WriteDouble(Participant^.PosY);
      Stream.WriteDouble(Participant^.VelX);
      Stream.WriteDouble(Participant^.VelY);
      Participant := Participant^.Next;
   end;
   Stream.WriteObject(nil);
   Stream.WriteDouble(FCurrentTickInterval);
end;

function TGravitationalSystem.GetRadius(): TGravitationalSystemValue;
begin
   Result := FCurrentRadius;
end;

const
   G = 6.67384E-11; // Gravitational Constant (m^3 kg^-1 s^-2)

procedure TGravitationalSystem.Tick(Interval: TDateTime);

   procedure ApplyThrust();
   begin
      XXX;
   end;

type
   TSystemStateEntry = record
      X, Y, VX, VY: TGravitationalSystemValue;
   end;
   TSystemState = array of TSystemStateEntry;
   TMasses = array of TGravitationalSystemValue;
   TRadii = array of TGravitationalSystemValue;
   TAccelerationEntry = record
      X, Y: TGravitationalSystemValue;
   end;
   TAccelerations = array of TAccelerationEntry;
   
   procedure GetSystemData(out State: TSystemState; out Masses: TMasses; out Radii: TRadii);
   var
      Index: Cardinal;
      Participant: PSystemParticipant;
   begin
      SetLength(State, FCachedParticipantCount);
      SetLength(Masses, FCachedParticipantCount);
      SetLength(Radii, FCachedParticipantCount);
      Participant := FParticipants;
      Assert(Length(State) > 0);
      for Index := Low(State) to High(State) do
      begin
         Assert(Assigned(Participant));
         Assert(Assigned(Participant^.Actor));
         State[Index].X := Participant^.PosX;
         State[Index].Y := Participant^.PosY;
         State[Index].VX := Participant^.VelX;
         State[Index].VY := Participant^.VelY;
         Masses[Index] := Participant^.Actor.Mass;
         Radii[Index] := Participant^.Actor.Radius;
         Participant := Participant^.Next;
      end;
   end;

   function ComputeAccelerations(const SystemData: TSystemState; const Masses: TMasses): TAccelerations; // returns accelerations
   var
      Index, SubIndex: Cardinal;
      DX, DY, R, FoverR, FX, FY: TGravitationalSystemValue;
   begin
      Assert(Length(SystemData) = FCachedParticipantCount);
      SetLength(Result, FCachedParticipantCount);
      Assert(Length(SystemData) > 0);
      for Index := Low(SystemData) to High(SystemData) do
      begin
         if (Index < High(SystemData)) then
            for SubIndex := Index+1 to High(SystemData) do
            begin
               DX := SystemData[SubIndex].X - SystemData[Index].X;
               DY := SystemData[SubIndex].Y - SystemData[Index].Y;
               R := SqRt(DX*DX+DY*DY);
               Assert(R > 0);
               FoverR := G * Masses[Index] * Masses[SubIndex] / (R*R*R); // F = G*m1*m2/(R^2), but we predivide by R to avoid having to do it twice below
               FX := FoverR * DX; // FX = F*DX/R
               Result[Index].X := Result[Index].X + FX;
               Result[SubIndex].X := Result[SubIndex].X - FX;
               FY := FoverR * DY; // FY = F*DY/R
               Result[Index].Y := Result[Index].Y + FY;
               Result[SubIndex].Y := Result[SubIndex].Y - FY;
            end;
         Assert(Masses[Index] > 0, 'Massless particles are not yet supported in the gravitational simulation');
         Result[Index].X := Result[Index].X / Masses[Index];
         Result[Index].Y := Result[Index].Y / Masses[Index];
      end;
   end;

   function ApplyAcceleration(const Accelerations: TAccelerations; const DT: TGravitationalSystemValue; const SystemData: TSystemState): TSystemState;
   var
      Index: Cardinal;
      DTSquaredOverTwo: TGravitationalSystemValue;
   begin
      Assert(Length(SystemData) = FCachedParticipantCount);
      Assert(Length(Accelerations) = FCachedParticipantCount);
      SetLength(Result, FCachedParticipantCount);
      DTSquaredOverTwo := DT*DT/2.0;
      Assert(Length(Result) > 0);
      for Index := Low(Result) to High(Result) do
      begin
         Result[Index].X  := Accelerations[Index].X * DTSquaredOverTwo + SystemData[Index].VX * DT + SystemData[Index].X;
         Result[Index].VX := Accelerations[Index].X * DT               + SystemData[Index].VX;
         Result[Index].Y  := Accelerations[Index].Y * DTSquaredOverTwo + SystemData[Index].VY * DT + SystemData[Index].Y;
         Result[Index].VY := Accelerations[Index].Y * DT               + SystemData[Index].VY;
      end;
   end;

type
   PCollision = ^TCollision;
   TCollision = record
      ParticipantA, ParticipantB: Cardinal;
      Next: PCollision;
   end;

   procedure CocktailRangeSort(var List: TCollisionDetectionData);

      function Compare(IndexA, IndexB: Cardinal): Boolean; inline;
      begin
         // Returns True if Value(IndexA) > Value(IndexB)
         Assert(ptStart < ptEnd);
         if (List[IndexA].ParticipantIndex = List[IndexB].ParticipantIndex) then
         begin
            // always make sure an object ends after it starts
            Result := List[IndexA].PointType > List[IndexB].PointType;
            // make sure the order is still sane
            Assert((Result and (List[IndexA].Position >= List[IndexB].Position)) xor ((not Result) and (List[IndexA].Position <= List[IndexB].Position)));
         end
         else
            Result := (List[IndexA].Position > List[IndexB].Position) or
                      (List[IndexA].PointType < List[IndexB].PointType); // sort ptEnd before ptStart, so that touching isn't a collision
      end;

      procedure Swap(IndexA, IndexB: Cardinal);
      var
         Temp: TCollisionDetectionDatum;
      begin
         Temp := List[IndexA];
         List[IndexA] := List[IndexB];
         List[IndexB] := Temp;
      end;

   var
      Left, Right, Index: Cardinal;
      FlipFlop: Boolean;
   begin
      Assert(Length(List) > 1);
      Left := Low(List);
      Right := High(List);
      Assert(Left < Right);
      FlipFlop := False;
      repeat
         Assert(not FlipFlop); // In this first section, FlipFlop becomes true if we swapped anything
         Assert(Left < Right);
         for Index := Left to Right-1 do
         begin
            if (Compare(Index, Index+1)) then
            begin
               Swap(Index, Index+1);
               FlipFlop := True;
            end;
         end;
         Inc(Left);
         if ((not FlipFlop) or (Left >= Right)) then // if still false, didn't swap anything, abort; also abort if we're out of things to sort
            Break;
         Assert(FlipFlop); // In this section section, FlipFlop becomes false if we swapped anything
         Assert(Left < Right);
         for Index := Right downto Left+1 do
         begin
            if (Compare(Index-1, Index)) then
            begin
               Swap(Index-1, Index);
               FlipFlop := False;
            end;
         end;
         Dec(Right);
      until FlipFlop; // if still true, didn't swap anything, abort
   end;

   function DetectCollisions(const SystemData: TSystemState; const Radii: TRadii; out Collisions: PCollision): Boolean;

      procedure UpdateCollisionDetectionData(); inline;
      var
         Index: Cardinal;
      begin
         for Index := Low(FCollisionFilterListX) to High(FCollisionFilterListX) do
         begin
            case (FCollisionFilterListX[Index].PointType) of
             ptStart: FCollisionFilterListX[Index].Position := SystemData[FCollisionFilterListX[Index].ParticipantIndex].X
                                                                       - Radii[FCollisionFilterListX[Index].ParticipantIndex];
             ptEnd: FCollisionFilterListX[Index].Position := SystemData[FCollisionFilterListX[Index].ParticipantIndex].X
                                                                     + Radii[FCollisionFilterListX[Index].ParticipantIndex];
             else Assert(False, 'Bogus PointType');
            end;
         end;
         for Index := Low(FCollisionFilterListY) to High(FCollisionFilterListY) do
         begin
            case (FCollisionFilterListY[Index].PointType) of
             ptStart: FCollisionFilterListY[Index].Position := SystemData[FCollisionFilterListY[Index].ParticipantIndex].Y
                                                                       - Radii[FCollisionFilterListY[Index].ParticipantIndex];
             ptEnd: FCollisionFilterListY[Index].Position := SystemData[FCollisionFilterListY[Index].ParticipantIndex].Y
                                                                     + Radii[FCollisionFilterListY[Index].ParticipantIndex];
             else Assert(False, 'Bogus PointType');
            end;
         end;
      end;


      // Consider this setup when working on these algorithms:
      //
      //  A |--------------------------|
      //  B                   |----------|
      //  C                     |----------|
      //  D                       |---------------------------------------|
      //  E                        |-|
      //  F                               |-------------------------------------------------------|
      //
      // or:
      //
      //  A    |-------------------|
      //  B                      |----------|
      //  C                              |-------------------|

   const
      kDeleted = High(Cardinal);
   var
      PotentialCollisions, CurrentCollision, NewCollision: PCollision;
      LastPotentialCollisionSplicePoint: ^PCollision;
      Index, SubIndex, IndexA, IndexB, OpenCount, MaxCount: Cardinal;
      OpenParticipantsStack: array of Cardinal;
      ParticipantToStackMap: array of Cardinal;
   begin
      Assert(Length(FCollisionFilterListX) = Length(FCollisionFilterListY));
      Assert(Length(FCollisionFilterListX) = Length(SystemData) * 2);
      Assert(Length(SystemData) = Length(Radii));
      Collisions := nil;
      if (Length(SystemData) < 2) then
      begin
         Result := False;
         Exit; // collisions obviously require at least two objects
      end;
      Assert(Length(FCollisionFilterListX) > 2);
      UpdateCollisionDetectionData();
      CocktailRangeSort(FCollisionFilterListX);
      CocktailRangeSort(FCollisionFilterListY);
      // First, look for potential collisions in the x axis
      PotentialCollisions := nil;
      OpenCount := 0;
      MaxCount := 0;
      // there should be a better way of doing this, even if it's just keeping them around across ticks...
      SetLength(OpenParticipantsStack, Length(SystemData));
      SetLength(ParticipantToStackMap, Length(SystemData));
      // The way this works is we walk the list of open/close points counting how many opens we've seen that haven't had a matching close.
      // (We assume that you'll never get a close before an open, so any close we see is always implicitly matching one of the pending opens.)
      // When we see an open, we add it to the list of open ones, OpenParticipantsStack, and increase our OpenCount.
      // That array points to the entry in the SystemData lists.
      // When we see a close, we wipe that entry from the OpenParticipantsStack by setting it to kDeleted.
      // Since that leaves holes, we actually also track which slot we're using in the OpenParticipantsStack using MaxCount
      // When we delete an entry from the stack, we check if that entry was the last one.
      // If it was, we reduce the MaxCount so that all the trailing kDeleted items are off our stack.
      // WARNING: There's a lot of code duplication here; when fixing bugs here, see also below
      for Index := Low(FCollisionFilterListX) to High(FCollisionFilterListX) do
      begin
         case (FCollisionFilterListX[Index].PointType) of
          ptStart:
             begin
                Assert(OpenCount <= MaxCount);
                Assert(MaxCount < Length(SystemData));
                OpenParticipantsStack[MaxCount] := FCollisionFilterListX[Index].ParticipantIndex;
                ParticipantToStackMap[FCollisionFilterListX[Index].ParticipantIndex] := MaxCount;
                if (OpenCount > 0) then
                begin
                   // new potential collision detected
                   // collide with all the open actors
                   for SubIndex := Low(OpenParticipantsStack) to MaxCount-1 do
                   begin
                      if (OpenParticipantsStack[SubIndex] <> kDeleted) then
                      begin
                         New(NewCollision);
                         Assert(FCollisionFilterListX[Index].ParticipantIndex <> OpenParticipantsStack[SubIndex]);
                         if (FCollisionFilterListX[Index].ParticipantIndex < OpenParticipantsStack[SubIndex]) then
                         begin
                            NewCollision^.ParticipantA := FCollisionFilterListX[Index].ParticipantIndex;
                            NewCollision^.ParticipantB := OpenParticipantsStack[SubIndex];
                         end
                         else
                         begin
                            NewCollision^.ParticipantA := OpenParticipantsStack[SubIndex];
                            NewCollision^.ParticipantB := FCollisionFilterListX[Index].ParticipantIndex;
                         end;
                         NewCollision^.Next := PotentialCollisions;
                         PotentialCollisions := NewCollision;
                      end;
                   end;
                end;
                Inc(OpenCount);
                Inc(MaxCount);
                Assert(OpenCount <= MaxCount);
                Assert(MaxCount <= Length(SystemData));
             end;
          ptEnd:
             begin
                Assert(OpenCount > 0);
                Dec(OpenCount);
                if (ParticipantToStackMap[FCollisionFilterListX[Index].ParticipantIndex] = MaxCount-1) then
                begin
                   repeat
                      Dec(MaxCount);
                   until ((MaxCount = 0) or (OpenParticipantsStack[MaxCount-1] <> kDeleted));
                end
                else
                begin
                   OpenParticipantsStack[ParticipantToStackMap[FCollisionFilterListX[Index].ParticipantIndex]] := kDeleted;
                end;
             end;
          else Assert(False, 'Bogus PointType');
         end;
      end;      
      if (Assigned(PotentialCollisions)) then
      begin
         // Now we walk the FCollisionFilterListY list looking for potential clashes, when there is one, check if it's a known one, if so, note it in TaggedCollisions
         // WARNING: There's a lot of code duplication here; when fixing bugs here, see also above
         OpenCount := 0;
         MaxCount := 0;
         // This is similar to the case above for X, except that when we get a potential clash, we look for it in the array and then
         // put it in the TaggedCollisions list instead
         for Index := Low(FCollisionFilterListY) to High(FCollisionFilterListY) do
         begin
            case (FCollisionFilterListY[Index].PointType) of
             ptStart:
                begin
                   Assert(OpenCount <= MaxCount);
                   Assert(MaxCount < Length(SystemData));
                   OpenParticipantsStack[MaxCount] := FCollisionFilterListY[Index].ParticipantIndex;
                   ParticipantToStackMap[FCollisionFilterListY[Index].ParticipantIndex] := MaxCount;
                   if (OpenCount > 0) then
                   begin
                      // new potential collision detected
                      // for each collision, check if we have a match
                      for SubIndex := Low(OpenParticipantsStack) to MaxCount-1 do
                      begin
                         if (OpenParticipantsStack[SubIndex] <> kDeleted) then
                         begin
                            CurrentCollision := PotentialCollisions;
                            LastPotentialCollisionSplicePoint := @PotentialCollisions;
                            while (Assigned(CurrentCollision)) do
                            begin
                               Assert(FCollisionFilterListY[Index].ParticipantIndex <> OpenParticipantsStack[SubIndex]);
                               if (FCollisionFilterListY[Index].ParticipantIndex < OpenParticipantsStack[SubIndex]) then
                               begin
                                  IndexA := FCollisionFilterListY[Index].ParticipantIndex;
                                  IndexB := OpenParticipantsStack[SubIndex];
                               end
                               else
                               begin
                                  IndexA := OpenParticipantsStack[SubIndex];
                                  IndexB := FCollisionFilterListY[Index].ParticipantIndex;
                               end;
                               if (((CurrentCollision^.ParticipantA = IndexA) and (CurrentCollision^.ParticipantB = IndexB))) then
                               begin
                                  // remove this entry from the potential list, move it to the tagged list
                                  LastPotentialCollisionSplicePoint^ := CurrentCollision^.Next;
                                  CurrentCollision^.Next := Collisions;
                                  Collisions := CurrentCollision;
                                  Break;
                               end;
                               LastPotentialCollisionSplicePoint := @CurrentCollision^.Next;
                               CurrentCollision := CurrentCollision^.Next;
                            end;
                            if (not Assigned(PotentialCollisions)) then
                               Break;
                         end;
                      end;
                   end;
                   Inc(OpenCount);
                   Inc(MaxCount);
                   Assert(OpenCount <= MaxCount);
                   Assert(MaxCount <= Length(SystemData));
                end;
             ptEnd:
                begin
                   Assert(OpenCount > 0);
                   Dec(OpenCount);
                   if (ParticipantToStackMap[FCollisionFilterListY[Index].ParticipantIndex] = MaxCount-1) then
                   begin
                      repeat
                         Dec(MaxCount);
                      until ((MaxCount = 0) or (OpenParticipantsStack[MaxCount-1] <> kDeleted));
                   end
                   else
                   begin
                      OpenParticipantsStack[ParticipantToStackMap[FCollisionFilterListY[Index].ParticipantIndex]] := kDeleted;
                   end;
                end;
             else Assert(False, 'Bogus PointType');
            end;
            if (not Assigned(PotentialCollisions)) then
               Break;
         end;
      end;
      // Free the non-matching potentials
      while (Assigned(PotentialCollisions)) do
      begin
         CurrentCollision := PotentialCollisions;
         PotentialCollisions := PotentialCollisions^.Next;
         Dispose(CurrentCollision);
      end;
      Result := Assigned(Collisions);
   end;

   procedure ApplyCollisions(var Collisions: PCollision);
   var
      CurrentCollision, SubCollision, TempCollision: PCollision;
      SubCollisionSplicePoint: ^PCollision;
      Participant: PSystemParticipant;
      ActorA, ActorB, TempActor: TChildActor;
      CrashArguments: TCrashArguments;
      Index, VictimIndex: Cardinal;
   begin
      Assert(Assigned(Collisions));
      CurrentCollision := Collisions;
      while (Assigned(CurrentCollision)) do
         CurrentCollision := CurrentCollision^.Next;
      while (Assigned(Collisions)) do
      begin
         CurrentCollision := Collisions;
         Collisions := CurrentCollision^.Next;
         Participant := FParticipants;
         Index := 0;
         while (Index < CurrentCollision^.ParticipantA) do
         begin
            Assert(Assigned(Participant));
            Participant := Participant^.Next;
            Inc(Index);
         end;
         Assert(Assigned(Participant));
         ActorA := Participant^.Actor;
         Participant := Participant^.Next;
         Inc(Index);
         while (Index < CurrentCollision^.ParticipantB) do
         begin
            Assert(Assigned(Participant));
            Participant := Participant^.Next;
            Inc(Index);
         end;
         Assert(Assigned(Participant));
         ActorB := Participant^.Actor;
         if (ActorA.Mass < ActorB.Mass) then
         begin
            TempActor := ActorA;
            ActorA := ActorB;
            ActorB := TempActor;
            VictimIndex := CurrentCollision^.ParticipantA;
         end
         else
         begin
            VictimIndex := CurrentCollision^.ParticipantB;
         end;
         CrashArguments.Signature := ctCrash;
         CrashArguments.Victim := ActorB;
         CrashArguments.Destination := ActorA;
         ActorA.RunCallbacks(CrashArguments);
         Assert(Assigned(CrashArguments.Destination));
         Assert(CrashArguments.Destination is TPhysicalActor);
         // XXX need to get the location of the crash and pass it as the second argument here
         (CrashArguments.Destination as TPhysicalActor).Crash(ActorB, nil);
         // Now remove any collisions that collided with the now departed actor
         SubCollisionSplicePoint := @Collisions;
         SubCollision := Collisions;
         while (Assigned(SubCollision)) do
         begin
            if ((SubCollision^.ParticipantA = CurrentCollision^.ParticipantB) or
                (SubCollision^.ParticipantB = CurrentCollision^.ParticipantB)) then
            begin
               SubCollisionSplicePoint^ := SubCollision^.Next;
               TempCollision := SubCollision;
               SubCollision := SubCollision^.Next;
               Dispose(TempCollision);
            end
            else
            begin
               if (SubCollision^.ParticipantA > VictimIndex) then
                  Dec(SubCollision^.ParticipantA);
               if (SubCollision^.ParticipantB > VictimIndex) then
                  Dec(SubCollision^.ParticipantB);
               SubCollisionSplicePoint := @SubCollision^.Next;
               SubCollision := SubCollision^.Next;
            end;
         end;
         Dispose(CurrentCollision);
      end;
   end;

var
   BiggestErrorFactor: TGravitationalSystemValue;

   procedure ConsiderStepError(RKF5, RKF4, DT, ToleratedError: TGravitationalSystemValue); inline;
   var
      ErrorFactor: TGravitationalSystemValue;
   begin
      if (RKF4 <> RKF5) then
      begin
         ErrorFactor := ToleratedError*DT / Abs(RKF4-RKF5);
         if (ErrorFactor < BiggestErrorFactor) then
            BiggestErrorFactor := ErrorFactor;
      end;
   end;

const // http://mathfaculty.fullerton.edu/mathews/n2003/RungeKuttaFehlbergMod.html (Runge-Kutta-Fehlberg method, RKF54)
      // XXX switch to http://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method table
   A2a = 1.0/4.0;
   A3a = 3.0/32.0;
   A3b = 9.0/32.0;
   A4a = 1932.0/2197.0;
   A4b = -7200.0/2197.0;
   A4c = 7296.0/2197.0;
   A5a = 439.0/216.0;
   A5b = -8.0;
   A5c = 3680.0/513.0;
   A5d = -845.0/4104.0;
   A6a = -8.0/27.0;
   A6b = 2.0;
   A6c = -3544.0/2565.0;
   A6d = 1859.0/4104.0;
   A6e = -11.0/40.0;
   C4_A1 = 25.0/216.0;
   C4_A3 = 1408.0/2565.0;
   C4_A4 = 2197.0/4104.0;
   C4_A5 = -1.0/5.0;
   C5_A1 = 16.0/135.0;
   C5_A3 = 6656.0/12825.0;
   C5_A4 = 28561.0/56430.0;
   C5_A5 = -9.0/50.0;
   C5_A6 = 2.0/55.0;

const
   FourthRootOfTwo = SqRt(SqRt(2));
   // XXX maybe make the error tolerances relative to the object radius?
   ToleratedErrorPos = 1000.0/(1000.0*60.0); // max error: 1km per 1000 minute tick (in m/s)
   ToleratedErrorVel = 0.01/(1000.0*60.0); // max error: 1m/s per 1000 minute tick (in m/(s^2))
   TickShrinkingFactor = 0.8; // we tick smaller than we call ideal, so as to avoid just narrowly missing the ideal window

var
   InitialConditions, RKF4, RKF5: TSystemState;
   Masses: TMasses;
   Radii: TRadii;
   A1, A2, A3, A4, A5, A6: TAccelerations;
   DT, RemainingTimeToSimulate, RemainingTimeToTick, IdealTick, CandidateTick: TGravitationalSystemValue;
   TickTimeInterval: TDateTime;
   Index: Cardinal;
   Participant: PSystemParticipant;
   Collisions: PCollision;
begin
   inherited;
   if (FCachedParticipantCount > 0) then
   begin
      RemainingTimeToSimulate := Interval*kSecondsPerDay; // days to seconds conversion
      RemainingTimeToTick := RemainingTimeToSimulate;
      Assert(Assigned(FParticipants));
      ApplyThrust();
      {$IFDEF GRAVITY_VERBOSE} Writeln('N-body simulation: ', (RemainingTimeToSimulate/60.0):5:2, ' minutes need ticking in system with ', FCachedParticipantCount, ' bodies; max tick interval is ', (FCurrentTickInterval/60.0):5:2, ' minutes per tick.'); {$ENDIF}
      // This computes the Runge-Kutta-Fehlberg values of K1-K6,
      // except instead of working out the actual K1-K6 values, we work
      // out the corresponding accelerations and apply those.
      // (Because for position, the K values here would actually be
      // functions of DT, so there's no single value we can calculate.)
      // To reduce confusion, I've used A1-A6 instead of K1-K6.
      GetSystemData(InitialConditions, Masses, Radii);
      A1 := ComputeAccelerations(InitialConditions, Masses);
      repeat
         Assert(FCurrentTickInterval >= MinimumTickInterval);
         Assert(RemainingTimeToSimulate > 0);
         if (Length(InitialConditions) > 1) then
         begin
            if (RemainingTimeToSimulate > FCurrentTickInterval) then
               DT := RemainingTimeToSimulate / Ceil(RemainingTimeToSimulate / FCurrentTickInterval)
            else
               DT := RemainingTimeToSimulate;
            A2 := ComputeAccelerations(ApplyAcceleration(A1, DT*A2a, InitialConditions), Masses);
            A3 := ComputeAccelerations(ApplyAcceleration(A1, DT*A3a,
                                       ApplyAcceleration(A2, DT*A3b, InitialConditions)), Masses);
            A4 := ComputeAccelerations(ApplyAcceleration(A1, DT*A4a,
                                       ApplyAcceleration(A2, DT*A4b,
                                       ApplyAcceleration(A3, DT*A4c, InitialConditions))), Masses);
            A5 := ComputeAccelerations(ApplyAcceleration(A1, DT*A5a,
                                       ApplyAcceleration(A2, DT*A5b,
                                       ApplyAcceleration(A3, DT*A5c,
                                       ApplyAcceleration(A4, DT*A5d, InitialConditions)))), Masses);
            A6 := ComputeAccelerations(ApplyAcceleration(A1, DT*A6a,
                                       ApplyAcceleration(A2, DT*A6b,
                                       ApplyAcceleration(A3, DT*A6c,
                                       ApplyAcceleration(A4, DT*A6d,
                                       ApplyAcceleration(A5, DT*A6e, InitialConditions))))), Masses);
            RKF4 := ApplyAcceleration(A1, DT*C4_A1,
                    ApplyAcceleration(A3, DT*C4_A3,
                    ApplyAcceleration(A4, DT*C4_A4,
                    ApplyAcceleration(A5, DT*C4_A5, InitialConditions))));
            RKF5 := ApplyAcceleration(A1, DT*C5_A1,
                    ApplyAcceleration(A3, DT*C5_A3,
                    ApplyAcceleration(A4, DT*C5_A4,
                    ApplyAcceleration(A5, DT*C5_A5,
                    ApplyAcceleration(A6, DT*C5_A6, InitialConditions)))));
            Assert(Length(RKF4) = Length(InitialConditions));
            Assert(Length(RKF5) = Length(InitialConditions));
            {$PUSH} {$IEEEERRORS OFF} BiggestErrorFactor := Infinity; {$POP}
            Assert(Length(InitialConditions) > 0);
            for Index := Low(InitialConditions) to High(InitialConditions) do
            begin
               ConsiderStepError(RKF5[Index].X, RKF4[Index].X, DT, ToleratedErrorPos);
               ConsiderStepError(RKF5[Index].Y, RKF4[Index].Y, DT, ToleratedErrorPos);
               ConsiderStepError(RKF5[Index].VX, RKF4[Index].VX, DT, ToleratedErrorVel);
               ConsiderStepError(RKF5[Index].VY, RKF4[Index].VY, DT, ToleratedErrorVel);
            end;
            if (not IsInfinite(BiggestErrorFactor)) then
            begin
               // IdealTick is how long we should have ticked at once, according to this...
               IdealTick := DT*FourthRootOfTwo*Sqrt(Sqrt(BiggestErrorFactor)); // ...dark magic
               CandidateTick := TickShrinkingFactor*IdealTick-2; // (but we tick a little less than that, just in case it's optimistic)
               if (CandidateTick < MinimumTickInterval) then
               begin
                  // However optimistic it is, this is getting out of hand. Let's just agree to disagree and
                  // aim for what is the smallest tick that it sanely makes sense to actually run at a time
                  IdealTick := MinimumTickInterval;
                  FCurrentTickInterval := MinimumTickInterval;
               end
               else
               begin
                  // The suggested tick length is reasonable, let's plan to do this next time
                  // (we check if we used it _this_ time lower down)
                  FCurrentTickInterval := CandidateTick;
               end;
            end
            else
            begin
               // if there's no conclusion to be drawn from the error data, assume that we
               // were doing fine last tick
               IdealTick := FCurrentTickInterval;
            end;
         end
         else
         begin
            // if there's only one thing to simulate, just do it all in one step and assume that's the right amount of time to do
            DT := RemainingTimeToSimulate;
            IdealTick := RemainingTimeToSimulate;
            RKF5 := ApplyAcceleration(A1, RemainingTimeToSimulate, InitialConditions);
         end;
         // if we simulated no more than the dark magic maths says we should have, then let's move on
         if (DT <= IdealTick) then
         begin
            InitialConditions := RKF5;
            RemainingTimeToSimulate := RemainingTimeToSimulate - DT;
            if (DetectCollisions(InitialConditions, Radii, Collisions) or (RemainingTimeToSimulate <= 0)) then
            begin
               TickTimeInterval := (RemainingTimeToTick - RemainingTimeToSimulate)/kSecondsPerDay;
               Participant := FParticipants;
               Assert(Length(InitialConditions) > 0);
               for Index := Low(InitialConditions) to High(InitialConditions) do
               begin
                  Assert(Assigned(Participant));
                  Participant^.PosX := InitialConditions[Index].X;
                  Participant^.PosY := InitialConditions[Index].Y;
                  Participant^.VelX := InitialConditions[Index].VX;
                  Participant^.VelY := InitialConditions[Index].VY;
                  Participant := Participant^.Next;
               end;
               Participant := FParticipants;
               Assert(Length(InitialConditions) > 0);
               for Index := Low(InitialConditions) to High(InitialConditions) do
               begin
                  Assert(Assigned(Participant));
                  Assert(Assigned(Participant^.Actor));
                  Participant^.Actor.Tick(TickTimeInterval);
                  Participant := Participant^.Next;
               end;
               RemainingTimeToTick := RemainingTimeToSimulate;
               if (Assigned(Collisions)) then
               begin
                  repeat
                     Assert(Assigned(Collisions));
                     ApplyCollisions(Collisions);
                     GetSystemData(InitialConditions, Masses, Radii);
                  until not DetectCollisions(InitialConditions, Radii, Collisions);
               end;
            end;
            Assert(not Assigned(Collisions));
            if (RemainingTimeToSimulate > 0) then
            begin
               // Prime the data structures for the next round
               A1 := ComputeAccelerations(InitialConditions, Masses);
            end;
         end
         else
         begin
            // we simulated too much, try again
            {$IFDEF GRAVITY_VERBOSE}
              Writeln('  Redoing step. Was using ', (DT/60.0):5:2, ' minutes per tick; calculated suggestion of ',
                                             (IdealTick/60.0):5:2, ' minutes per tick; now trying ',
                                  (FCurrentTickInterval/60.0):5:2, ' minutes per tick.');
            {$ENDIF}
         end;
      until RemainingTimeToSimulate <= 0;
      // XXX determine if we need to split into different solar systems
      // XXX recenter the system's center of mass
      // this will require some mechanism to tell the parent solar system, if any, what we're doing
      // (no need to propagate velocities, that's taken care of by moving our position)
      // XXX update our radius (if we're centered, that's the distance to object furthest from 0,0)
      MarkDirty([dfNeedNotifications]);
   end;
end;

procedure TGravitationalSystem.ReportState(Dynasty: TAbstractDynasty; Stream: TStringStreamWriter);
var
   Participant: PSystemParticipant;
begin
   inherited;
   Stream.WriteCardinal(FCachedParticipantCount);
   Participant := FParticipants;
   while (Assigned(Participant)) do
   begin
      Stream.WriteDouble(Participant^.PosX);
      Stream.WriteDouble(Participant^.PosY);
      Stream.WriteDouble(Participant^.VelX);
      Stream.WriteDouble(Participant^.VelY);
      Participant^.Actor.ReportIdentity(Stream);
      Participant := Participant^.Next;
   end;
end;

function TGravitationalSystem.ParseLocationFromMessage(Stream: TStringStreamReader): TAbstractLocation;
begin
   XXX;
   Result := nil;
end;

procedure TGravitationalSystem.RemoveChild(const OldChild: TAbstractTreeNode);

   procedure RemoveActorFromFilterList(var List: TCollisionDetectionData; const VictimIndex: Cardinal);
   var
      Index, FirstIndex, SecondIndex: Cardinal;
   begin
      FirstIndex := 0;
      while (List[FirstIndex].ParticipantIndex <> VictimIndex) do
         Inc(FirstIndex);
      SecondIndex := FirstIndex + 1;
      while (List[SecondIndex].ParticipantIndex <> VictimIndex) do
         Inc(SecondIndex);
      if (FirstIndex+1 < SecondIndex) then
      begin
         // move middle chunk of list down one
         Move(List[FirstIndex+1], List[FirstIndex], ((SecondIndex-FirstIndex)-1)*SizeOf(TCollisionDetectionDatum));
      end;
      if (SecondIndex < High(List)) then
      begin
         // move last chunk of list down two
         Move(List[SecondIndex+1], List[SecondIndex-1], (High(List)-SecondIndex)*SizeOf(TCollisionDetectionDatum));
      end;
      SetLength(List, Length(List)-2);
      if (Length(List) > 0) then
         for Index := Low(List) to High(List) do
         begin
            if (List[Index].ParticipantIndex > VictimIndex) then
               Dec(List[Index].ParticipantIndex);
         end;
   end;

var
   Participant: PSystemParticipant;
   ParticipantSplicePoint: ^PSystemParticipant;
   ParticipantIndex: Cardinal;
begin
   inherited;
   ParticipantIndex := 0;
   ParticipantSplicePoint := @FParticipants;
   Participant := FParticipants;
   Assert(Assigned(Participant));
   while (Participant^.Actor <> OldChild) do
   begin
      Inc(ParticipantIndex);
      ParticipantSplicePoint := @Participant^.Next;
      Participant := Participant^.Next;
      Assert(Assigned(Participant));
   end;
   ParticipantSplicePoint^ := Participant^.Next;
   Dispose(Participant);
   Dec(FCachedParticipantCount);
   RemoveActorFromFilterList(FCollisionFilterListX, ParticipantIndex);
   RemoveActorFromFilterList(FCollisionFilterListY, ParticipantIndex);
end;

function TGravitationalSystem.AddChild(const NewChild: TAbstractTreeNode; const Location: TAbstractLocation): Boolean;
var
   NewRecord: PSystemParticipant;
begin
   Assert(NewChild is TChildActor);
   Assert(NewChild.Independent);
   New(NewRecord);
   NewRecord^.Actor := NewChild as TChildActor;
   NewRecord^.PosX := XXX; // base this on Location
   NewRecord^.PosY := XXX; // base this on Location
   NewRecord^.VelX := XXX; // base this on Location
   NewRecord^.VelY := XXX; // base this on Location
   InternalAddChild(NewRecord);
   Result := True;
end;

procedure TGravitationalSystem.InternalAddChild(const NewRecord: PSystemParticipant);
{$IFOPT C+}
var
   Result: Boolean;
{$ENDIF}
begin
   Assert(Assigned(NewRecord));
   NewRecord^.Next := FParticipants;
   SetLength(FCollisionFilterListX, Length(FCollisionFilterListX)+2);
   SetLength(FCollisionFilterListY, Length(FCollisionFilterListY)+2);
   // XXX should insert into sorted position, and shift
   FCollisionFilterListX[High(FCollisionFilterListX)-1].Position := NewRecord^.PosX-NewRecord^.Actor.Radius;
   FCollisionFilterListX[High(FCollisionFilterListX)-1].PointType := ptStart;
   FCollisionFilterListX[High(FCollisionFilterListX)-1].ParticipantIndex := FCachedParticipantCount;
   FCollisionFilterListX[High(FCollisionFilterListX)].Position := NewRecord^.PosX+NewRecord^.Actor.Radius;
   FCollisionFilterListX[High(FCollisionFilterListX)].PointType := ptEnd;
   FCollisionFilterListX[High(FCollisionFilterListX)].ParticipantIndex := FCachedParticipantCount;
   FCollisionFilterListY[High(FCollisionFilterListY)-1].Position := NewRecord^.PosY-NewRecord^.Actor.Radius;
   FCollisionFilterListY[High(FCollisionFilterListY)-1].PointType := ptStart;
   FCollisionFilterListY[High(FCollisionFilterListY)-1].ParticipantIndex := FCachedParticipantCount;
   FCollisionFilterListY[High(FCollisionFilterListY)].Position := NewRecord^.PosY+NewRecord^.Actor.Radius;
   FCollisionFilterListY[High(FCollisionFilterListY)].PointType := ptEnd;
   FCollisionFilterListY[High(FCollisionFilterListY)].ParticipantIndex := FCachedParticipantCount;
   FParticipants := NewRecord;
   Inc(FCachedParticipantCount);
   {$IFOPT C+} Result := {$ENDIF} inherited AddChild(NewRecord^.Actor, nil); {$IFOPT C+} Assert(Result); {$ENDIF}
end;

procedure TGravitationalSystem.TeleportToOrbit(Actor: TChildActor; Orbit: TGravitationalSystemValue; Clockwise: Boolean; Around: TChildActor);
var
   NewRecord, AroundRecord: PSystemParticipant;
   Angle, Speed: TGravitationalSystemValue;
begin
   Assert(Assigned(Actor));
   New(NewRecord);
   NewRecord^.Actor := Actor;
   Angle := Random() * 2 * Pi;
   if (Assigned(Around)) then
   begin
      AroundRecord := FParticipants;
      while (Assigned(AroundRecord) and (AroundRecord^.Actor <> Around)) do
        AroundRecord := AroundRecord^.Next;
      Assert(Assigned(AroundRecord));
      Assert(Orbit > Around.Radius);
      Writeln('Orbitting ', Actor.Name, ' around ', Around.Name, '; ', Actor.Mass, 'kg around ', Around.Mass, 'kg');
      //Assert(Around.Mass > Actor.Mass * 1000);
      Speed := SqRt((G*Around.Mass)/Orbit);
      if (Clockwise) then
         Speed := -Speed;
      NewRecord^.PosX := AroundRecord^.PosX + Orbit * Cos(Angle);
      NewRecord^.PosY := AroundRecord^.PosY + Orbit * Sin(Angle);
      NewRecord^.VelX := AroundRecord^.VelX + Speed * Sin(Angle);
      NewRecord^.VelY := AroundRecord^.VelY + Speed * Cos(Angle);
   end
   else
   begin
      NewRecord^.PosX := Orbit * Cos(Angle);
      NewRecord^.PosY := Orbit * Sin(Angle);
      NewRecord^.VelX := 0.0;
      NewRecord^.VelY := 0.0;
   end;
   InternalAddChild(NewRecord);
end;

procedure TGravitationalSystem.Launch(Actor: TChildActor; Around: TChildActor; Angle: TGravitationalSystemValue);
const
   Speed: TGravitationalSystemValue = 0.0;
var
   NewRecord, AroundRecord: PSystemParticipant;
begin
   Assert(Assigned(Actor));
   Assert(Assigned(Around));
   AroundRecord := FParticipants;
   while (Assigned(AroundRecord) and (AroundRecord^.Actor <> Around)) do
     AroundRecord := AroundRecord^.Next;
   Assert(Assigned(AroundRecord));

   New(NewRecord);
   NewRecord^.Actor := Actor;
   NewRecord^.PosX := AroundRecord^.PosX + AroundRecord^.Actor.Radius * Cos(Angle);
   NewRecord^.PosY := AroundRecord^.PosY + AroundRecord^.Actor.Radius * Sin(Angle);
   NewRecord^.VelX := AroundRecord^.VelX + Speed * Cos(Angle);
   NewRecord^.VelY := AroundRecord^.VelY + Speed * Sin(Angle);
   InternalAddChild(NewRecord);
end;

function TGravitationalSystem.GetRendererClass(): UTF8String;
begin
   Result := 'space';
end;

procedure TGravitationalSystem.KillChildren();
var
   Participant: PSystemParticipant;
begin
   inherited;
   Participant := FParticipants;
   while (Assigned(Participant)) do
   begin
      Participant^.Actor.KillChildren();
      Participant := Participant^.Next;
   end;
end;


{$IFDEF DEBUG}
function TGravitationalSystem.GetNameByIndex(Index: Cardinal): UTF8String;
var
   Temp: PSystemParticipant;
begin
   Temp := FParticipants;
   while (Index > 0) do
   begin
      Assert(Assigned(Temp));
      Temp := Temp^.Next;
      Dec(Index);
   end;
   Assert(Assigned(Temp));
   Result := Temp^.Actor.Name;
end;

procedure TGravitationalSystem.DebugPrintFilterLists();
var
   Index: Cardinal;
begin
   System.Write('X: ');
   for Index := Low(FCollisionFilterListX) to High(FCollisionFilterListX) do
   begin
      if (FCollisionFilterListX[Index].PointType = ptStart) then
         System.Write(' <', FCollisionFilterListX[Index].ParticipantIndex, '>')
      else
         System.Write(' </', FCollisionFilterListX[Index].ParticipantIndex, '>')
   end;
   System.Writeln();
   System.Write('Y: ');
   for Index := Low(FCollisionFilterListY) to High(FCollisionFilterListY) do
   begin
      if (FCollisionFilterListY[Index].PointType = ptStart) then
         System.Write(' <', FCollisionFilterListY[Index].ParticipantIndex, '>')
      else
         System.Write(' </', FCollisionFilterListY[Index].ParticipantIndex, '>')
   end;
   System.Writeln();
end;

{$ENDIF}

initialization
{$INCLUDE registrations/gravity.inc}
end.