unit uNBody;
{
  COSBI: Comprehensive Open Source Benchmarking Initiative
  Copyright (c) 2004 Van Smith

  This program is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License
  as published by the Free Software Foundation; either version 2
  of the License, or (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

  The current web address of the GNU General Public License is:
  http://www.gnu.org/licenses/gpl.html

  You can contact the authors of this software at:
  cosbi@vanshardware.com
  See www.vanshardware.com or www.cosbi.org for more contact details.
}
//==============================================================================
// Unit name: uNBody
// Unit description: OpenGL version of the nBody test.
// Author: Van Smith
// Date: March 16, 2005
// OS dependent: Yes: Windows
// Resolution dependent: No, but resolution and color depth will impact scores.
// External unit dependencies: GLScene, GLObjects,
//   GLWin32Viewer, GLMisc, OpenGL1x, GLTexture, VectorGeometry, GLGraph,
//   GLGeomObjects, ExtCtrls, StdCtrls, GLFireFX, GLCadencer;
//==============================================================================
// Modification History
// ------------ -------
// Ver  Date   Who    Description
// ==== ====== ====== ==========================================================
// 1.0 050316 Van     Created.
//==============================================================================

interface

uses
  Windows, Messages, SysUtils, Classes, Controls, Forms, GLScene, GLObjects,
  GLWin32Viewer, GLMisc, OpenGL1x, GLTexture, VectorGeometry, GLGraph,
  GLGeomObjects, ExtCtrls, StdCtrls, GLFireFX, GLCadencer;

type
  TfrmNBody = class(TForm)
    GLScene1: TGLScene;
    SceneViewer: TGLSceneViewer;
    GLCamera: TGLCamera;
    GLPlane: TGLPlane;
    GLPoints: TGLPoints;
    DirectOpenGL: TGLDirectOpenGL;
    GLXYZGrid1: TGLXYZGrid;
    Timer1: TTimer;
    GLSphere1: TGLSphere;
    GLSphere2: TGLSphere;
    GLSphere3: TGLSphere;
    GLSphere4: TGLSphere;
    GLLightSource2: TGLLightSource;
    GLLightSource1: TGLLightSource;
    GLFireFXManager1: TGLFireFXManager;
    GLCadencer1: TGLCadencer;
    procedure FormCreate(Sender: TObject);
//    procedure DirectOpenGLRender(Sender : TObject; var rci: TRenderContextInfo);
    procedure Timer1Timer(Sender: TObject);
    procedure FormKeyPress(Sender: TObject; var Key: Char);
  private
    { Private declarations }
    fIndex : integer;
    fAbort : Boolean;
    procedure BallGo;
    procedure RunTest;
    procedure IntializeFPU;
  public
    { Public declarations }
    mx, my : Integer;
  end;
  // required for TnBodyTest below
  TBody = record
    m:   extended;
    x:   extended;
    y:   extended;
    Vx:  extended;
    Vy:  extended;
    Fx:  extended;
    Fy:  extended;
  end; // TBody

var
  frmNBody: TfrmNBody;

implementation

{$R *.dfm}

const

  L_STEPS = 10000;

procedure TfrmNBody.IntializeFPU;
asm
  finit
  wait
end;


procedure TfrmNBody.FormCreate(Sender: TObject);
var
   j : Integer;
   x0, y0, z0, x1, y1, z1 : double;
   h, a, b, c : double;

begin
  fIndex := 0;
  fAbort := FALSE;
  h := 0.01;
  a := 10;
  b := 28;
  c := 8 / 3;
  x0 := 0.1;
  y0 := 0;
  z0 := 0;

  IntializeFPU;

  for j := 0 to L_STEPS do begin
    x1 := x0 + h * a * (y0 - x0);
    y1 := y0 + h * (x0 * (b - z0) - y0);
    z1 := z0 + h * (x0 * y0 - c * z0);
    x0 := x1;
    y0 := y1;
    z0 := z1;
//      glColor3f( 1, 1, 0);
    GLPoints.Positions.Add(x0, y0, z0);
  end; // for j
//    generate a bunch of random points
  for j := 1 to 400 do begin
    GLPoints.Positions.Add( (Random-0.5) * 12, (Random-0.5) * 10 - 2, 4 + (Random-0.5)* 20);
  end; // for j
end;

//procedure TfrmNBody.DirectOpenGLRender(Sender : TObject; var rci: TRenderContextInfo);
//var
//   i : Integer;
//   mat : TMatrix;
//   p, pProj : TVector;
//   planePoint, planeNormal : TVector;
//   plane : THmgPlane;
////   planePoint1, planeNormal1 : TVector;
////   plane1 : THmgPlane;
//begin
////   // Here we recover our plane point and normal...
////   planePoint:=GLPlane.Position.AsVector;
////   planeNormal:=GLPlane.Direction.AsVector;
////   // ...which we use to create a plane (equation)
////   plane:=PlaneMake(planePoint, planeNormal);
////   // from that plane equation and our pojection direction
////   // (which is here the plane normal)
////   mat:=MakeParallelProjectionMatrix(plane, planeNormal);
//{
//   // Here we recover our plane point and normal...
//   planePoint1:=GLPlane1.Position.AsVector;
//   planeNormal1:=GLPlane1.Direction.AsVector;
//   // ...which we use to create a plane (equation)
//   plane1:=PlaneMake(planePoint1, planeNormal1);
//   // from that plane equation and our pojection direction
//   // (which is here the plane normal)
//   mat:=MakeParallelProjectionMatrix(plane1, planeNormal1);
//}
//   // save state, turn off lighting and specify the lines color
////   glPushAttrib(GL_ENABLE_BIT);
////   glDisable(GL_LIGHTING);
//
//   // we'll be drawing a bunch of lines, to specify a line in OpenGL,
//   // you only need to specify the line start and end vertices
////   glBegin(GL_LINES);
//{      for i:=0 to GLPoints.Positions.Count-2 do begin
//         glColor3f( i / L_STEPS , i /L_STEPS, 1);
//         // read the point coordinates, directly from the TGLPoints list
//         MakePoint(p, GLPoints.Positions.List[i]);
//         // project this point on the plane with the matrix
//         MakePoint(pProj, GLPoints.Positions.List[i + 1]);
////         pProj:=VectorTransform(p, mat);
//         // specify the two vertices
//         glVertex3fv(@p);
//         glVertex3fv(@pProj);
//      end;
//}
////   glEnd;
//
//   // restore state
////   glPopAttrib;
//end;

procedure TfrmNBody.Timer1Timer(Sender: TObject);
begin
  timer1.Enabled := FALSE;
  RunTest;
end;

procedure TfrmNBody.BallGo;
var
  i : integer;
begin
  for i := 0 to L_STEPS do begin;
    GLCamera.MoveAroundTarget(0.00, 0.04);
    inc( fIndex );
    if fIndex < L_STEPS then begin
      GLSphere1.Position.AsAffineVector := GLPoints.Positions.List[ i ];
    end else begin
      fIndex := 0;
    end;
    Application.ProcessMessages;
  end; // for
  Application.Terminate;
end; // procedure TForm1.BallGo;


procedure TfrmNBody.RunTest;
const
  G             = 0.01;
  NUM_BODIES    = 3;
  END_TIME      = 30;
  BLEND_FREQ    = 100000;
var
  i     : integer;
  j     : integer;
  timesteps : integer;
  n     : integer;
  p     : array[0..NUM_BODIES] of TBody;
  dt    : extended;
  t     : extended;
  theta : extended;
  F     : extended;
  x,y,r : extended;
  dVx, dVy : extended;

  procedure MoveBodies;
  begin
    GLSphere1.Position.X := p[0].x;
    GLSphere1.Position.z := p[0].Y;
    GLCamera.Position.X := p[0].x - 2;
    GLCamera.Position.z := p[0].y - 2;
    GLCamera.Position.y := 2;
    GLCamera.Direction.X := p[0].x;
    GLCamera.Direction.z := p[0].y;
    GLCamera.Direction.y := 0;
    GLLightSource2.Position.X := p[0].x;
    GLLightSource2.Position.z := p[0].y;
    GLLightSource2.Position.y := 0;
    GLSphere2.Position.X := p[1].x;
    GLSphere2.Position.z := p[1].Y;
    GLSphere3.Position.x := p[2].x;
    GLSphere3.Position.z := p[2].Y;
    GLSphere4.Position.X := p[3].x;
    GLSphere4.Position.z := p[3].Y;
    GLCadencer1.Progress;
    if fAbort then begin
      Application.Terminate;
      Abort;
    end;
    Application.ProcessMessages;
  end; // procedure MoveBodies;

begin
  // intitial conditions
  //mass 0:
  p[0].m  := 2000;
  p[0].x  := 0;
  p[0].y  := 0;
  p[0].Vx := 0;
  p[0].Vy := 0.15; //-0.084;
  //mass 1:
  p[1].m  := 0.02 * p[0].m;
  p[1].x  := 2;
  p[1].y  := 0;
  p[1].Vx := 0;
  p[1].Vy := 0.9 * sqrt( G * p[0].m / 2 );
  //mass 2:
  p[2].m  := 0.002 * p[0].m; //0.00001 * p[0].m;//0.002 * p[0].m; //
  p[2].x  := 2.1;
  p[2].y  := 0;
  p[2].Vx := 0;
  p[2].Vy := p[1].Vy + sqrt( G * p[1].m / 0.1);
  //mass 3:
  p[3].m  := 0.06 * p[0].m;
  p[3].x  := -2.5;
  p[3].y  := 0;
  p[3].Vx := 0;
  p[3].Vy := sqrt( G * p[0].m / 2.5 );

  MoveBodies;

  dt := 0.00001;
  t := 0;
  timesteps := round( END_TIME / dt );
  for j := 1 to timesteps do begin
    t := t + dt;
    // zero out forces for all bodies:
    for n := 0 TO NUM_BODIES do begin
      p[n].Fx := 0;
      p[n].Fy := 0;
    end;
    for n := 0 TO NUM_BODIES do begin
      for i := n to NUM_BODIES do begin
        if n <> i then begin
          r := sqrt( sqr( p[i].x - p[n].x ) + sqr( p[i].y - p[n].y ) );
          F := G * p[i].m * p[n].m / (r * r);
          y := p[i].y - p[n].y;
          x := p[i].x - p[n].x;
          theta := ArcTan(y / x);
          if x < 0 then theta := theta + PI;
          // the force on the body by another
          p[n].Fx := p[n].Fx + F * cos(theta);
          p[n].Fy := p[n].Fy + F * sin(theta);
          // the force that the body exerts on another
          p[i].Fx := p[i].Fx - F * cos(theta);
          p[i].Fy := p[i].Fy - F * sin(theta);
        end; // if n <> i
      end; // next i;
      dVx := p[n].Fx * dt / p[n].m;
      dVy := p[n].Fy * dt / p[n].m;
      p[n].Vx := p[n].Vx + dVx;
      p[n].Vy := p[n].Vy + dVy;
      p[n].x := p[n].x + p[n].Vx * dt;
      p[n].y := p[n].y + p[n].Vy * dt;
      if ( j mod 800 = 0 ) and ( n = NUM_BODIES ) then begin
        MoveBodies;
      end; // if j mod 100 = 0 then
    end; // NEXT n
  end; // NEXT j
  Application.Terminate;
end; // procedure TfrmNBody.RunTest;

procedure TfrmNBody.FormKeyPress(Sender: TObject; var Key: Char);
begin
  if Key = #27 then begin
    fAbort := TRUE;
  end;
end;

end.
