/*******************************************************************
*        Differential equations with p variables of order 1        *
*             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++ in API style Version By J-P Moreau, Paris  *
*                      (use with teqdifp1.mak and graph2d.cpp)     *
*                                 (www.jpmoreau.fr)                *
* ---------------------------------------------------------------- *
* SAMPLE DATA:                                                     *
*                                                                  *
* Example: integrate system of equations from x=0 to x=3:          *
*              y1' = y2 + y3 - 3*y1                                *
*              y2' = y1 + y3 - 3*y2                                *
*              y3' = y1 + y2 - 3*y3                                *               
* with initial conditions: y1(0)=1, y2(0)=2 and y3(0)=-1           *
*                                                                  *
*                                                                  *
*        DIFFERENTIAL EQUATION WITH P VARIABLE OF ORDER 1          *
*              of type yi' = f(y1,y2,...,yn), i=1..n               *
*                                                                  *
*   number of variables: 3                                         *
*   begin value x      : 0                                         *
*   end value x        : 3                                         *
*   y1 value at x0     : 1                                         *
*   y2 value at x0     : 2                                         *
*   y3 value at x0     : -1                                        *
*   number of points   : 7                                         *
*   finesse            : 100                                       *
*                                                                  *
*       X            Y1           Y2                               *
* --------------------------------------                           *
*   0.000000     1.000000     2.000000                             *
*   0.500000     0.449466     0.584801                             *
*   1.000000     0.251358     0.269674                             *
*   1.500000     0.149580     0.152058                             *
*   2.000000     0.090335     0.090671                             *
*   2.500000     0.054738     0.054784                             *
*   3.000000     0.033193     0.033200                             *
*                                                                  *
*******************************************************************/
#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,i,ndata,p;
double xi,xf;
double yi[10];
float  v1[SIZE],v2[SIZE];

//Example: y1'=y2+y3-3y1, y2'=y1+y3-3y2, y3'=y1+y2-3y3
double fp(int k, double x, double *y) {
  switch(k) {
    case 0: return y[1]+y[2]-3*y[0];
    case 1: return y[0]+y[2]-3*y[1];
    case 2: return y[0]+y[1]-3*y[2];
    default: return 0;
  }
}

/***************************************************************************
*         SOLVING DIFFERENTIAL SYSTEMS WITH P VARIABLES OF ORDER 1         *
*                 of type yi' = f(y1,y2,...,yn), i=1..n                    *
* ------------------------------------------------------------------------ *
*  INPUTS:                                                                 *
*    m         number of points to display                                 *
*    xi, xf    begin, end values of variable x                             *
*    Yi        table of begin values of functions at xi                    *
*    p         number of independant variables                             *
*    fi        finesse (number of intermediate points)                     *
*                                                                          *
*  OUTPUTS:                                                                *
*    t1,t2     real vectors storing the results for first two functions,   *
*              y1 and y2.                                                  *
* ------------------------------------------------------------------------ *      
*  EXAMPLE:    y1'=y2+y3-3y1, y2'=y1+y3-3y2, y3'=y1+y2-3y3                 *
*              Exact solution :  y1 = 1/3 (exp(-4x)  + 2 exp(-x))          *
*                                y2 = 1/3 (4exp(-4x) + 2 exp(-x))          *
*                                y3 = 1/3 (-5exp(-4x)+ 2 exp(-x))          *
***************************************************************************/
void Eqdifp(float *t1,float *t2,double xi,double xf,double *Yi,int m,int p,int fi)  {
  double h,x;
  double ta[10],tb[10],tc[10],td[10],y[10],z[10];
  int    i,j,k,ni;
  char   s1[20],s2[20],s3[20],s4[20];

  extern MaxX, MaxY;  //defined in graph2d.h

  //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);
   p--;
   t1[1]=(float) Yi[0];
   t2[1]=(float) Yi[1];
   for (k=0; k<p+1; 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<p+1; k++)  y[k]=z[k];
       for (k=0; k<p+1; k++)  ta[k]=h*fp(k,x,y);
       for (k=0; k<p+1; k++)  y[k]=z[k]+ta[k]/2;
       x=x+h/2;
       for (k=0; k<p+1; k++)  tb[k]=h*fp(k,x,y);
       for (k=0; k<p+1; k++)  y[k]=z[k]+tb[k]/2;
       for (k=0; k<p+1; k++)  tc[k]=h*fp(k,x,y);
       for (k=0; k<p+1; k++)  y[k]=z[k]+tc[k];
       x=x+h/2;
       for (k=0; k<p+1; k++)  td[k]=h*fp(k,x,y);
       for (k=0; k<p+1; k++)
         z[k]=z[k]+(ta[k]+2*tb[k]+2*tc[k]+td[k])/6;
     }
     t1[i+1]=(float) z[0];
     t2[i+1]=(float) z[1];
   }

  //Draw two first solutions with legends
  strcpy(s1,"y1'=y2+y3-3y1");
  strcpy(s2,"y2'=y1+y3-3y2");
  strcpy(s3,"y3'=y1+y2-3y3");
  strcpy(s4,"from x=0 to x=3");
  CourbeXY(m,5,t1,xi,xf);
  Legendes(5," SOLUTION Y1 = F(X) ","X","Y",MaxX-150,50,s1,
	                                        MaxX-150,70,s2);
  CourbeXY(m,3,t2,xi,xf);
  Legendes(3," SOLUTION Y2 = F(X) ","X","Y",MaxX-150,90,s3,
	                                        MaxX-160,110,s4);
}

void Run_eqdifp() {
  //data section
  p=3;         // number of variables
  xi=0; xf=3;  // x begin, x end
  yi[0]=1;     // y1[xi]
  yi[1]=2;     // y2[xi]
  yi[2]=-1;    // y3[xi]
  ndata=50;    // number of tabulated points of solutions
  fi=20;       // number of intermediate calculated points

  Eqdifp(v1,v2,xi,xf,yi,ndata,p,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_eqdifp();
               
            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[] = "Teqdifp1";
    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 P 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 teqdifp1.cpp