// SurfCylinder.cpp

#include "SurfCylinder.h"
#include <cassert>
#include "objstream/ObjData.hpp"
#include "objstream/ObjTable.hpp"
#include "spacegeom/SpacePath.h"
#include "trfutil/trfstream.h"
#include "trfutil/TRFMath.h"
#include "trfbase/CrossStat.h"
#include "trfbase/VTrack.h"

using std::fabs;

using namespace trf;

//**********************************************************************
// Free functions.
//**********************************************************************

namespace {

// Creator.

ObjPtr create(const ObjData& data) {
  assert( data.get_object_type() == "SurfCylinder" );
  double radius = data.get_double("radius");
  return ObjPtr( new SurfCylinder(radius) );
}

}

//**********************************************************************
// Member functions.
//**********************************************************************

// Output stream operator.

void SurfCylinder::ostr(ostream& stream) const {
  stream << begin_object;
  raw_ostr(stream);
  stream << end_object;
}

//**********************************************************************

// Raw output stream operator.

void SurfCylinder::raw_ostr(ostream& stream) const {
  stream << "Cylinder with radius " << _radius;
}

//**********************************************************************

// Equality.

bool SurfCylinder::safe_pure_equal(const Surface& srf) const {
  return is_equal(_radius, ((SurfCylinder*)&srf)->_radius);
}

//**********************************************************************

// ordering.

bool SurfCylinder::safe_pure_less_than(const Surface& srf) const {
  return _radius < ((SurfCylinder*)&srf)->_radius;
}

//**********************************************************************

// Constructor.

SurfCylinder::SurfCylinder(double radius) : _radius(radius) { }

//**********************************************************************

// Destructor.

SurfCylinder::~SurfCylinder() { }

//**********************************************************************

// Virtual constructor.

#ifdef DEFECT_NO_VIRTUAL_COVARIANCE
Surface*
#else
SurfCylinder*
#endif
              SurfCylinder::new_pure_surface() const {
  return new SurfCylinder(_radius);
}

//**********************************************************************

// Return the creator.

ObjCreator SurfCylinder::get_creator() {
  return create;
}

//**********************************************************************

// Write the object data.

ObjData SurfCylinder::write_data() const {
  ObjData data( get_type_name() );
  data.add_double( "radius", get_radius() );
  return data;
}

//**********************************************************************

// Return the crossing status for a track vector without error.

CrossStat SurfCylinder::pure_status(const VTrack& trv) const {
  // If the track surface is the same as this, return at.
  SurfacePtr psrf = trv.get_surface();
  if ( psrf==this || pure_equal(*psrf) ) return CrossStat(CrossStat::AT);
  // Otherwise extract the space point and set flags using r.
  double rtrk = trv.space_point().rxy();
  double rsrf = _radius;
  double prec = CrossStat::get_static_precision();
  if ( abs(rtrk-rsrf) < prec ) return CrossStat(CrossStat::ON);
  if ( rtrk > rsrf ) return CrossStat(CrossStat::OUTSIDE);
  return CrossStat(CrossStat::INSIDE);
}

//**********************************************************************

// Return the surface parameter.
double SurfCylinder::get_parameter(int ipar) const {
  assert( ipar == RADIUS);
  if ( ipar == RADIUS ) return _radius;
  return 0.0;
}

//**********************************************************************

// Return the difference between two track vectors.

TrackVector SurfCylinder::vec_diff(const TrackVector& vec1, 
                                   const TrackVector& vec2) const {
 TrackVector diff = vec1;
 diff -= vec2;
 int iphi = 0;
 diff(iphi) = fmod2( diff(iphi), TWOPI );
 return diff;
} 

//**********************************************************************

// Return the space point for a track.

SpacePoint SurfCylinder::space_point(const TrackVector& vec) const {
  return CylindricalPoint( _radius, vec(0), vec(1) );
}

//**********************************************************************

// Return the space vector for a track.
//     dr/ds = cos(lambda)*cos(alpha)
// r*dphi/ds = cos(lambda)*sin(alpha)
//     dz/ds = sin(lambda)

SpacePath SurfCylinder::space_vector(const TrackVector& vec,
                                     TrackSurfaceDirection dir) const {
  double r = _radius;
  double phi = vec(0);
  double z = vec(1);
  double alf = vec(2);
  double tlam = vec(3);
  double salf = sin(alf);
  double calf = cos(alf);
  double slam = tlam/sqrt(1.0+tlam*tlam);
  double clam = 1.0;
  if ( tlam != 0.0 ) clam = slam/tlam;
  double dr_ds = clam*calf;
  double r_dphi_ds = clam*salf;
  double dz_ds = slam;
  return CylindricalPath(r, phi, z, dr_ds, r_dphi_ds, dz_ds);
}

//**********************************************************************

// Return q/p.

double SurfCylinder::qoverp(const TrackVector& vec) const {
  double tlam = vec(ITLM);
  double clam = 1.0/sqrt(1.0+tlam*tlam);
  assert( clam != 0.0 );
  return vec(IQPT)*clam;
}

//**********************************************************************

// Return the direction.
// Forward is radially outward.

TrackSurfaceDirection
SurfCylinder::get_direction(const TrackVector& vec) const {
  double aalf = fabs( fmod2( vec(IALF), TWOPI ) );
  assert( aalf <= PI );
  if ( aalf <= PI2 ) return TSD_FORWARD;
  else return TSD_BACKWARD;
}

//**********************************************************************