#include <iostream>
using std::cout;
using std::cerr;
using std::cin;
using std::endl;

#include <valarray>
using std::valarray;

#include <fstream>
using std::fstream;

#include <string>
using std::string;

#include <list>
using std::list;

#include <vector>
using std::vector;

#include <utility>
using std::pair;

#include "TVector3.h"
#include "TRotation.h"
#include "TRandom.h"
#include "TMatrixD.h"
#include "TVectorD.h"
#include "TDecompSVD.h"
#include "Math/Vector3D.h"
using ROOT::Math::XYZVector;

#include "Math/Point3D.h"
using ROOT::Math::XYZPoint;

#include "Math/RotationX.h"
using ROOT::Math::RotationX;
#include "Math/RotationY.h"
using ROOT::Math::RotationY;
#include "Math/RotationZ.h"
using ROOT::Math::RotationZ;
#include "Math/Rotation3D.h"
using ROOT::Math::Rotation3D;

#include "Math/Transform3D.h"
using ROOT::Math::Transform3D;

#include "PndDrcPhoton.h"
#include "PndDrcOptReflAbs.h"
#include "PndDrcSurfAbs.h"
#include "PndDrcSurfPolyFlat.h"
#include "PndDrcSurfPolyAsphere.h"
#include "PndDrcOptMatAbs.h"
#include "PndDrcOptDev.h"
#include "PndDrcUtil.h"

//----------------------------------------------------------------------
PndDrcSurfPolyAsphere::PndDrcSurfPolyAsphere()
{
  m_radius   = -999;
  m_checked  = false;
  m_conConst = 0;        // sphere
}
//----------------------------------------------------------------------
PndDrcSurfPolyAsphere* PndDrcSurfPolyAsphere::clone() const
{
  return new PndDrcSurfPolyAsphere(*this); 
}
//----------------------------------------------------------------------
void PndDrcSurfPolyAsphere::copy(const PndDrcSurfPolyAsphere& s)
{
  m_radius       = s.m_radius;
  m_conConst     = s.m_conConst;
  m_checked      = s.m_checked;
  m_trans        = s.m_trans;
  m_transInv     = s.m_transInv;
}
//----------------------------------------------------------------------
PndDrcSurfPolyAsphere::PndDrcSurfPolyAsphere(const PndDrcSurfPolyAsphere& s) : PndDrcSurfPolyFlat(s)
{
  if (s.m_verbosity>=1) cout<<"  PndDrcSurfPolyAsphere::PndDrcSurfPolyAsphere"
			    <<"(const PndDrcSurfPolyAsphere&) name,copy: "  
			    <<s.m_name<<" "<<s.m_copyNumber<<endl;
  copy(s);
}
//----------------------------------------------------------------------
PndDrcSurfPolyAsphere&  PndDrcSurfPolyAsphere::operator=(const PndDrcSurfPolyAsphere& s)
{
  if (s.m_verbosity>=1) cout<<"  PndDrcSurfPolyAsphere::operator="
			    <<"(const PndDrcSurfPolyAsphere&) name,copy: "  
			    <<s.m_name<<" "<<s.m_copyNumber<<endl;
  if (&s != this)
    {
      static_cast<PndDrcSurfPolyFlat&>((*this)) = s; // assignment of base class part.
      copy(s);
    }

  return *this;
}
//----------------------------------------------------------------------
XYZVector PndDrcSurfPolyAsphere::normal(const XYZPoint& point) const
{
  if (!m_checked) m_checked = check();// check dimensions

  // transform point back, where sphere was defined.
  XYZPoint p1(m_transInv(point));
  // p1 point on surface in unrotated and unshifted system 

  double h2   = p1.X()*p1.X()+p1.Y()*p1.Y();
  double h,root,dzdh,dhdx,dhdy;
  if (h2>1e-12)
    {
      h    = sqrt(h2); 
      root = sqrt(m_radius*m_radius - (1.0+m_conConst)*h2);

      dzdh = (2*h * (m_radius+root) + 
	      (1.0+m_conConst)*h/root * h2)/((m_radius+root)*(m_radius+root));
      dhdx = p1.X()/h;
      dhdy = p1.Y()/h;
    }
  else
    {
      dhdx=0; // eg axis parallel beam throug zero
      dhdy=0;
      dzdh=0;
    }


  // real z is R-z
  dzdh *= -1;

 
  XYZVector pa(1,0,dzdh*dhdx); // vector along slope in x dir.
  XYZVector pb(0,1,dzdh*dhdy); // vector along slope in y dir.

  XYZVector pc(pa.Cross(pb));


  pc = m_trans(pc);

  return -pc.Unit();
}
//----------------------------------------------------------------------
bool PndDrcSurfPolyAsphere::surfaceHit(PndDrcPhoton& ph, 
				    XYZPoint&  pos_new, 
				    double&    path_length) const
{
  const double eps = 0.001;
  if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit: test of "
			  <<name()<<endl;  

  if (!m_checked) m_checked = check();// check dimensions

  // transform back in sphere definitions space
  XYZPoint  pos(ph.position());
  XYZVector dir(ph.direction());

  if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit: original photon: \n"
  		  <<" pos: "<<pos<<" dir: "<<dir<<endl;

  dir = m_transInv*dir;
  pos = m_transInv*pos;

  if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit: transformed photon: \n"
  		  <<" pos: "<<pos<<" dir: "<<dir<<endl;

  // compute hit position
  //
  // line : x = pos + lambda*dir
  //
  // h2 = x*x+y*y = (posx+lambda dirx)^2 + (posy+lambda diry)^2
  //
  // asphere: 
  //
  // R - h2/(R+sqrt[R2-(1+k)h2]) = pz + lambda dz
  //
  // R - h2*(R-sqrt[])/(R2-sqrt[]^2) = pz + lambda dz
  //
  // R - pz = h2*(R-sqrt[])/(R2-R2+(1+k)h2) + lambda dz
  //
  // R - pz = (R-sqrt[])/(1+k) + lambda dz      
  //
  // (1+k)(R - pz - lambda dz) = R-sqrt
  //
  // kR - (1+k) pz - (1+k) lambda dz = -sqrt
  //
  // with a=kR - (1+k) pz and b=(1+k) dz
  //
  // a - b lambda = -sqrt[]
  //
  // b2 lambda2 - 2 ab lambda + a2 = R2 - (1+k)h2
  //
  // insert h2
  //
  // b2 lambda2 - 2ab lambda + a2 - R2 = 
  // -(1+k)*([px2 + 2 lambda px dx + lambda2 dx2]+[py2 + 2 lambda py dy + lambda2 dy2])
  // 
  //
  // with c1 = -(1+k)px2
  //      c2 = -(1+k)2 px dx
  //      c3 = -(1+k)dx2
  // with d1 = -(1+k)py2
  //      d2 = -(1+k)2 py dy
  //      d3 = -(1+k)dy2
  //
  // b2 lambda2 - 2ab lambda + a2 - R2 =
  // c1 + lambda c2 + lambda2 c3 + d1 + lambda d2 + lambda2 d3
  //
  // lambda2 (b2-c3-d3) - lambda (2ab+c2+d2) +(a2-R2-c1-d1) = 0
  //

  long double a  = m_conConst*m_radius - (1.0+m_conConst)*pos.Z();
  long double b  = (1.0+m_conConst)*dir.Z();
  long double c1 = -(1.0+m_conConst)*pos.X()*pos.X();
  long double c2 = -(1.0+m_conConst)*2*pos.X()*dir.X();
  long double c3 = -(1.0+m_conConst)*dir.X()*dir.X();
  long double d1 = -(1.0+m_conConst)*pos.Y()*pos.Y();
  long double d2 = -(1.0+m_conConst)*2*pos.Y()*dir.Y();
  long double d3 = -(1.0+m_conConst)*dir.Y()*dir.Y();

  if (m_verbosity>=4) 
    {
      cout<<" PndDrcSurfPolyAsphere::surfaceHit: a  ="<<a<<endl; 
      cout<<" PndDrcSurfPolyAsphere::surfaceHit: b  ="<<b<<endl; 
      cout<<" PndDrcSurfPolyAsphere::surfaceHit: c1 ="<<c1<<endl; 
      cout<<" PndDrcSurfPolyAsphere::surfaceHit: c2 ="<<c2<<endl; 
      cout<<" PndDrcSurfPolyAsphere::surfaceHit: c3 ="<<c3<<endl; 
      cout<<" PndDrcSurfPolyAsphere::surfaceHit: d1 ="<<d1<<endl; 
      cout<<" PndDrcSurfPolyAsphere::surfaceHit: d2 ="<<d2<<endl; 
      cout<<" PndDrcSurfPolyAsphere::surfaceHit: d3 ="<<d3<<endl;     
    }
  
  //
  // for k=-1 , parabola: c1,c2,c3,d1,d2,d3,b=0 a=-R
  //
  // R - h2/(R+sqrt[R2]) = pz + lambda dz
  // 2R2 - h2 = 2R pz + lamda 2R dz
  // 2R2 - [(posx+lambda dirx)^2 + (posy+lambda diry)^2] = 2R pz + lambda 2R dz
  // 2R2 - [px2+py2 + lambda 2 (pxdx+pydy) + lambda2 (dx2 + dy2) ] 
  //                      = 2R pz + lambda 2R dz
  // [px2+py2 + lambda 2 (pxdx+pydy) + lambda2 (dx2 + dy2) ] 
  //                      = -2R pz - lambda 2R dz +2R2
  // 
  // lambda2 (dx2+dy2) + lambda 2 [(pxdx+pydy)+R dz] + px2+py2+2Rpz-2R2 = 0
  //            h2                    a                  b

  long double lambda;

  // special treatment for parabola
  if (fabs(m_conConst+1)<1.0e-4)
    {// parabola
      long double h2 = dir.X()*dir.X()+dir.Y()*dir.Y();
      long double a = 2*(pos.X()*dir.X()+pos.Y()*dir.Y()+m_radius*dir.Z());
      long double b = 2*m_radius*(pos.Z()-m_radius) + pos.X()*pos.X() + pos.Y()*pos.Y();
      if (fabs(h2)<1.0e-9)
	{
	  if (a==0)
	    {
	      if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit no lambda"
				      <<endl;  
	      return false;
	    }
	  else
	    {
	      lambda = -b/a;
	      XYZPoint pos_new1     = pos + lambda*dir;
	      XYZPoint pos_check(pos_new1);
	      pos_check.SetZ(0);
	      bool hit = withinSurface(pos_check);  

	      if (hit && lambda>0)
		{
		  pos_new = pos_new1;
		  path_length = lambda;
		  if (!onSurface(pos_new)) 
		    {
		      if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit"
					      <<" calc. point not on surface,"
					      <<" wrong lambda"<<endl;
		      return false; // wrong lambda
		    }
		  pos_new = m_trans*pos_new;
		  if (m_pixel) ph.setFate(Drc::PhotMeasured);
		  if (m_verbosity>=4) 
		    cout<<" PndDrcSurfPolyAsphere::surfaceHit hit lambda="
				      <<lambda<<endl;  
		  return true;
		}
	      else
		{
		  if (m_verbosity>=4) 
		    cout<<" PndDrcSurfPolyAsphere::surfaceHit hit,lambda="
			<<hit<<" "<<lambda<<endl;  
		  return false;;
		}
	    }
	}
      else
	{// h2>0
	  long double p = a/h2;
	  long double q = b/h2;
	  long double root2 = p*p/4-q;

	  if (m_verbosity>=4)
	    cout<<" PndDrcSurfPolyAsphere::surfaceHit p,q,root2: "
		<<p<<" "<<q<<" "<<root2<<endl;  
	  if (root2<0)
	    {
	      if (m_verbosity>=4)
		cout<<" PndDrcSurfPolyAsphere::surfaceHit root2<0"<<endl;  
	      return false;
	    }
	  else
	    {// root2>0
	      long double lambda1 = -p/2 + sqrt(root2);
	      long double lambda2 = -p/2 - sqrt(root2);
	      if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit lambda_1,2 ="
				      <<lambda1<<" "<<lambda2<<endl;  
	      if (lambda1>eps && lambda2>eps)
		{ // check both solutions
		  XYZPoint pos_check     = pos + lambda1*dir;
		  pos_check.SetZ(0);
		  bool hit1 = withinSurface(pos_check);  
		  pos_check              = pos + lambda2*dir;
		  pos_check.SetZ(0);
		  bool hit2 = withinSurface(pos_check);
		  if (hit1 && hit2)
		    {
		      // take closest
		      lambda = (lambda1<lambda2) ? lambda1 : lambda2;
		      pos_new = pos + lambda*dir;
		      if (!onSurface(pos_new)) 
			{
			  if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit"
						  <<" calc. point not on surface,"
						  <<" wrong lambda"<<endl;
			  return false; // wrong lambda
			}
		      path_length = lambda;
		      pos_new = m_trans*pos_new;
		      if (m_pixel) ph.setFate(Drc::PhotMeasured);
		      if (m_verbosity>=4) 
			{
			  cout<<" PndDrcSurfPolyAsphere::surfaceHit lambda_1,2="
			      <<lambda1<<" "<<lambda2<<endl;  
			  cout<<"                                took lambda="
			      <<lambda<<endl;
			}
		      return true;
		    }
		  else if (hit1)
		    {
		      pos_new = pos + lambda1*dir;
		      if (!onSurface(pos_new)) 
			{
			  if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit"
						  <<" calc. point not on surface,"
						  <<" wrong lambda"<<endl;
			  return false; // wrong lambda
			}
		      path_length = lambda1;
		      pos_new = m_trans*pos_new;
		      if (m_pixel) ph.setFate(Drc::PhotMeasured);
		      if (m_verbosity>=4) 
			cout<<" PndDrcSurfPolyAsphere::surfaceHit hit1 lambda_1="
					      <<lambda1<<endl;  
		      return true;
		    }
		  else if (hit2)
		    {
		      pos_new = pos + lambda2*dir;
		      if (!onSurface(pos_new)) 
			{
			  if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit"
						  <<" calc. point not on surface,"
						  <<" wrong lambda"<<endl;
			  return false; // wrong lambda
			}
		      path_length = lambda2;
		      pos_new=m_trans*pos_new;
		      if (m_pixel) ph.setFate(Drc::PhotMeasured);
		      if (m_verbosity>=4) 
			cout<<" PndDrcSurfPolyAsphere::surfaceHit hit2 lambda_1="
			    <<lambda2<<endl;  
		      return true;
		    }
		  else
		    {
		      if (m_verbosity>=4) 
			cout<<" PndDrcSurfPolyAsphere::surfaceHit "
			    <<"no hits for both lambdas"<<endl;  
		      return false;
		    }
		}
	      else if (lambda1>eps && lambda2<=eps)
		{
		  lambda = lambda1;
		}
	      else if (lambda1<=eps && lambda2>eps)
		{
		  lambda = lambda2;
		}
	      else
		{
		  //cout<<" both lambdas <0"<<endl;//###
		  if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit "
					  <<"both lambdas < 0 "
					  <<lambda1<<" "<<lambda2<<endl;  
		  return false;
		}
	    }  
	}
    } // end parabola
  else if (fabs(b*b-c3-d3)<1.0e-12)
    { // otherwise the too small argument will cause trouble in
      //                        the next clause
    
      lambda = (a*a-m_radius*m_radius-c1-d1)/(2*a*b+c2+d2);
      if (lambda<eps) 
	{
	  if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit lambda="
				  <<lambda<<" <eps no hit"<<endl;  
	  return false;
	}
      if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit lambda="
			      <<lambda<<endl;
    }
  else
    {// most general case
      long double p  = -(2*a*b+c2+d2)/(b*b-c3-d3);
      long double q  = (a*a-m_radius*m_radius-c1-d1)/(b*b-c3-d3);
      long double root2 = p*p/4 - q;
      if (root2<0) 
	{
	  return false; // photon does not hit asphere 
	  // eg. photon behind asphere perpendicular to z-axis.
	  cerr<<" *** PndDrcSurfPolyAsphere::surfaceHit: root<0 : "<<root2<<endl;
	}
      long double lambda1 = -p/2 + sqrt(root2);
      long double lambda2 = -p/2 - sqrt(root2);

      if (m_verbosity>=4) 
	{
	  cout<<" PndDrcSurfPolyAsphere::surfaceHit lambda_1,2="
			      <<lambda1<<" "<<lambda2<<endl; 
	  cout<<"                     dir1 = "<< lambda1*dir<<endl;
	  cout<<"                     dir2 = "<< lambda2*dir<<endl;
	}


      if (lambda1>eps && lambda2>eps)
	{ // check both solutions
	  XYZPoint pos_check   = pos + lambda1*dir;
	  pos_check.SetZ(0);
	  bool hit1 = withinSurface(pos_check);  
	  pos_check            = pos + lambda2*dir;
	  pos_check.SetZ(0);
	  bool hit2 = withinSurface(pos_check);
	  if (hit1 && hit2)
	    {
	      // take closest
	      lambda = (lambda1<lambda2) ? lambda1 : lambda2;
	      pos_new = pos + lambda*dir;
	      if (!onSurface(pos_new)) 
		{
		  if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit"
					  <<" calc. point not on surface,"
					  <<" wrong lambda"<<endl;
		  return false; // wrong lambda
		}
	      path_length = lambda;
	      pos_new = m_trans*pos_new;
	      if (m_pixel) ph.setFate(Drc::PhotMeasured);
	      if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit lambda_1,2="
				      <<lambda1<<" "<<lambda2<<endl;  
	      return true;
	    }
	  else if (hit1)
	    {
	      pos_new = pos + lambda1*dir;
	      if (!onSurface(pos_new)) 
		{
		  if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit"
					  <<" calc. point not on surface,"
					  <<" wrong lambda"<<endl;
		  return false; // wrong lambda
		}
	      path_length = lambda1;
	      pos_new = m_trans*pos_new;
	      if (m_pixel) ph.setFate(Drc::PhotMeasured);
	      if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit lambda_1,2="
				      <<lambda1<<" "<<lambda2<<endl;  
	      return true;
	    }
	  else if (hit2)
	    {
	      pos_new = pos + lambda2*dir;
	      if (!onSurface(pos_new)) 
		{
		  if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit"
					  <<" calc. point not on surface,"
					  <<" wrong lambda"<<endl;
		  return false; // wrong lambda
		}
	      path_length = lambda2;
	      pos_new = m_trans*pos_new;
	      if (m_pixel) ph.setFate(Drc::PhotMeasured);
	      if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit lambda_1,2="
				      <<lambda1<<" "<<lambda2<<endl;  
	      return true;
	    }
	  else
	    {
	      if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit "
				      <<"no hits for both lambdas"<<endl;  
	      return false;
	    }
	}
      else if (lambda1>eps && lambda2<=eps)
	{
	  lambda = lambda1;
	}
      else if (lambda1<=eps && lambda2>eps)
	{
	  lambda = lambda2;
	}
      else
	{
	  //cout<<" both lambdas <0"<<endl;//###
	      if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit "
				      <<"both lambdas < 0 "
				      <<lambda1<<" "<<lambda2<<endl;  
	  return false;
	}
    }

  path_length = lambda;
  pos_new     = pos + path_length*dir;

  if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit:"<<name()
			  <<" only one lambda survived"
			  <<endl;  

  if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit pos_new before trans. ="
			  <<pos_new<<endl;  
  if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit zVal of pos_new ="
			  <<zVal(pos_new.X(),pos_new.Y())<<endl;  



  if (!onSurface(pos_new)) 
    {
      if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit"
			      <<" calc. point not on surface,"
			      <<" wrong lambda"<<endl;
      return false; // wrong lambda
    }


  // check if pos_new is within polynom
  XYZPoint pos_check(pos_new);
  pos_check.SetZ(0);

 if (withinSurface(pos_check)) 
    {
      pos_new = m_trans*pos_new;
      if (m_pixel) ph.setFate(Drc::PhotMeasured);
      if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit: lambda -> hit"
			      <<endl;  
      return true;
    }

  if (m_verbosity>=4) cout<<" PndDrcSurfPolyAsphere::surfaceHit "
			  <<name()<<": false, not within surface bounds"<<endl;  

  return false;

} 
//----------------------------------------------------------------------
bool PndDrcSurfPolyAsphere::check() const
{
  if (m_radius<=0)
    {
      cerr<<" ***  PndDrcSurfPolyAsphere::check: "
	  <<m_name<<": wrong or not set radius="
	  <<m_radius<<endl;
      exit(EXIT_FAILURE);
    }
  for (unsigned int i=0; i<m_p.size(); i++)
    {
      if (zVal(m_p[i].X(),m_p[i].Y())<0)
	{
	  cerr<<" ***  PndDrcSurfPolyAspher::check: "
	      <<m_name<<": point outside sphere projection:"
	      <<" m_p["<<i<<"]="
	      <<m_p[i].X()<<" "
	      <<m_p[i].Y()<<" "
	      <<m_p[i].Z()<<endl;
	  cerr<<"      z = "<<zVal(m_p[i].X(),m_p[i].Y())<<endl;
	  exit(EXIT_FAILURE);
	  //  return false;
	}
    }
  for (unsigned int i=0; i<m_p.size(); i++)
    {
      if (m_p[i].Z() != 0)
	{
	  cerr<<" ***  PndDrcSurfPolyAsphere::check: "
	      <<m_name<<": z component of point != 0 "
	      <<" m_p["<<i<<"]="
	      <<m_p[i].X()<<" "
	      <<m_p[i].Y()<<" "
	      <<m_p[i].Z()<<endl;
	  exit(EXIT_FAILURE);
	  //return false;
	}
    }

  return true;
}
//----------------------------------------------------------------------
XYZPoint PndDrcSurfPolyAsphere::limitingPoint(unsigned int i)
{
  if (i>(m_p.size()-1) || i<0)
    {
      cerr<<" *** PndDrcSurfPolyAsphere::limitingPoint: index out of range, i="<<i<<endl;
      exit(EXIT_FAILURE);
    }
  XYZPoint p1(m_p[i]);
  
  p1.SetZ(zVal(p1.X(),p1.Y()));

  p1 = m_trans*p1;
  return p1;  
}
//----------------------------------------------------------------------
XYZPoint PndDrcSurfPolyAsphere::centerPoint()
{
  return m_trans*XYZPoint(0,0,m_radius);
}
//----------------------------------------------------------------------
void PndDrcSurfPolyAsphere::print(fstream& stream) const
{
  if (!m_checked) m_checked = check();// check dimensions

  // print points
  unsigned int isize = m_p.size();

  if (isize < 3)
    {
      cerr<<" *** PndDrcSurfPolyAsphere::surfaceHit: surface has only "<<m_p.size()
	  <<" points. Aborting. Name="<<name()<<endl;
      exit(EXIT_FAILURE);
    }



  for (double frac=1.0; frac>0.19; frac-=0.2)
    {// several (5) lines
      

      stream<<" TPolyLine3D *l = new TPolyLine3D("<<20*isize+1<<");"<<endl;

      XYZPoint m_p1;
      int icnt=0;

      for (unsigned int i=0; i<isize; i++)
	{
	  for (unsigned int j=0; j<20; j++)
	    {
	      unsigned int next = i+1;
	      if (next==isize) next=0;


	      m_p1 = (m_p[i]*frac) + j/20.0*(m_p[next]-m_p[i])*frac;
	      m_p1.SetZ(zVal(m_p1.X(),m_p1.Y())); // m_p1.Z()=0...
	  
	      m_p1 = m_trans*m_p1;

	      stream<<" l->SetPoint("<<(icnt++)<<","
		    <<m_p1.X()<<","
		    <<m_p1.Y()<<","
		    <<m_p1.Z()<<");"<<endl;
	    }
	}
  
      m_p1 = m_p[0]*frac;
      m_p1.SetZ(zVal(m_p1.X(),m_p1.Y())); // m_p1.Z()=0...

      m_p1 = m_trans*m_p1;  
      stream<<" l->SetPoint("<<(icnt++)<<","
	    <<m_p1.X()<<","
	    <<m_p1.Y()<<","
	    <<m_p1.Z()<<");"<<endl;
      stream<<" l->SetLineColor("<<m_printColor<<");"<<endl;
      stream<<" l->Draw();"<<endl;
  
    }

}
//----------------------------------------------------------------------
double PndDrcSurfPolyAsphere::zVal(double x, double y) const
{
  double h2  = x*x + y*y;
  double den = m_radius*(1.0+sqrt(1.0-(1.0+m_conConst)*h2/(m_radius*m_radius)));
  return m_radius - (h2/den);
}
//----------------------------------------------------------------------
bool PndDrcSurfPolyAsphere::onSurface(const XYZPoint& p) const
{
  double x = p.x();
  double y = p.y();
  double z = p.z();

  double zz = zVal(x,y);

  return (fabs(z-zz)<0.1) ? true : false;

}
//----------------------------------------------------------------------
void PndDrcSurfPolyAsphere::addTransform(const Transform3D& trans)
{
  if (verbosity()>=3) cout<<"    PndDrcSurfCyl::addTransform() name="<<m_name<<endl;
  m_trans    = trans * m_trans;
  m_transInv = m_trans.Inverse();
}