/*******************************************************************
*        Differential equations of order 1 with 1 variable         *
*                by Runge-Kutta method of order 4                  *
*                                                                  *
*                   C++ version with graphic result by J-P Moreau  *
*                      (use with teqdif1a.mak and graph2d.cpp)     *
*                                (www.jpmoreau.fr)                 *
* ---------------------------------------------------------------- *
* Reference: "Analyse en Turbo Pascal versions 5.5 et 6.0 By Marc  *
*             DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991"    *
*             [BIBLI 03].                                          *
* ---------------------------------------------------------------- *
* SAMPLE DATA:                                                     *
*                                                                  *
* Example: integrate y'=4x*(y+sqrt(y))/1+x^2 from x=0 to x=1       *
*                                                                  *
*        DIFFERENTIAL EQUATION WITH 1 VARIABLE OF ORDER 1          *
*                     of type y" = f(x,y)                          *
*                                                                  *
*   begin value x   : 0                                            *
*   end value x     : 1                                            *
*   y value at x0   : 1                                            *
*   number of points: 11                                           *
*   finesse         : 10                                           *
*                                                                  *
*       X            Y                                             *
* -------------------------                                        *
*   0.000000     1.000000                                          *
*   0.100000     1.040400                                          *
*   0.200000     1.166400                                          *
*   0.300000     1.392400                                          *
*   0.400000     1.742400                                          *
*   0.500000     2.250000                                          *
*   0.600000     2.958400                                          *
*   0.700000     3.920400                                          *
*   0.800000     5.198400                                          *
*   0.900000     6.864400                                          *
*   1.000000     9.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);
}

#define SIZE  101

int fi,m;
double xi,xf,yi;
float  T[SIZE];


double fp(double x, double y) {
  //Example: y'=4x(y+root(y))/(1+x²)
  return 4*x*(y+sqrt(y))/(1+x*x);
}


/**************************************************************************
*        SOLVING DIFFERENTIAL EQUATIONS WITH 1 VARIABLE OF ORDER 1        *
*                       of type y' = f(x,y)                               *
* ----------------------------------------------------------------------- *
*  INPUTS:                                                                *
*    fp        Given equation to integrate (see test program)             *
*    xi, xf    Begin, end values of variable x                            *
*    yi        Begin Value of y at x=xi                                   *
*    m         Number of points to calculate                              *
*    fi        finesse (number of intermediate points)                    *
*                                                                         *
*  OUTPUTS:                                                               *
*    t         real vector storing m results for function y               *
**************************************************************************/
void Equadif1(float *t,double xi,double xf, double yi,int m,int fi)  {
  int i,j,ni;
  double a,b,c,d,h,x,y;
  extern MaxX, MaxY;

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

  if (fi < 1) return;
  h  = (xf - xi) / fi / (m-1);
  y  = yi;
  t[1] = (float) yi;

  for (i = 1; i<m+1; i++)  {
    ni = (i - 1) * fi - 1;
    for (j = 1; j<fi+1; j++) {
      x = xi + h * (ni + j);
      a = h * fp(x,y);
      b = h * fp(x+h/2,y+a/2);
      c = h * fp(x+h/2,y+b/2);
      x = x + h;
      d = h * fp(x,y+c);
      y = y + (a + b + b + c + c + d) / 6;
	}
    t[i+1] = (float) y;
  }
  //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'=4x*(y+sqrt(y))/1+x^2",MaxX-275,70,"from x=0 to x=1");
}


void Run_equadif1() {

  // data
  xi=0; xf=1;  // x begin, x end
  yi=1;        // y'(xi)
  m=50;        // number of tabulated points of solution Y
  fi=10;       // number of intermediate calculated points
     
  Equadif1(T,xi,xf,yi,m,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_equadif1();
               
            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[] = "Teqdif1a";
    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 1", // 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 teqdif1a.cpp