/* ***********************************************************************
   *   									 *
   *	         	Three Particle Orbit Simulator			 *
   *									 *
   *									 *
   *				Program by				 *
   *			   Christopher D. Watkins			 *
   *									 *
   *		              'C' conversion by				 *
   *				 Larry Sharp				 *
   *									 *
   ***********************************************************************

   increase time step "dt" to greater than zero to increase speed
*/

#include "stdio.h"
#include "dos.h"
#include "conio.h"
#include "math.h"
#include "math.inc"
#include "graph2.inc"

int MSize;
int x,y;
int oVx1, oVy1;
int oAx1, oAy1;
int oVx2, oVy2;
int oAx2, oAy2;
int oVx3, oVy3;
int oAx3, oAy3;

void InitMeters()
{
  MSize=30;
  Circle(MSize, MSize, MSize-4, 35);
  oVx1=MSize;
  oVy1=MSize;
  oAx1=MSize;
  oAy1=MSize;
  oVx2=MSize;
  oVy2=MSize;
  oAx2=MSize;
  oAy2=MSize;
  oVx3=MSize;
  oVy3=MSize;
  oAx3=MSize;
  oAy3=MSize;
}

void Meter(float dx, float dy, int *xt, int *yt, int Radius, int Color)
{
  float r;

  Draw(MSize, MSize, *xt, *yt, 0);
  r=Radius/sqrt(SqrFP(dx)+SqrFP(dy))-10;
  x=Round(dx*r)+MSize;
  y=-Round(dy*r)+MSize;
  *xt=x;
  *yt=y;
  Draw(MSize, MSize, x, y, 15);
}

float X1, Y1, Vx1, Vy1, Ax1, Ay1;
float X2, Y2, Vx2, Vy2, Ax2, Ay2;
float X3, Y3, Vx3, Vy3, Ax3, Ay3;
float D12, D23, D31, dt;
float Dx12, Dx23, Dx31;
float Dy12, Dy23, Dy31;
float Tx12, Tx23, Tx31;
float Ty12, Ty23, Ty31;
int   Xp1, Yp1, Xp2, Yp2, Xp3, Yp3;
int   M1, M2, M3;

void main()
{
  M1=1;
  M2=10;
  M3=3;

  X1=-40.0;	Y1=0.0;
  X2= 00.0;	Y2=0.0;
  X3= 90.0;	Y3=0.0;

  Vx1= 0.1010;	Vy1= 0.2500;
  Vx2= 0.0010;	Vy2= 0.0010;
  Vx3=-0.0200;	Vy3=-0.1010;

  dt=0.4;

  InitGraphics(18);
  InitMeters();
  do
  {
    Xp1=160+Trunc(X1);
    Yp1=100-Trunc(Y1);
    Xp2=160+Trunc(X2);
    Yp2=100-Trunc(Y2);
    Xp3=160+Trunc(X3);
    Yp3=100-Trunc(Y3);

    Draw(Xp1, Yp1, Xp2, Yp2, 35);
    Draw(Xp2, Yp2, Xp3, Yp3, 35);
    Draw(Xp3, Yp3, Xp1, Yp1, 35);

    Circle(Xp1, Yp1, M1, 143);
    Circle(Xp2, Yp2, M2, 71);
    Circle(Xp3, Yp3, M3, 251);

    X1+=Vx1*dt;
    Y1+=Vy1*dt;
    X2+=Vx2*dt;
    Y2+=Vy2*dt;
    X3+=Vx3*dt;
    Y3+=Vy3*dt;

    Dx12=X1-X2;
    Dy12=Y1-Y2;
    Dx23=X2-X3;
    Dy23=Y2-Y3;
    Dx31=X3-X1;
    Dy31=Y3-Y1;

    D12=sqrt(SqrFP(Dx12)+SqrFP(Dy12));
    D12=1.0/(D12*D12*D12);
    D23=sqrt(SqrFP(Dx23)+SqrFP(Dy23));
    D23=1.0/(D23*D23*D23);
    D31=sqrt(SqrFP(Dx31)+SqrFP(Dy31));
    D31=1.0/(D31*D31*D31);

    Tx31=Dx31*D31;
    Ty31=Dy31*D31;
    Tx12=Dx12*D12;
    Ty12=Dy12*D12;
    Tx23=Dx23*D23;
    Ty23=Dy23*D23;

    Ax1=(M3*Tx31-M2*Tx12);
    Ay1=(M3*Ty31-M2*Ty12);
    Ax2=(M1*Tx12-M3*Tx23);
    Ay2=(M1*Ty12-M3*Ty23);
    Ax3=(M2*Tx23-M1*Tx31);
    Ay3=(M2*Ty23-M1*Ty31);

    Vx1+=Ax1*dt;
    Vy1+=Ay1*dt;
    Vx2+=Ax2*dt;
    Vy2+=Ay2*dt;
    Vx3+=Ax3*dt;
    Vy3+=Ay3*dt;

    Meter(Vx1, Vy1, &oVx1, &oVy1, 25, 123);
    Meter(Ax1, Ay1, &oAx1, &oAy1, 12, 143);
    Meter(Vx2, Vy2, &oVx2, &oVy2, 25, 51);
    Meter(Ax2, Ay2, &oAx2, &oAy2, 12, 71);
    Meter(Vx3, Vy3, &oVx3, &oVy3, 25, 231);
    Meter(Ax3, Ay3, &oAx3, &oAy3, 12, 251);

    Draw(Xp1, Yp1, Xp2, Yp2, 0);
    Draw(Xp2, Yp2, Xp3, Yp3, 0);
    Draw(Xp3, Yp3, Xp1, Yp1, 0);
    delay(5);

    Circle(Xp1, Yp1, M1, 0);
    Circle(Xp2, Yp2, M2, 0);
    Circle(Xp3, Yp3, M3, 0);
  }
  while(!(kbhit()));
  ExitGraphics();
}