// -------------------------------------------
//  euler.cpp   a very simple ode engine
// -------------------------------------------

#include "w.h"

 // de engine

	  class euler
      {

		public:

		euler(fptype2 fpi): fp(fpi),h(0.1) {}
		~euler() {}


		inline void initxy(double xi, double yi)
		  { euler::x = x0 = xi; euler::y = yi; euler::i = 0; }

	   inline long geti()
		  { return(i); }

	   inline double getx()
	      { return(x); }

	   inline double gety()
	      { return(y); }

		inline void seth(double hi)  //  sets h,h2
		  { h = hi; }

		virtual void iterate(int ni = 1); // set default here


		protected:  // allows access by derived classes
		  fptype2 fp;
		  double x0,x,y,h;
		  long i;

       
		};

  // -------------------------------------------------------------
	void euler::iterate(int ni)   // revise (i x y) ni times
	 {
		 for (int j = 0; j < ni; j++)
		 {
			y += h*fp(x,y);
			i++;
			x = x0+i*h;
		 }
	 }

 // -----------------------------------   end of class euler  ----

 // a slightly improved ode engine derived from euler
 class eulerII:  public euler
 {
    public:
	 eulerII(fptype2 fpi): euler(fpi) { }
	 void iterate(int ni);  // declaration of revised iterate
 };

  void eulerII::iterate(int ni = 1) // redefinition of euler::iterate
	 {

		for (int j = 0; j < ni; j++)
		{
		  double k1,k2,y1;
		  k1 = fp(x,y);       // near slope
		  y1 = y + h*k1;      // eulerI result
		  i++;
		  x = x0+i*h;         // new x
		  k2 = fp(x,y1);      // far slope
		  y += 0.5*h*(k1+k2); // use average slope
		}
	 }


 // ------------------------------ end of class eulerII -----------

 //   Top level test function (its name is an fptype2)
		double A = 0.2;  // global parameter

		double f(double x, double y)
		 {
			 double z = x*y;
			 return z*z*A;
		 }

 //    Exact solution starting at (x,y) = (1,1)
 //    for the ode  y'  =  A (xy)^2

	  double g(double x)  // x < ((3+A)/A)^(1/3)
		{
			 double z = x*x*x;
			 z = 1.0 + A*(1-z)/3.0;
			 return 1/z;
		}

// ---------------------------------------------------


 void main()

	{
	  nl(0);
	  banner("euler.cpp");


	  euler   *e1 = new euler(f);
	  eulerII *e2 = new eulerII(f);
	  eulerII *e3 = new eulerII(f);



	  int n = 10;   // number of iterations
	  int n3 = 100; // extra iterations for e3
	  double x0 = 1, x1 = 2; // range of integration
	  double y0 = 1; // initial y-value
	  double h = (x1-x0)/n;
	  double h3 = h /n3; // step size for e3


	  e1->seth(h);
	  e2->seth(h);
	  e3->seth(h3);

	  double x = x0,y = y0;    // 'exact' values
	  double y1,y2,y3;         // approximations
     y1 = y2 = y3 = y0;

	  e1->initxy(x0,y0);
	  e2->initxy(x0,y0);
	  e3->initxy(x0,y0);

		pl("----------------------------------------------------------");
		pl("     x         y1           y2           y3         y");
		pl("----------------------------------------------------------");

	  
	  int i;
	  
	  for (i = 0; i < n; i++)
	  {
		  e1->iterate();   y1 = e1->gety();
		  e2->iterate();   y2 = e2->gety();
		  e3->iterate(n3); y3 = e3->gety();
		  x = e1->getx();   y = g(x);  // 'exact'

		  pfl("  %.6f   %.6f    %.6f   %.6f   %.6f",x,y1,y2,y3,y);
		}
		pl("----------------------------------------------------------");

		
		wait();

}

 /* output -----------------------------------------------------
	-----------
    euler.cpp
   -----------
   ----------------------------------------------------------
        x         y1           y2           y3         y
   ----------------------------------------------------------
     1.100000   1.020000    1.022589   1.022565   1.022565
     1.200000   1.045178    1.051054   1.051009   1.051009
     1.300000   1.076639    1.086779   1.086720   1.086720
     1.400000   1.115818    1.131621   1.131563   1.131563
     1.500000   1.164624    1.188145   1.188119   1.188119
     1.600000   1.225660    1.260015   1.260081   1.260081
     1.700000   1.302574    1.352653   1.352936   1.352936
     1.800000   1.400644    1.474444   1.475216   1.475216
     1.900000   1.527769    1.639075   1.640958   1.640958
     2.000000   1.696289    1.870462   1.875000   1.875000
   ----------------------------------------------------------

  */