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

#include <math.h>
#include "circle.h"

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) 
{
  UnitVector temp=c.Omega();
  SetCircle(temp,c._angouv);
}

void Circle::SetCircle(const UnitVector& temp, double aperture) 
{
  _spinunitaxis=temp;
  _angouv=aperture;
  _spinaxis=_spinunitaxis*fabs(cos(_angouv));
  _theta=_spinunitaxis.Theta();
  _phi=_spinunitaxis.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 
{
  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::ConvCircleSphere(double psi) const 
{
  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 
{
  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 
{
  psi=mod(psi,pi2);
  return ConvCircleSphere(psi).EPhi();
}

UnitVector Circle::ETheta(double psi) const 
{
  psi=mod(psi,pi2);
  return ConvCircleSphere(psi).ETheta();
}

double Circle::SepAngleTanEPhi02PI(double psi) const
{
  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;
}

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)));
}

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;
}