/***************************************************************
* View the solutions of a differential equations system having *
* the form:   y'=f(x,y,z),  z'=g(x,y,z)                        *
*                                                              *
*               Visual C++ in API style version by J-P Moreau  *
* ------------------------------------------------------------ *
* REFERENCE:                                                   *
* "Graphisme dans le plan et dans l'espace avec Turbo Pascal   *
*  4.0 by R. Dony - MASSON, Paris 1990, pages 254-260 ".       *
* ------------------------------------------------------------ *
* TUTORIAL:                                                    *
* It can be proved that any differential equation of order > 1 *
* may be replaced by a 1st order differential equations system.*
* For example, the 2nd order equation y" = f(x,y,y') is equi-  *
* valent to the system:                                        *
*                         y' = z                               *
*                         z' = f(x,y,z)                        *
*                                                              *
* This program draws the solutions (integral curves) of a sys- *
* tem with the form:                                           *
*                        y' = f(x,y,z)                         *
*                        z' = g(x,y,z)                         *
*                                                              *
* The proposed example is:    y' = x*x*y*y - 1                 *
*                             z' = x*x+y*y - 4                 *
*                                                              *
* The steps to follow are:                                     *
*                                                              *
*     1. choose a physical coordinates window x1,x2,y1,y2      *
*     2. choose an increment step, dx                          *
*     3. for each point of plane (x,y) inside previous window, *
*        calculate by Runge-Kutta method and draw the solution *
*        curve y=f(x) passing through this point.              *
*                                                              *
* Obviously, these curves will eventually go out of screen. It *
* can be managed by the clipping capability of module Gr2D.    *
* However to limit computation time, a curve is stopped when:  *
*                                                              *
*     - the following point will exit graph window,            *
*     - the number of segments already drawn > 300,            *
*     - the curve is approaching a critical point where        *
*     - f(x,y) and g(x,y) ==> 0.                               *
*                                                              *
* Other possible examples:                                     *
*                                                              *
*     1) y' = x + y            z' = x*y                        *
*     2) y' = y(y*y-1)         z' = x(x*x-1)                   *
*     3) y' = -x               z' = x*x + y*y - 1              *
*     4) y' = cos(x)-x*sin(y)  z' = x*sin(y) + y*sin(x)        *
*                                                              *
* ------------------------------------------------------------ *
* NOTE: the variable switch (in Runge-Kutta procedure) is used *
*       to "switch" direction. For each start point in plane   *
*       (x,y), the integral curve can be drawn in two opposite *
*       directions. When switch=1, the curve is drawn in one   *
*       direction until a stop condition is met; at this moment*
*       switch is forced to -1 and the other part of the curve *
*       is drawn until a new stop condition is met. Then switch*
*       is again put to 1 and the process will continue with a *
*       new starting point.                                    *
***************************************************************/
//Use with gr2d.cpp and rk.mak
#include <windows.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

HDC  hdc;
RECT rect;
HPEN hpen, hpenOld;

//Functions used here of module Gr2D.cpp
void Fenetre(double,double,double,double); 
void Cloture(int,int,int,int); 
void Bordure();
void Axes();
void Gradue(double,double);
void Grille(double,double);
void MoveXY(double,double);
void LineXY(double,double);

//"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-5; rct.bottom=iy+15;
  DrawText(hdc, s, -1, &rct,
                DT_SINGLELINE | DT_CENTER | DT_VCENTER);
}


int    nbreseg = 300;
double f1,f2,f3,f4;
double i,j,h,xpos,ypos;


// Example:  y' = x*x*y*y-1,  z' = x*x+y*y-4
// Other functions can be defined by user
double F(double x, double y)  {
  return x*x*y*y-1.0;
}

double G(double x, double y)  {
  return x*x+y*y-4.0;
}

//These values can be changed by user
void Data()  {
  f1=-2.5; f2=2.5; f3=-2.5; f4=2.5;
  h=0.05;
}

// Integration by Runge-Kutta method
void RungeKutta(double y0, double z0)  {
   
  double y,z,y1,z1;
  int    seg,iswitch;
  double pointcrit;
  double k1,k2,k3,k4,l1,l2,l3,l4;

  iswitch=1;
  do {
    y=y0; z=z0;
    MoveXY(y,z);
    y1=y; z1=z; seg=0;
    do {
       LineXY(y1,z1);
       seg++;
       y=y1; z=z1;
       k1=iswitch*h*F(y,z);
       l1=iswitch*h*G(y,z);
       y=y1+k1/2.0; z=z1+l1/2.0;
       k2=iswitch*h*F(y,z);
       l2=iswitch*h*G(y,z);
       y=y1+k2/2.0; z=z1+l2/2.0;
       k3=iswitch*h*F(y,z);
       l3=iswitch*h*G(y,z);
       y=y1+k3/2.0; z=z1+l3/2.0;
       k4=iswitch*h*F(y,z);
       l4=iswitch*h*G(y,z);
       y1=y+iswitch*h/6.0*(k1+2*k2+2*k3+k4);
       z1=z+iswitch*h/6.0*(l1+2*l2+2*l3+l4);
       pointcrit=sqrt(F(y,z)*F(y,z)+G(y,z)*G(y,z));
       y1=y+iswitch*h*F(y,z);
       z1=z+iswitch*h*G(y,z);
    } while (y1>=f1 && y1<=f2 && z1>=f3 && z1<=f4 && seg<nbreseg && pointcrit>1e-4);
    iswitch=-iswitch;
  } while (iswitch!=1);
} //RungeKutta()


void Initgraph() {
  extern MaxX,MaxY;
  HPEN redpen;
  redpen = CreatePen(PS_SOLID, 1, RGB(255, 0, 0));
  Fenetre(f1,f2,f3,f4);
  Cloture(60,MaxX-10,40,MaxY-10);
  Gradue((f2-f1)/10,(f4-f3)/10);
  SelectObject(hdc,redpen);
  Grille((f2-f1)/10,(f4-f3)/10);
  SelectObject(hdc,hpen);
  Axes();
  Bordure();
}

//visualize Runge-Kutta solutions
void Draw_RungeKutta() {

  extern MaxX, MaxY;

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

  Data();
  Initgraph();
  i=0.0;
  while (i < 1) {
    i+=h;
    xpos=f1+i*(f2-f1);
    j=0;
    while (j < 1) {
      j+=h;
      ypos=f3+j*(f4-f3);
      RungeKutta(xpos,ypos);
    }
  }

} //Draw_RungeKutta()


//Handle window Paint and Destroy messages
long FAR PASCAL /*_export*/ WndProc(HWND hwnd,
    UINT message, UINT wParam, LONG lParam) 
{
   
    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_RungeKutta();
               
            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[] = "RK";
    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
        " RUNGE KUTTA",             // 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 rk.cpp