double THelixTrack::Eval(double step, double *xyz, double *dir) const
{

  double ztep = step*fCosL;
  double teta = ztep*(fRho+0.5*ztep*fDRho);

  sgCX1   = TComplex(fX[0]  ,fX[1]);
  sgCD1   = TComplex(fP[0],fP[1])/fCosL;
  sgImTet = TComplex(0,teta);
  sgCOne  = expOne(sgImTet);			//(exp(I*Rho*L)-1)/(I*Rho*L)
  sgCf1   = sgImTet*sgCOne;
  sgCD2   = sgCD1*sgCf1+sgCD1; 			// exp(I*Fi0+I*Rho*L)
  sgCX2   = sgCD1*sgCOne*ztep;			// exp(I*Fi0)*(exp(I*Rho*L)-1)/(I*Rho)

  if (xyz) {
    xyz[0] = sgCX2.Re()+sgCX1.Re();
    xyz[1] = sgCX2.Im()+sgCX1.Im();
    xyz[2] = fX[2]+fP[2]*step;
  }
  if (dir) {    
    sgCD2/= TComplex::Abs(sgCD2);
    dir[0] = sgCD2.Re()*fCosL;
    dir[1] = sgCD2.Im()*fCosL;
    dir[2] = fP[2];
  }
  return step;
}

   X1 =        x + i*y
   D1 = cos(phi) + i*sin(phi) = Exp(i*phi)
   T  =        + i*step*curvature
   C  = (Exp(T) - 1)/T
   Cf1= T*C = Exp(T) - 1
   D2 = D1*Cf1 + D1 = D1 * Exp(T) = Exp(i*(phi+teta));
   X2 = D1*C*zstep  = D1*(Exp(T) - 1)/T *zstep

   x2 = X2 + x1