///////////////////////////////////////////////////////////
// eulerg.cpp   euler with graphics
///////////////////////////////////////////////////////////

char prog[] = "eulerg";
 
// help string
extern const char helpstr[] = 
"eulerg  v(0.8) \n\n\
Num  Text\n\n\
Plot  Graphics\n\n\
(c) Richard Hall  October 2002";

#include "gwiz.h"
#include "gwiz.rc"





// ------------------------------------------------------
//  minor external additions to gwtiny

	void axes(gwin &w)
	{
		 double x0 = w.getscale(1);
		 double x1 = w.getscale(2);
		 double y0 = w.getscale(3);
		 double y1 = w.getscale(4);
		 if ((y0 < 0)&&(y1 > 0)) // x axis
		 {
			 w.move(x0,0);
			 w.line(x1,0);
		 }

		 if ((x0 < 0)&&(x1 > 0)) // y axis
		 {
			 w.move(0,y0);
			 w.line(0,y1);
		 }
	 }

    // add scale list at bottom of graph
	void pscale1(gwin &w)
	{
		 double x0 = w.getscale(1);
		 double x1 = w.getscale(2);
		 double y0 = w.getscale(3);
		 double y1 = w.getscale(4);
		 w.unclip();
		 double yp = y0 - (y1-y0)/10;
		 w.pf(x0, yp,"[%3.2f,%3.2f,%3.2f,%3.2f]",x0,x1,y0,y1);
		 w.clip();

	 }

	 // add scale labels
	 // more complicated code is required to suit all window sizes
	 // this is OK for a large window
	void pscale(gwin &w, char *xs, char* ys)
	{
		 double x0 = w.getscale(1);
		 double x1 = w.getscale(2);
		 double y0 = w.getscale(3);
		 double y1 = w.getscale(4);
		 int scan = w.getscan();
		 // convenient temporary scale
		 w.scale(0,100,0,100,100);
		 w.unclip();
		 w.pf(0,-10,"%3.2f", x0);
		 w.pf(90,-10,"%3.2f", x1);
		 w.p(45,-10,xs);
		 w.pf(-10,0,"%3.2f",y0);
		 w.pf(-10,98,"%3.2f",y1);
		 w.p(-10,50,ys);
		 // restore original scale
		 w.scale(x0,x1,y0,y1,scan);
		 w.clip();
  }
  
// ------------------------------------------------------
// a very simple ode solver

   class euler
      {

	   protected:  // allows access by derived classes
		  fptype2 fp;
		  double x0,x,y,h;
		  long i;
		  int iF; // iterate flag [1:Euler-I, 2:Euler-II]


		public:

		euler(fptype2 fpi): fp(fpi),h(0.1),iF(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; }

		inline setiF(int iFi = 1) {iF = iFi;}

		virtual void iterate(int ni = 1); 


		
       
		};

  // -------------------------------------------------------------
	void euler::iterate(int ni)   // revise (i x y) ni times
	 {
		if (iF ==1) // Euler-I
		{
		   for (int j = 0; j < ni; j++)
		   {
			y += h*fp(x,y);
			i++;
			x = x0+i*h;
		   }
		}
		else  // Euler-II
		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 euler  ----


// derived class from euler to adjoin graphics
  class eulerg : public euler

{

public:
	  eulerg(fptype2 fpi): euler(fpi) {}
	  void xplotd(gwin &w, double xa, double xb, double yi, int m);
	  void xplot (gwin &w, double yi);
	  void xplots (gwin &w, double ya, double yb, int nd);
};


void eulerg::xplotd(gwin &w, double xa, double xb, double yi, int m)

	{
		x = xa;
		h=(xb-xa)/m;
		initxy(x, yi); w.move(x, yi);
		for (int j=1;j<=m;j++)
		{iterate (); w.line(x, y);}
	};

 // eulerIIg::xplot calls eulerIIg::xplotd
 void eulerg::xplot (gwin &w, double yi)
	{
		double xa = w.getscale(1);
		double xb = w.getscale(2);
		int scan = w.getscan();
		eulerg::xplotd(w,xa,xb,yi,scan);
	};

// ------------------------------------------------------
//  variables and functions needed for this program
// ------------------------------------------------------
 
 //   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;
		}

// ----------------------------------------------------------
//  The definitions of num() and plot() declared in gwiz.h
// ----------------------------------------------------------

    void num(gwin &w, int p)

	   { w.open(p);
		 w.locate(10,90,10,90);
		 w.preframe();
		 w.scale(0,10,20,0,100);

		 euler   *e1 = new euler(f);
	     euler *e2 = new euler(f);
		 e2->setiF(2); // Euler-II
	     euler *e3 = new euler(f);
		 e3->setiF(2);



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

		w.p(1,1,"----------------------------------------------------------");
		w.p(1,2,"     x         y1           y2           y3         y");
		w.p(1,3,"----------------------------------------------------------");

	  
	    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'

		  w.pf(1,4+i,"  %.6f   %.6f    %.6f   %.6f   %.6f",x,y1,y2,y3,y);
		}
		w.p(1,15,"----------------------------------------------------------");

		delete e1,e2,e3;
		 w.close();
		 }

	 // ------------------------------------------------------------------
    
	
	void plot(gwin &w, int p)

		{ w.open(p);
	      int wb = 20;          // screen window bottom
		  if (p == 1) wb = 40;  // higher for printer
		  w.locate(20,80,wb,80);

		  w.preframe();
		  w.scale(1,2.5,0,2.5,200);
		  axes(w);
		  pscale(w,"x","y");

		  double x,y;

		  eulerg *e = new eulerg(f);
		  e->xplotd(w,1,2,1,10); // Euler-I
		  x = 2; y = e->gety();
		  w.p(2,y,"Euler-I");
		  e->setiF(2);
          e->xplotd(w,1,2,1,10); // Euler-II
		  x = 2; y = e->gety();
		  w.p(2,y,"Euler-II");
		  w.xplotd(g,1,2,200);   // Exact
		  x = 2; y = g(x);
		  w.move(x,y);
		  x = 1.8; y = y + 0.2;
		  w.line(x,y);
		  w.p(x,y, "exact");
		  
		  w.unclip();
		  w.scale(0,100,0,100,100); // for title
		  w.p(5,110,"eulerg    y' = (x^2 y^2)/5    y(1) = 1    h = 0.1");
		  delete e;
		  w.close();
		 }

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