/*******************************************************************
*             Differential equations of order N                    *
*             by Runge-Kutta method of order 4                     *
* ---------------------------------------------------------------- *
* Reference: "Analyse en Turbo Pascal versions 5.5 et 6.0" By Marc *
*             DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991.    *
*             [BIBLI 03].                                          *
*                                                                  *
*                             C++ Version By J-P Moreau, Paris     *
*                         (use with teqdifn1.mak and graph2d.cpp)  *
*                                    (www.jpmoreau.fr)             *
* ---------------------------------------------------------------- *
* SAMPLE DATA:                                                     *
*                                                                  *
* Example: integrate y"=(2y'y'+yy)/y from x=4 to x=6               *
*          with initial conditions: y(4)=2 and y'(4)=-2tan(1)      *
*                                                                  *
* Exact solution is:   y = 2cos(1)/cos(x-5)                        *
*                                                                  *
*        DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER N          *
*              of type y(n) = f(y(n-1),...,y"",y,x)                *
*                                                                  *
*   order of equation: 2                                           *
*   begin value x    : 4                                           *
*   end value x      : 6                                           *
*   y0 value at x0   : 2                                           *
*   y1 value at x0   : -3.114815                                   *
*   number of points : 11                                          *
*   finesse          : 20                                          *
*                                                                  *
*       X            Y                                             *
* ---------------------------                                      *
*   4.000000     2.000000                                          *
*   4.200000     1.551018                                          *
*   4.400000     1.309291                                          *
*   4.600000     1.173217                                          *
*   4.800000     1.102583                                          *
*   5.000000     1.080605                                          *
*   5.200000     1.102583                                          *
*   5.400000     1.173217                                          *
*   5.600000     1.309291                                          *
*   5.800000     1.551018                                          *
*   6.000000     2.000000                                          *
*                                                                  *
*******************************************************************/
#include <windows.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>


HDC  hdc;
RECT rect;
HPEN hpen, hpenOld;

//Functions used here of Graph2D.cpp
void CourbeXY (int,int,float *,double,double);
void Legendes (int,char *,char *,char *,int,int,char *,int,int,char *);


//"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;
  double dx,dy;
  if (ix2<ix1) {
	Swap(&ix1,&ix2); Swap(&iy1,&iy2);
  }
  il=(int) sqrt((ix2-ix1)*(ix2-ix1)+(iy2-iy1)*(iy2-iy1));
  dx=(double) (ix2-ix1)/il; 
  dy=(double) (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    fi,i,ndata,ordre;
double xi,xf;
double yi[10];
char   s[2];

// Example: y"=(2y'y'+yy)/y
double fp(double x, double *y) {
  if (fabs(y[0])<1e-12) y[0]=1e-12;
  return (2*y[1]*y[1]+y[0]*y[0])/y[0];
}

/**************************************************************************
*        SOLVING DIFFERENTIAL EQUATIONS WITH 1 VARIABLE OF ORDER N        *
*                of type y(n) = f(y(n-1),y(n-2),...,y',y,x)               *
* ----------------------------------------------------------------------- *
*  INPUTS:                                                                *
*    f         Equation to integrate                                      *
*    xi, xf    Begin, end values of variable x                            *
*    Yi        Begin values at xi (of f and derivatives)                  *
*    m         Number of points to calculate                              *
*    n         Order of differential equation                             *
*    fi        finesse (number of intermediate points)                    *
*                                                                         *
*  OUTPUTS:                                                               *
*    t         real vector storing m results for function y               *
* ----------------------------------------------------------------------- *      
*  EXAMPLE:    y" = (2 y'y' + yy) / y with y(4) = 2, y'(4) = -2*tan(1)    *
*              Exact solution:  y = (2 cos(1))/cos(x-5)                   *
**************************************************************************/
void equadiffn(double xi,double xf,double *yi,int m,int n,int fi)  {
  double h,x,a,b,c,d;
  double ta[10],tb[10],tc[10],y[10],z[10];
  float  t[101];
  int    i,j,k,ni,n1,n2;
  extern MaxX;
   if (fi<1) return;
   h = (xf - xi) / fi / (m-1);
   n1=n-1; n2=n-2;
   t[1]=(float) yi[0];
   for (k=0; k<n; k++) {
     y[k]=yi[k]; z[k]=yi[k];
   }
   for (i=1; i<m+1; i++) {
     ni=(i-1)*fi-1;
     for (j=1; j<fi+1; j++) {
       x=xi+h*(ni+j);
       for (k=0; k<n; k++) y[k]=z[k];
       a=h*fp(x,y);
       for (k=0; k<n1; k++) {
         ta[k]=h*y[k+1]; y[k]=z[k]+ta[k]/2;
       }
       y[n1]=z[n1]+a/2;
       x=x+h/2;
       b=h*fp(x,y);
       for (k=0; k<n1; k++) {
         tb[k]=h*y[k+1]; y[k]=z[k]+tb[k]/2;
       }
       y[n1]=z[n1]+b/2;
       c=h*fp(x,y);
       for (k=0; k<n1; k++) {
         tc[k]=h*y[k+1]; y[k]=z[k]+tc[k];
       }
       y[n1]=z[n1]+c;
       x=x+h/2;
       d=h*fp(x,y);
       for (k=0; k<n1; k++) 
         z[k]=z[k]+(ta[k]+2*tb[k]+2*tc[k]+h*y[k+1])/6;
       z[n1]=z[n1]+(a+b+b+c+c+d)/6;
	 }
     t[i+1]=(float) z[0];
   }

   //Draw result y(x) stored in table t
   for (i=0; i<m; i++)  t[i]=t[i+1]; //shift index of t
   CourbeXY(m,10,t,xi,xf);
   Legendes(10," SOLUTION Y = F(X) ","X","Y",MaxX-350,50,
	           "Equation: y'' = (2 y'y' + yy) / y",MaxX-275,70,"from x=4 to x=6");
}

void Run_equadiffn()  {
  extern MaxX, MaxY;

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

  //data section
  ordre=2;            //order of equation 
  ndata=100;          //number of tabulated points of solution
  fi=20;              //number of intermediate calculated points
  xi=4; xf=6;         //x begin, x end
  yi[0]=2;            //initial value of y = y(xi)
  yi[1]=-3.114815;    //initial value of y'= y'(xi) 

  equadiffn(xi,xf,yi,ndata,ordre,fi);
}


//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);

            Run_equadiffn();
               
            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[] = "Teqdifn";
    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
        " DIFFERENTIAL EQUATIONS WITH 1 VAR OF ORDER N", // 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 teqdifn1.cpp