source: Sophya/trunk/SophyaLib/Samba/circle.cc@ 1758

#include <math.h>
#include "circle.h"
//++ 
// Class        Circle
//
// include      circle.h math.h
//--
//++
//
// Links        Parents
//
//    Geometry
//    
//-- 
//++
// Titre        Constructors
//--
//++
Circle::Circle() 
//
//--
{
  UnitVector temp;
  SetCircle(temp,M_PI/2.);
}
//++
Circle::Circle(double theta, double phi, double aperture) 
//
//--
{
  UnitVector temp(theta,phi);
  SetCircle(temp,aperture);
}
//++
Circle::Circle(double x, double y, double z, double aperture) 
//
//--
{
  UnitVector temp(x,y,z);
  SetCircle(temp,aperture);
}
//++
Circle::Circle(const Vector3d& v, double aperture) 
//
//--
{
  UnitVector temp=v;
  SetCircle(temp,aperture);
}
//++
Circle::Circle(const Circle& c) 
//
//    copy constructor
//--
{
  UnitVector temp= c.Omega();
  SetCircle(temp,c._angouv);
}
//++
// Titre        Public Methods
//--
//++
void Circle::SetCircle(const UnitVector& temp, double aperture) 
//
//--
{
  _spinunitaxis= temp;

  _angouv = aperture;
  _cangouv= cos(_angouv);
  _sangouv= sin(_angouv);

  _spinaxis=_spinunitaxis*fabs(_cangouv);

  _theta =_spinunitaxis.Theta();
  _ctheta= cos(_theta);
  _stheta= sin(_theta);

  _phi =_spinunitaxis.Phi();
  _cphi= cos(_phi);
  _sphi= sin(_phi);

  _x= _spinunitaxis.X();
  _y= _spinunitaxis.Y();
  _z= _spinunitaxis.Z();
}
//++
void Circle::SetSpinAxis(double theta, double phi)
//
//--
{ 
  UnitVector temp(theta,phi);
  SetCircle(temp,_angouv);
}
//++
void Circle::SetSpinAxis(const Vector3d& u) 
//
//--
{
  UnitVector temp=u;
  SetCircle(temp,_angouv);
}
//++
void Circle::SetSpinAxis(double x, double y, double z) 
//
//--
{
  UnitVector temp(x,y,z);
  SetCircle(temp,_angouv);
}
//++
void Circle::SetApertureAngle(double aperture) 
//
//--
{ 
  SetCircle(_spinunitaxis,aperture);
}
//++
void Circle::SetApertureAngle(const Circle& c) 
//
//--
{
  SetCircle(_spinunitaxis,c._angouv);
}

//++
bool Circle::Intersection(const Circle& c, double* psi) const 
//
//    psi contains  4 values of the intersection  angles. 
//    -1 if  circles do not intersect
//    psi[0]=psi(i,j,0)
//    psi[1]=psi(i,j,1)
//    psi[2]=psi(j,i,0)
//    psi[3]=psi(j,i,1)
//--
{
  double alphak=_angouv;
  double alphal=c._angouv;
  Vector3d ok=_spinaxis;
  Vector3d ol=c._spinaxis;
  double gamma=ok.SepAngle(ol);

  if( fabs(alphak-alphal) < gamma && gamma <= (alphak+alphal) && this != &c ) 
    { 
      // then the 2 circles intersect
      double sg=sin(gamma),cg=cos(gamma);
      double sak=sin(alphak),cak=cos(alphak);
      double sal=sin(alphal),cal=cos(alphal);
      double st=sin(_theta),ct=cos(_theta);
      double stc=sin(c._theta),ctc=cos(c._theta);
      double dphi=_phi-c._phi;
      double sdphi=sin(dphi),cdphi=cos(dphi);
      double sinusk=stc*sdphi/sg,cosinusk=(ctc*st-stc*ct*cdphi)/sg;
      double sinusl=-st*sdphi/sg,cosinusl=(ct*stc-st*ctc*cdphi)/sg;
      double gammaik=scangle(sinusk,cosinusk);
      double gammail=scangle(sinusl,cosinusl);
      double omegak=acos((cal-cak*cg)/sg/sak);
      double omegal=acos((cak-cal*cg)/sg/sal);
      psi[0]=fmod(gammaik-omegak+pi2,pi2);
      psi[1]=fmod(gammaik+omegak+pi2,pi2);
      psi[2]=fmod(gammail-omegal+pi2,pi2);
      psi[3]=fmod(gammail+omegal+pi2,pi2);
      if( psi[0] > psi[1] ) 
        {
          // psi[0]=psi(i,j,0)
          // psi[1]=psi(i,j,1)
          // psi[2]=psi(j,i,0)
          // psi[3]=psi(j,i,1)
          swap(psi[0],psi[1]); 
          swap(psi[2],psi[3]); 
        }                      
      return true;
    }
  else 
    {
      psi[0] = -1.;
      psi[1] = -1.;
      psi[2] = -1.;
      psi[3] = -1.;
      return false;
    }
}
//++
UnitVector Circle::ConvToSphere(double psi) const 
//
//    Return UnitVector corresponding to a given position donnee on the circle
//--
{
  psi=mod(psi,pi2);
  double xout, yout, zout;
  double cosa=cos(_angouv);
  double sina=sin(_angouv);
  double cost=cos(_theta);
  double sint=sin(_theta);
  double cosphi=cos(_phi);
  double sinphi=sin(_phi);
  double cosp=cos(psi);
  double sinp=sin(psi);
  xout = cosa*sint*cosphi+sina*(sinphi*sinp-cost*cosphi*cosp);
  yout = cosa*sint*sinphi-sina*(cosphi*sinp+cost*sinphi*cosp);
  zout = cosa*cost+sina*sint*cosp;
  return UnitVector(xout,yout,zout);
}
//++
UnitVector Circle::TanOnCircle(double psi) const 
//
//    Return UnitVector corresponding to the tangent to the circle 
//    at given position on the circle.
//--
{
  psi=mod(psi,pi2);
  double xout, yout, zout;
  double cost=cos(_theta);
  double sint=sin(_theta);
  double cosphi=cos(_phi);
  double sinphi=sin(_phi);
  double cosp=cos(psi);
  double sinp=sin(psi);
  xout = cosp*sinphi+sinp*sint*cosphi;
  yout = -cosp*cosphi+sinp*sint*sinphi;
  zout = -sinp*cost;
  return UnitVector(xout,yout,zout);
}
//++
UnitVector Circle::EPhi(double psi) const 
//
//    Return the  vector  tangent to the sphere in the plane (xy) 
//    at a given position on the circle.
//--
{
  psi=mod(psi,pi2);
  return ConvToSphere(psi).EPhi();
}
//++
UnitVector Circle::ETheta(double psi) const 
//
//    Return the other tangent  vector( orthogonal to EPhi)-- 
//    see previous method
//--
{
  psi=mod(psi,pi2);
  return ConvToSphere(psi).ETheta();
}
//++
double Circle::SepAngleTanEPhi02PI(double psi) const
//
//    Return separation angle in [0,2Pi] at a given position on the 
//    circle and EPhi
//--
{
  psi=mod(psi,pi2);
  UnitVector pol=this->TanOnCircle(psi);
  UnitVector ephi=this->EPhi(psi);
  double angle=pol.SepAngle(ephi);
  if( pol.Z() <= 0 ) angle=pi2-angle;
  return angle;
}
//++
void Circle::Print(ostream& os) const 
//
//--
{
  os << "1 - Circle - Axe de Spin Unitaire : " << _spinunitaxis << endl;
  os << "1 - Circle - Axe de Spin : " << _spinaxis << endl;
  os << "2 - Circle - Angle d'ouverture : " << _angouv << endl;
  os << "3 - Circle - Theta,Phi : " << _theta << "," << _phi << endl;
  os << "4 - Circle - x,y,z : " << _x << "," << _y << "," << _z << endl;
}
//++
//
// inline double Theta() const
// inline double Phi() const 
// inline double ApertureAngle() const 
// inline Vector3d Omega() const 
//--
//++
// Titre        Operators
//--

Circle& Circle::operator=(const Circle& c) 
{
  if( this != &c ) 
    {
      UnitVector temp(c.Omega());
      SetCircle(temp,c.ApertureAngle());
    } 
  return *this;
}
//++
bool Circle::operator==(const Circle& c) const 
//
//--
{
  bool flag;
  if( this == &c ) flag=true;
  else flag=false;
  return flag;
}
//++
bool Circle::operator!=(const Circle& c) const 
//
//--
{
  return (bool)(1-(this->operator==(c)));
}
//
bool Circle::intersection(const Circle* c) const
{
  double alphak= _angouv;
  double alphal= c->_angouv;
  Vector3d ok= _spinaxis;
  Vector3d ol= c->_spinaxis;
  double gamma= ok.SepAngle(ol);

  if(fabs(alphak-alphal) < gamma && gamma <= (alphak+alphal) && this != c) {
    return true;
  } else {
    return false;
  }
}
//
bool Circle::intersection(const Circle& c, double* psi) const
{
  double alphak= _angouv;
  double alphal= c._angouv;
  Vector3d ok= _spinaxis;
  Vector3d ol= c._spinaxis;
  double gamma= ok.SepAngle(ol);

  if(fabs(alphak-alphal) < gamma && gamma <= (alphak+alphal) && this != &c) {

    double sgamma= sin(gamma); 
    double cgamma= cos(gamma);

    double sdphi= _sphi*c._cphi - _cphi*c._sphi; 
    double cdphi= _cphi*c._cphi + _sphi*c._sphi;

    double ssk= c._stheta*sdphi/sgamma;
    double csk= (c._ctheta*_stheta-c._stheta*_ctheta*cdphi)/sgamma;

    double ssl= -_stheta*sdphi/sgamma;
    double csl= (_ctheta*c._stheta-_stheta*c._ctheta*cdphi)/sgamma;

    double ak= atan2(ssk,csk);
    double al= atan2(ssl,csl);
    double omegak= acos((c._cangouv-_cangouv*cgamma)/sgamma/_sangouv);
    double omegal= acos((_cangouv-c._cangouv*cgamma)/sgamma/c._sangouv);

    psi[0]= fmod(ak-omegak+pi2,pi2);
    psi[1]= fmod(ak+omegak+pi2,pi2);
    psi[2]= fmod(al-omegal+pi2,pi2);
    psi[3]= fmod(al+omegal+pi2,pi2);

    if(psi[0] > psi[1]) {
      swap(psi[0],psi[1]); 
      swap(psi[2],psi[3]); 
    }                      
    return true;
    
  } else {
    psi[0] = -1.;
    psi[1] = -1.;
    psi[2] = -1.;
    psi[3] = -1.;
    return false;
  }
}