// --------------------------------------------
//  nest.cpp     functions and nested classes
// --------------------------------------------

#include "w.h"
/*
	For the purpose of encapsulation we may choose to define
	'nested classes' within a class.
*/

 class hq  // head quarters
 { 
  public:
  hq();
  ~hq();
  inline void seta(double ai) { g->seta(ai); }
  double hg(double x);
  double hdg(double x);
  double hddg(double x);
 // ----------------------------------------------
 //  first nested class
 class gfun : public fun
	{ private:
		  double a;
	  public:
		  gfun(double ai = 1): a(ai) {}
		  void seta(double ai) { a = ai; }
		  double geta() { return a; }
		  double f(double x)
			 { return exp(-a*x*x); } };
 // ----------------------------------------------
 //  second nested class
	class dif: public fun         // derivative
		  { private: fun *f1;
				double h,h2;
			 public:
				dif(fun *f1i, double hi = 0.01)
				  { f1 = f1i; h = hi; h2 = 0.5/h; }
				double f(double x)
					{ double y =  (f1->f)(x+h);
					  y -= (f1->f)(x-h);
					  return y*h2; } };
// -----------------------------------------------

  private:
	 gfun *g;
	 dif *dg;
  };

 hq::hq(): g(new gfun()), dg(new dif(g,0.001)) {}

 hq::~hq() { delete g,dg; }

 double hq::hg(double x) { return g->f(x); }

 double hq::hdg(double x) { return dg->f(x); }

 double hq::hddg(double x)  // this is illustration is inefficient
 { dif *ddg = new dif(dg,0.001);
	double y = ddg->f(x);
	delete ddg;
	return y;
 }


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




 void main()

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

	  hq *h = new hq();

	  double a,x,y,dy,ddy;
	  x = 1;

	  a = 0.1; pl("a = ",a);
	  h->seta(a); y = h->hg(x); dy = h->hdg(x);
	  p("x = ",x);p("   y = ",y); pl("   dy = ",dy);

     nl();

	  a = 1; pl("a = ",a);
	  h->seta(a); y = h->hg(x); dy = h->hdg(x); ddy = h->hddg(x);
	  p("x = ",x);p("   y = ",y); p("   dy = ",dy); pl("   ddy = ",ddy);

	  delete h;

	  nl();

	  // we can also use the public nested classes inside hq
	  hq::gfun *g = new hq::gfun(); 
	  hq::dif *dg = new hq::dif(g);
	  // some compilers (TCW45) allow this without hq::  
	  // dif *dg = new dif(g)
	  // But it's better not to rely on this.

	  x = 2; a = 0.2;
	  g->seta(a); pl("a = ",a);
	  double z = g->f(x);
	  double dz = dg->f(x);
     p("x = ",x);p("   z = ",z); pl("   dz = ",dz);

	  delete g,dg;
	}

 /* output

	----------
    nest.cpp
   ----------
   a = 0.1
   x = 1   y = 0.904837   dy = -0.180967

   a = 1
   x = 1   y = 0.367879   dy = -0.735759   ddy = 0.735756

   a = 0.2
   x = 2   z = 0.449329   dz = -0.35946

  */