/************************************************************************
*     PLANETS.CPP        Version 1.0 (Visual C++) 09/05/2000            *
* --------------------------------------------------------------------- *
* This program calculates the motions of N bodies submitted to their    *
* gravity attractions. The integration of motions is made by a Runge-   *
* Kutta method with controled time step to satisfy a required precision.*
*                                                                       *
*                                     C++ version by J-P Moreau         *
*                         (from a DOS program in pascal by J-P Dumont)  *
*                                         (www.jpmoreau.fr)             *
*                                                                       *
* Uses: Planets.mak and Gr2d.cpp                                        *
************************************************************************/
#include <windows.h>
#include <stdio.h>
#include <math.h>

HDC  hdc;
RECT rect;

//Functions of module Gr2D.cpp
void Fenetre(double,double,double,double); 
void Cloture(int,int,int,int); 
void Bordure();
void Cercle(double,double,double,int);


//"home made" graphic commands for hdc environment used by above functions
void DrawPixel(int ix,int iy) {
  //sorry, no other available command found
  Rectangle(hdc,rect.left+ix,rect.top+iy,
	            rect.left+ix+2,rect.top+iy+1);
}

void Swap(int *i1,int *i2) {
  int it;
  it=*i1; *i1=*i2; *i2=it;
}

void DrawLine(int ix1,int iy1,int ix2,int iy2) {
  int i,il,ix,iy;
  float dx,dy;
  if (ix2<ix1) {
	Swap(&ix1,&ix2); Swap(&iy1,&iy2);
  }
  il=(int) sqrt((ix2-ix1)*(ix2-ix1)+(iy2-iy1)*(iy2-iy1));
  dx=(float) (ix2-ix1)/il; 
  dy=(float) (iy2-iy1)/il;
  ix=ix1; iy=iy1;
  for (i=1;i<il+1; i++) {
	DrawPixel(ix,iy);
	ix=(int) (ix1+i*dx); 
	iy=(int) (iy1+i*dy);
  }
}

void OutText(int ix,int iy,char *s) {
  RECT rct;
  rct.left=ix; rct.right=ix+10*strlen(s);
  rct.top=iy; rct.bottom=iy+10;
  DrawText(hdc, s, -1, &rct,
                DT_SINGLELINE | DT_CENTER | DT_VCENTER);
}

// ==== Begin planets section ====

int NcorpsMx = 8,          // Maximal number of bodies
    Nddl     = 4,          // Number of degrees of freedom per body
    NsmbMx   = 32;         // Maximal number of second members

double stab  = 0.05;

double ystart[33];         // Instantaneous state vector (index 0 not used)

double dt,dtmin,dtsave,eps,h1,Tmax,time,t2,XRayon,YRayon;

long   istep,nstep,NstepMax;
int    Done,Ncorps,Nsmb,nok,nbad;


////////////////////////////////////////////////////////////////////////
void initpointers (int *ipt)  {
//======================================================================
//Calculate pointers to store state variables in Table
//----------------------------------------------------------------------
  int i;

  ipt[0] = - Nddl;
  for (i=1; i<=Ncorps; i++)
    ipt[i] = ipt[i-1] + Nddl;
}

//----------------------------------------------------------------------
void ReadData ( double *Masse, double *q1, double *q2,
			    double *q1d,   double *q2d,
                double *Tmax,  double *dtsave,
                int    *NCorps )  {
    int i;
    double sign;

    *Tmax     = 100;
    *dtsave   = 0.002;
    *NCorps   = 8;        // number of planets (max 8)
    sign      = 1.0;

    Masse[1] = 10000;
    q1[1]    = 0;
    q2[1]    = 0;
    q1d[1]   = 0;
    q2d[1]   = 0;

    for (i=2; i<=*NCorps; i++) {
      sign     = -sign;
      Masse[i] = 1.0*(9.0 - i);
      q1[i]    = 1.0*(i-1.0);
      q2[i]    = 0;
      // Factor 1.0 circular trajectory
      // Factor 1.414 parabolic trajectory
      q2d[i]   = sign*1.38*sqrt(Masse[1]/q1[i]);
      q1d[i]   = 0;
      // centrifugal force = gravity force
      q1[i]    = q1[1]+sign*q1[i];
    }
} // ReadData()

//-----------------------------------------------------------------------
void InitialConditions(double *XN,  double *q1, double *q2,
					   double *q1d, double *q2d)  {

  int i,ipti;
  extern int ipt[32];

  for (i=1; i<= Ncorps; i++) {
    ipti = ipt[i];
    XN[ipti+1] = q1d[i];
    XN[ipti+2] = q2d[i];
    XN[ipti+3] = q1[i];
    XN[ipti+4] = q2[i];
  }
} // InitialConditions()

//------------------------------------------------------------------------
void Display(double *XN)  {

  double xw,yw;
  int i,ipti;

  extern int ipt[32];

  for (i=1; i<=Ncorps; i++) {
    ipti = ipt[i];
    xw = XN[ipti+3];  //xi
    yw = XN[ipti+4];  //yi

    // Draw central star at first time
    if (i==1) {
      if (Done==0) {
        Cercle(0,0,XRayon,0);
        Done=1;
      }
    }
    //draw planets around central star
    else
      Cercle(xw,yw,YRayon,0);
  } // i loop
}

// Function of unit equa_dif.cpp
void odeint1(double *, int,
             double, double, double, double, double,
             int *, int *);

void Draw_planets()  {

// LABEL: fin;

    long   tt;

    extern double Masse[33],q1[33],q1d[33],q2[33],q2d[33];
	extern int ipt[32];

    extern MaxX, MaxY;
    //client drawing zone in pixels
    MaxX=(rect.right-rect.left);
    MaxY=(rect.bottom-rect.top);

    Done=0;

    //define graph window in physical coordinates
    Fenetre(-50,50,-40,40);
    //define graph window in pixels
    Cloture(10,MaxX-10,10,MaxY-10);
    //draw a rectangle around graph zone 
    Bordure();

    ReadData(Masse,q1,q2,q1d,q2d,&Tmax,&dtsave,&Ncorps);

    Nsmb = Nddl*Ncorps;

    NstepMax = (long) floor(Tmax/dtsave);

    initpointers (ipt);

    // Define useful constants
    XRayon = 0.5;                               //Radius of central star
    YRayon = 0.25;                                   //Radius of planets

    eps   = 1e-14;                          //Required relative precision
    h1    = 0.01*dtsave;
    dt    = h1;
    dtmin = 0.001*h1;                       //Smallest time step allowed
    time  = 0;
    istep = 0;

    InitialConditions(ystart,q1,q2,q1d,q2d);

    // main integration loop
    while ((time <= Tmax) && (istep < NstepMax)) {

      istep++;
      t2=time+dtsave;

      odeint1(ystart,Nsmb,time,t2,eps,dt,dtmin,&nok,&nbad);

      if (nok < 0) {
        istep--;
       	goto fin;          // no storage
      }
      time=t2;
      // display planets once every four times
	  tt = (long) floor(t2*100);
      if ((tt % 4)==0)  Display(ystart);
    }

    fin:  nstep=istep;     // actual number of time steps
    return;
}

// ==== End planets section ====


//Handle window Paint and Destroy messages
long FAR PASCAL /*_export*/ WndProc(HWND hwnd,
    UINT message, UINT wParam, LONG lParam) 
{
    HPEN        hpen, hpenOld;
    PAINTSTRUCT ps;

    switch (message) {
        case WM_PAINT:
	    hdc = BeginPaint(hwnd, &ps);
	    GetClientRect(hwnd, &rect);
            hpen = CreatePen(PS_SOLID, 1, RGB(0, 0, 255));
            hpenOld = SelectObject(hdc,hpen);

            Draw_planets();
               
            SelectObject(hdc,hpenOld);
            DeleteObject(hpen);
            EndPaint(hwnd, &ps);
            return 0;
        case WM_DESTROY:
            PostQuitMessage(0);
            return 0;
    }
    return DefWindowProc(hwnd, message, wParam, lParam);
}

// main program
int PASCAL WinMain(HANDLE hInstance, HANDLE hPrevInstance,
    LPSTR lpszCmdParam, int nCmdShow)
{
    static char szAppName[] = "Planets";
    HWND        hwnd;
    MSG         msg;
    WNDCLASS    wndclass;

    if (!hPrevInstance) {
        wndclass.style         = CS_HREDRAW | CS_VREDRAW;
        wndclass.lpfnWndProc   = WndProc;
        wndclass.cbClsExtra    = 0;
        wndclass.cbWndExtra    = 0;
        wndclass.hInstance     = hInstance;
        wndclass.hIcon         = LoadIcon(NULL, IDI_APPLICATION);
        wndclass.hCursor       = LoadCursor(NULL, IDC_ARROW);
        wndclass.hbrBackground = GetStockObject(WHITE_BRUSH);
        wndclass.lpszMenuName  = NULL;
        wndclass.lpszClassName = szAppName;

        RegisterClass(&wndclass);
    }

    hwnd = CreateWindow(szAppName,  // window class name
        " SEVEN  PLANETS",          // window caption
        WS_OVERLAPPEDWINDOW,        // window style                       
        CW_USEDEFAULT,              // initial x position
        CW_USEDEFAULT,              // initial y position
        CW_USEDEFAULT,              // initial x size
        CW_USEDEFAULT,              // initial y size
        NULL,                       // parent window handle
        NULL,                       // window menu handle
        hInstance,                  // program instance handle
        NULL);                      // creation parameters

   ShowWindow(hwnd, nCmdShow);
   UpdateWindow(hwnd);

   while (GetMessage(&msg, NULL, 0, 0)) {
       TranslateMessage(&msg);
       DispatchMessage(&msg);
   }
   return msg.wParam;
}

//end of file planets.cpp