Index: trunk/SophyaLib/Samba/Samba.o.list
===================================================================
--- trunk/SophyaLib/Samba/Samba.o.list	(revision 262)
+++ trunk/SophyaLib/Samba/Samba.o.list	(revision 262)
@@ -0,0 +1,12 @@
+bruit.o
+circle.o
+gorfph.o
+localmap.o
+longlat.o
+mlobe.o
+scan.o
+spheregorski.o
+spherethetaphi.o
+unitvector.o
+utilgeom.o
+vector3d.o
Index: trunk/SophyaLib/Samba/circle.cc
===================================================================
--- trunk/SophyaLib/Samba/circle.cc	(revision 262)
+++ trunk/SophyaLib/Samba/circle.cc	(revision 262)
@@ -0,0 +1,208 @@
+#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;
+}
Index: trunk/SophyaLib/Samba/circle.h
===================================================================
--- trunk/SophyaLib/Samba/circle.h	(revision 262)
+++ trunk/SophyaLib/Samba/circle.h	(revision 262)
@@ -0,0 +1,73 @@
+#ifndef CIRCLE_H_SEEN
+#define CIRCLE_H_SEEN
+
+#include <math.h>
+#include <iostream.h>
+#include "vector3d.h"
+#include "unitvector.h"
+#include "utilgeom.h"
+
+class Circle 
+{
+  
+ public:
+
+  Circle();
+  Circle(double theta, double phi, double aperture);
+  Circle(double x, double y, double z, double aperture);
+  Circle(const Vector3d& v, double aperture=M_PI/2.);
+  Circle(const Circle& c);
+  virtual ~Circle() {}
+
+  void SetCircle(const UnitVector&, double);
+  void SetSpinAxis(double theta, double phi);
+  void SetSpinAxis(const Vector3d&);
+  void SetSpinAxis(double x, double y, double z);
+  void SetApertureAngle(double aperture);
+  void SetApertureAngle(const Circle&);
+  
+  // psi contient les 4 valeurs des angles d'intersection. -1 si les cercles ne se coupent pas
+  // voir la numerotation dans le .cc
+  bool Intersection(const Circle&, double* psi) const;
+
+  // donne le UnitVector correspondant a une position donnee sur le cercle
+  UnitVector ConvCircleSphere(double psi) const;
+
+  // donne le UnitVector correspondant la tangente au cercle en une position donnee sur le cercle
+  UnitVector TanOnCircle(double psi) const;
+
+  // donne le vecteur tangent dans le plan (xy) a la sphere en une position donnee sur le cercle
+  UnitVector EPhi(double psi) const;
+
+  // donne l'autre vecteur tangent (orthogonal a EPhi)
+  UnitVector ETheta(double psi) const;
+
+   // donne l'angle de separation dans [0,2Pi] en une position donnee sur le cercle et EPhi
+  double SepAngleTanEPhi02PI(double psi) const;
+
+  double Theta() const {return _spinaxis.Theta();}
+  double Phi() const {return _spinaxis.Phi();}
+  double ApertureAngle() const {return _angouv;}
+  Vector3d Omega() const {return _spinaxis;}
+  virtual void Print(ostream&) const;
+  Circle& operator=(const Circle&);
+  bool operator==(const Circle&) const;
+  bool operator!=(const Circle&) const;
+
+ private:
+
+  Vector3d _spinaxis;
+  UnitVector _spinunitaxis;
+  double _angouv;
+  double _theta,_phi;
+  double _x,_y,_z;
+
+};
+
+inline ostream& operator<<(ostream& s, const Circle& c) 
+{  
+  c.Print(s);  
+  return s;  
+}
+
+#endif
Index: trunk/SophyaLib/Samba/longlat.cc
===================================================================
--- trunk/SophyaLib/Samba/longlat.cc	(revision 262)
+++ trunk/SophyaLib/Samba/longlat.cc	(revision 262)
@@ -0,0 +1,46 @@
+#include "longlat.h"
+
+LongLat::LongLat()
+{
+  _lon=0.;
+  _lat=0.;
+}
+
+LongLat::LongLat(double longitude, double latitude)
+{
+  _lon=mod(longitude,pi2);
+  if( latitude == pi_over_2 ) _lat=latitude;
+  else _lat=mod(latitude+pi_over_2,M_PI)-pi_over_2; // dans [-Pi/2,Pi/2]
+}
+
+LongLat::LongLat(double x, double y, double z)
+{
+  double norm=sqrt(x*x+y*y+z*z);
+  double theta;
+  if( norm != 0. ) 
+    {
+      theta=acos(z/norm); // dans [0,Pi]
+      if( mod(theta,M_PI) == 0. ) _lon=0.; // on est sur +-Oz, le vecteur z est en phi=0
+      //      else _lon=acos(x/sin(theta)/norm)+M_PI*(y<0);
+      else _lon=scangle(y/sin(theta)/norm,x/sin(theta)/norm);
+    }
+  else // vecteur nul
+    {
+      theta=0.;
+      _lon=0.;
+    }
+  _lat=pi_over_2-theta;
+}
+
+void LongLat::Set(double longitude, double latitude)
+{
+  _lon=mod(longitude,pi2);
+  if( latitude == pi_over_2 ) _lat=latitude;
+  else _lat=mod(latitude+pi_over_2,M_PI)-pi_over_2; // dans [-Pi/2,Pi/2]
+}
+
+void LongLat::Print(ostream& os) const
+{
+  os << "LongLat : longitude = " << _lon << " phi   = " << this->Phi() << endl;
+  os << "LongLat : latitude    = " << _lat << " theta = " << this->Theta() << endl;
+}
Index: trunk/SophyaLib/Samba/longlat.h
===================================================================
--- trunk/SophyaLib/Samba/longlat.h	(revision 262)
+++ trunk/SophyaLib/Samba/longlat.h	(revision 262)
@@ -0,0 +1,47 @@
+#ifndef LONGLAT_H_SEEN
+#define LONGLAT_H_SEEN
+
+#include <math.h>
+#include <iostream.h>
+#include <stdio.h>
+#include <string.h>
+#include "utilgeom.h"
+
+/*
+  passage des coordonnees longitude, latitude en theta,phi
+  longitude=phi
+  latitude=Pi/2-theta
+*/
+
+class LongLat
+{
+  
+ public:
+  
+  LongLat();
+  LongLat(double longitude, double latitude);
+  LongLat(double x, double y, double z);
+  LongLat(const LongLat&);
+  
+  void Set(double longitude, double latitude);
+  double Longitude() const {return _lon;}
+  double Phi() const {return _lon;}
+  double Latitude() const {return _lat;}
+  double Theta() const {return pi_over_2-_lat;}
+
+  void Print(ostream& os) const;
+  
+ private:
+  
+  double _lon; // lambda
+  double _lat; // beta
+  
+};
+
+inline ostream& operator<<(ostream& s, const LongLat& ll) 
+{  
+  ll.Print(s);  
+  return s;  
+}
+
+#endif
Index: trunk/SophyaLib/Samba/mlobe.cc
===================================================================
--- trunk/SophyaLib/Samba/mlobe.cc	(revision 254)
+++ trunk/SophyaLib/Samba/mlobe.cc	(revision 262)
@@ -1,4 +1,4 @@
 #include "mlobe.h"
-#include "unitvec.h"
+#include "unitvector.h"
 
 #include "timing.h"
@@ -113,6 +113,6 @@
   dx = -(mDx[k]*sphi+mDy[k]*ctet*cphi);
   dy = mDx[k]*cphi-mDy[k]*ctet*sphi;
-  vc.setxyz(pd.x+dx, pd.y+dy, pd.z+dz);
-  mTC[k] = vc.Teta();
+  vc.Setxyz(pd.X()+dx, pd.Y()+dy, pd.Z()+dz);
+  mTC[k] = vc.Theta();
   mPC[k] = vc.Phi();
   if (mTC[k] < 0.)  printf("ERREUR !!!! mTC[%d] = %g \n", k, mTC[k]);
Index: trunk/SophyaLib/Samba/unitvec.h
===================================================================
--- trunk/SophyaLib/Samba/unitvec.h	(revision 254)
+++ 	(revision )
@@ -1,57 +1,0 @@
-#ifndef UNITVEC_H_SEEN
-#define UNIITVEC_H_SEEN
-
-#include <math.h>
-#include <iostream.h>
-
-
-
-//  Classe pour faire des vecteurs normes x,y,z a partir de Teta,Phi
-//  et l inverse.  On peut calculer le produit scalaire et les vecteurs 
-//  perpendiculaires.
-//  Exemple:
-//  UnitVector u1(0.5, 0.7);   // Vers Teta=0.5 Phi = 0.7 Rad.
-//  UnitVector u2(1.,2.,1.);   // Vecteur norme parallele a (x=1., y=2., z=1.)
-//  u2.Teta();                 // Direction teta de u2
-//  u1.psc(u2);                // u1 prod. scalaire u2
-//  UnitVector vp1 = u2.vperp();   // Vecteur perpendiculaire a u1 , meme phi (ou a 180 deg)
-//  UnitVector vp2 = u2.vperp2();  // Vecteur perpendiculaire a u1 , a phi+90 deg
-
-//  
-class UnitVector {
-public:
-  UnitVector(double teta=0., double phi=0.) 
-    {  set(teta, phi); }
-  UnitVector(double vx, double vy, double vz) 
-      { setxyz(vx, vy, vz); }
-
-  inline void set(double teta, double phi) 
-    { t=teta; p=phi; z = cos(teta); x = sin(teta)*cos(phi);  y =  sin(teta)*sin(phi); }
-  inline void setxyz(double vx, double vy, double vz) 
-    { double r = sqrt(vx*vx+vy*vy+vz*vz); x = vx/r; y = vy/r; z = vz/r;
-    t = acos(z); p = atan2(y,x); }
-
-  inline double psc(UnitVector& u) { return(x*u.x+y*u.y+z*u.z); }
-
-  inline double Teta() { return(t); }
-  inline double Phi() { return(p); }
-  inline double X() { return(x); }
-  inline double Y() { return(y); }
-  inline double Z() { return(z); }
-
-  inline UnitVector vperp() {double tp, pp; pp = p;  
-     tp = t+0.5*M_PI; if (tp > M_PI) { tp = tp-M_PI; pp += M_PI*0.5; }
-     return(UnitVector(tp, pp)); }
-  inline UnitVector vperp2() {double tp, pp; tp = M_PI*0.5;  pp = p + 0.5*M_PI; 
-     return(UnitVector(tp, pp)); }
-
-  inline friend ostream& operator << (ostream& s, UnitVector const& u)
-      { s << "UnitVector: X,Y,Z= (" << u.x << "," << u.y << "," << u.z 
-          << ") Teta,Phi= " << u.t << "," << u.p << endl ;  return(s); }
-
-  double x,y,z;
-  double t,p;
-};
-
-
-#endif
Index: trunk/SophyaLib/Samba/unitvector.cc
===================================================================
--- trunk/SophyaLib/Samba/unitvector.cc	(revision 262)
+++ trunk/SophyaLib/Samba/unitvector.cc	(revision 262)
@@ -0,0 +1,59 @@
+#include <math.h>
+#include "unitvector.h"
+
+UnitVector::UnitVector() : Vector3d(1.,0.,0.) 
+{
+}
+
+UnitVector::UnitVector(double x, double y, double z) : Vector3d(x,y,z) 
+{
+  this->Normalize();
+}
+
+UnitVector::UnitVector(double theta, double phi) : Vector3d(theta,phi) 
+{
+  this->Normalize();
+}
+
+UnitVector::UnitVector(const Vector3d& v) : Vector3d(v) 
+{
+  this->Normalize();
+}
+
+Vector3d& UnitVector::operator=(const Vector3d& v)
+{
+  Setxyz(v.X(),v.Y(),v.Z());
+  this->Normalize();
+  return *this;
+}
+
+Vector3d& UnitVector::operator+=(const Vector3d& v)
+{
+  Setxyz(_x+v.X(),_y+v.Y(),_z+v.Z());
+  this->Normalize();
+  return *this;
+}
+
+Vector3d& UnitVector::operator-=(const Vector3d& v)
+{
+  Setxyz(_x-v.X(),_y-v.Y(),_z-v.Z());
+  this->Normalize();
+  return *this;
+}
+
+Vector3d UnitVector::operator+(const Vector3d& v) const
+{
+  return UnitVector(_x+v.X(),_y+v.Y(),_z+v.Z());
+}
+
+Vector3d UnitVector::operator-(const Vector3d& v) const
+{
+  return UnitVector(_x-v.X(),_y-v.Y(),_z-v.Z());
+}
+
+void UnitVector::Print(ostream& os) const 
+{
+  os << "UnitVector : (X,Y,Z)= (" << _x << ", " << _y << ", " << _z 
+     << ") Theta,Phi= " << _theta << ", " << _phi << "\n"
+     << "norme =" << this->Norm() << endl;
+}
Index: trunk/SophyaLib/Samba/unitvector.h
===================================================================
--- trunk/SophyaLib/Samba/unitvector.h	(revision 262)
+++ trunk/SophyaLib/Samba/unitvector.h	(revision 262)
@@ -0,0 +1,31 @@
+#ifndef UNITVECTOR_H_SEEN
+#define UNITVECTOR_H_SEEN
+
+#include <iostream.h>
+#include <stdio.h>
+#include <string.h>
+#include "vector3d.h"
+
+class UnitVector : public Vector3d 
+{
+  
+ public:
+  
+  UnitVector();
+  UnitVector(double x, double y, double z);
+  UnitVector(double theta, double phi=0);
+  UnitVector(const Vector3d&);
+  
+  Vector3d& operator=(const Vector3d&);
+  Vector3d& operator+=(const Vector3d&);
+  Vector3d& operator-=(const Vector3d&);
+  Vector3d operator+(const Vector3d&) const;
+  Vector3d operator-(const Vector3d&) const;
+
+  void Print(ostream& os) const;
+  
+};
+
+#endif
+
+
Index: trunk/SophyaLib/Samba/utilgeom.cc
===================================================================
--- trunk/SophyaLib/Samba/utilgeom.cc	(revision 262)
+++ trunk/SophyaLib/Samba/utilgeom.cc	(revision 262)
@@ -0,0 +1,19 @@
+#ifndef MYFUNC_H_SEEN
+#define MYFUNC_H_SEEN
+#include <math.h>
+#include <vector>
+
+const double pi2=2.*M_PI;
+const double dtor=M_PI/180.;
+const double pi_over_2=M_PI/2.;
+
+int sign(double);
+double absd(double);
+double mod(double, double=0);
+void swap(double&, double&);
+double min(double, double);
+double max(double, double);
+int arrondi(double);
+double scangle(double sinus, double cosinus);
+
+#endif
Index: trunk/SophyaLib/Samba/utilgeom.h
===================================================================
--- trunk/SophyaLib/Samba/utilgeom.h	(revision 262)
+++ trunk/SophyaLib/Samba/utilgeom.h	(revision 262)
@@ -0,0 +1,19 @@
+#ifndef MYFUNC_H_SEEN
+#define MYFUNC_H_SEEN
+#include <math.h>
+#include <vector>
+
+const double pi2=2.*M_PI;
+const double dtor=M_PI/180.;
+const double pi_over_2=M_PI/2.;
+
+int sign(double);
+double absd(double);
+double mod(double, double=0);
+void swap(double&, double&);
+double min(double, double);
+double max(double, double);
+int arrondi(double);
+double scangle(double sinus, double cosinus);
+
+#endif
Index: trunk/SophyaLib/Samba/vector3d.cc
===================================================================
--- trunk/SophyaLib/Samba/vector3d.cc	(revision 262)
+++ trunk/SophyaLib/Samba/vector3d.cc	(revision 262)
@@ -0,0 +1,301 @@
+#include <math.h>
+#include "vector3d.h"
+#include "utilgeom.h"
+
+Vector3d::Vector3d() 
+{
+  Setxyz(1.,0.,0.);
+}
+
+Vector3d::Vector3d(double x, double y, double z) 
+{
+  _x=x;
+  _y=y;
+  _z=z;
+  xyz2ThetaPhi();
+}
+
+Vector3d::Vector3d(double theta, double phi) 
+{
+  _theta=mod(theta,M_PI); // dans [0;pi]
+  _phi=mod(phi,pi2); // dans [0;2pi]
+  ThetaPhi2xyz();
+}
+
+Vector3d::Vector3d(const LongLat& ll)
+{
+  _theta=ll.Theta(); // dans [0;pi]
+  _phi=ll.Phi(); // dans [0;2pi]
+  ThetaPhi2xyz();
+}
+
+Vector3d::Vector3d(const Vector3d& v) 
+{
+  _x=v._x;
+  _y=v._y;
+  _z=v._z;
+  _theta=v._theta;
+  _phi=v._phi;
+}
+
+void Vector3d::SetThetaPhi(double theta, double phi) 
+{
+  _theta=mod(theta,M_PI);
+  _phi=mod(phi,pi2);
+  ThetaPhi2xyz();
+}
+
+void Vector3d::Setxyz(double x, double y, double z) 
+{
+  _x=x;
+  _y=y;
+  _z=z;
+  xyz2ThetaPhi();
+}
+
+void Vector3d::ThetaPhi2xyz()
+{
+  _x=sin(_theta)*cos(_phi);
+  _y=sin(_theta)*sin(_phi);
+  _z=cos(_theta);
+}
+
+void Vector3d::xyz2ThetaPhi()
+{
+  double norm=this->Norm();
+  if( norm != 0. ) 
+    {
+      _theta=acos(_z/norm); // dans [0,Pi]
+      if( mod(_theta,M_PI) == 0. ) _phi=0.; // on est sur +-Oz, le vecteur z est en phi=0
+      //      else _phi=acos(_x/sin(_theta)/norm)+M_PI*(_y<0);
+      else _phi=scangle(_y/sin(_theta)/norm,_x/sin(_theta)/norm);
+    }
+  else // vecteur nul
+    {
+      _theta=0.;
+      _phi=0.;
+    }
+}
+
+Vector3d& Vector3d::Normalize() 
+{
+  double norm=this->Norm();
+  if( norm != 0. )  (*this)/=norm;
+  else cerr << "Division par zero" << endl;
+  return *this;
+} 
+
+double Vector3d::Norm() const 
+{
+  return sqrt(_x*_x+_y*_y+_z*_z);
+}
+
+double Vector3d::Psc(const Vector3d& v) const 
+{
+  return _x*v._x+_y*v._y+_z*v._z;
+}
+
+double Vector3d::SepAngle(const Vector3d& v) const 
+{
+  double n1=this->Norm();
+  double n2=v.Norm();
+  double ret;
+  if( n1!=0. && n2!=0. ) ret=acos((this->Psc(v))/n1/n2);
+  else 
+    {
+      cerr << "Division par zero" << endl;
+      ret=0.;
+    }
+  return ret;
+}
+
+Vector3d Vector3d::Vect(const Vector3d& v) const 
+{
+  double xo=_y*v._z-_z*v._y;
+  double yo=_z*v._x-_x*v._z;
+  double zo=_x*v._y-_y*v._x;
+  return Vector3d(xo,yo,zo);
+}
+
+Vector3d Vector3d::VperpPhi() const // vecteur perpendiculaire de meme phi
+{
+  double theta;
+  if( _theta != pi_over_2 ) theta=_theta+(0.5-(_theta>pi_over_2))*M_PI; // on tourne theta de +-pi/2
+  else theta=0.;
+  return Vector3d(theta,_phi);
+}
+
+Vector3d Vector3d::VperpTheta() const // vecteur perpendiculaire de meme theta
+{
+  double phi=mod(_phi+pi_over_2,pi2); // on tourne phi de pi/2
+  return Vector3d(_theta,phi);
+}
+
+Vector3d Vector3d::EPhi() const 
+{
+  Vector3d temp;
+  if ( abs(_z) == 1. ) // si on est en +- Oz
+    {
+      temp=Vector3d(1.,0.,0.);
+    }
+  else 
+    {
+      Vector3d k(0,0,-1);
+      temp=this->Vect(k);
+      temp.Normalize();
+    }
+  return temp;
+}
+
+Vector3d Vector3d::ETheta() const 
+{
+  Vector3d temp=this->Vect(EPhi());
+  temp.Normalize();
+  return temp;
+}
+
+
+Vector3d Vector3d::Euler(double phi, double theta, double psi) const 
+{
+  double cpsi=cos(psi);
+  double ctheta=cos(theta);
+  double cphi=cos(phi);
+  double spsi=sin(psi);
+  double stheta=sin(theta);
+  double sphi=sin(phi);
+  double xnew=(cpsi*cphi-ctheta*sphi*spsi)*_x
+    +(cpsi*sphi+ctheta*cphi*spsi)*_y
+    +spsi*stheta*_z;
+  double ynew=(-spsi*cphi-ctheta*sphi*cpsi)*_x
+    +(-spsi*sphi+ctheta*cphi*cpsi)*_y
+    +cpsi*stheta*_z;
+  double znew=stheta*sphi*_x-stheta*cphi*_y+ctheta*_z;
+  return Vector3d(xnew,ynew,znew);
+}
+
+Vector3d Vector3d::InvEuler(double phi, double theta, double psi) const 
+{
+  double cpsi=cos(psi);
+  double ctheta=cos(theta);
+  double cphi=cos(phi);
+  double spsi=sin(psi);
+  double stheta=sin(theta);
+  double sphi=sin(phi);
+  double xnew=(cpsi*cphi-ctheta*sphi*spsi)*_x
+    -(spsi*cphi+ctheta*sphi*cpsi)*_y
+    +sphi*stheta*_z;
+  double ynew=(cpsi*sphi+ctheta*cphi*spsi)*_x
+    +(-spsi*sphi+ctheta*cphi*cpsi)*_y
+    -cphi*stheta*_z;
+  double znew=stheta*spsi*_x+stheta*cpsi*_y+ctheta*_z;
+  return Vector3d(xnew,ynew,znew);
+}
+
+Vector3d Vector3d::Rotate(const Vector3d& omega, double phi)
+{
+  Vector3d rotationaxis=omega;
+  rotationaxis.Normalize();
+  double n=this->Psc(rotationaxis);
+  Vector3d myself=*this;
+  Vector3d rotate=n*rotationaxis+(myself-n*rotationaxis)*cos(phi)-(myself^rotationaxis)*sin(phi);
+  return rotate;
+}   
+
+Vector3d& Vector3d::operator+=(const Vector3d& v) 
+{
+  *this=*this+v;
+  return *this;
+}
+
+Vector3d& Vector3d::operator-=(const Vector3d& v) 
+{
+  *this=*this-v;
+  return *this;
+}
+
+Vector3d& Vector3d::operator+=(double d) 
+{
+  Setxyz(_x+d,_y+d,_z+d);
+  return *this;
+}
+
+Vector3d& Vector3d::operator/=(double d) 
+{
+  if( d != 0. ) Setxyz(_x/d,_y/d,_z/d);
+  else cerr << "Division par zero." << endl;
+  return *this;
+}
+
+Vector3d& Vector3d::operator*=(double d) 
+{
+  Setxyz(_x*d,_y*d,_z*d);
+  return *this;
+}
+
+Vector3d Vector3d::operator^(const Vector3d& v) const 
+{
+  return this->Vect(v);
+}
+
+Vector3d Vector3d::operator+(double d) const 
+{
+  return Vector3d(_x+d,_y+d,_z+d);
+}
+
+Vector3d Vector3d::operator+(const Vector3d& v) const 
+{
+  return Vector3d(_x+v._x,_y+v._y,_z+v._z);
+}
+
+Vector3d Vector3d::operator-(double d) const 
+{
+  return *this+(-d);
+}
+
+Vector3d Vector3d::operator-(const Vector3d& v) const 
+{
+  return *this+(v*(-1.));
+}
+
+Vector3d Vector3d::operator*(double d) const 
+{
+  return Vector3d(d*_x,d*_y,d*_z);
+}
+
+double Vector3d::operator*(const Vector3d& v) const 
+{
+  return this->Psc(v);
+}
+
+Vector3d Vector3d::operator/(double d) const 
+{
+  Vector3d ret=*this;
+  if( d != 0. ) ret/=d;
+  else  cerr << "Division par zero." << endl;
+  return ret;
+}
+
+Vector3d& Vector3d::operator=(const Vector3d& v) 
+{
+  if( this != &v ) 
+    {
+      _x=v._x; 
+      _y=v._y; 
+      _z=v._z; 
+      _theta=v._theta; 
+      _phi=v._phi;
+    }
+  return *this;
+}
+
+bool Vector3d::operator==(const Vector3d& v) 
+{
+  return (this==&v);
+}
+
+void Vector3d::Print(ostream& os) const 
+{
+  os << "Vector3: (X,Y,Z)= (" << _x << ", " << _y << ", " << _z 
+     << ") --- Theta,Phi= " << _theta << ", " << _phi << "\n"
+     << "Norme = " << this->Norm() << endl;
+}
Index: trunk/SophyaLib/Samba/vector3d.h
===================================================================
--- trunk/SophyaLib/Samba/vector3d.h	(revision 262)
+++ trunk/SophyaLib/Samba/vector3d.h	(revision 262)
@@ -0,0 +1,119 @@
+#ifndef VECTOR3D_H_SEEN
+#define VECTOR3D_H_SEEN
+
+#include <math.h>
+#include <iostream.h>
+#include <stdio.h>
+#include <string.h>
+#include "longlat.h"
+
+/*
+  Geometrie en dimension 3. 
+  Tous les calculs sont faits en radians 
+  et en coordonnees spheriques theta,phi
+  pour les rotations (angles d'Euler) ma source est 
+  "Classical Mechanics" 2nd edition, H. Goldstein, Addison Wesley
+*/
+
+class Vector3d 
+{
+
+ public:
+  
+  Vector3d();
+  Vector3d(double x, double y, double z);
+  Vector3d(double theta, double phi);
+  Vector3d(const LongLat&);
+  Vector3d(const Vector3d&);
+
+  void Setxyz(double x, double y, double z);
+  void SetThetaPhi(double theta,  double phi);
+  void ThetaPhi2xyz();
+  void xyz2ThetaPhi();
+  double Theta() const {return _theta;}
+  double Phi() const {return _phi;}
+  double X() const {return _x;}
+  double Y() const {return _y;}
+  double Z() const {return _z;}
+
+  Vector3d& Normalize();
+  double Norm() const;
+
+  // produit scalaire
+  double Psc(const Vector3d&) const;
+
+  // ecart angulaire entre 2 vecteurs dans [0,Pi]
+  double SepAngle(const Vector3d&) const;
+
+  // produit vectoriel
+  Vector3d Vect(const Vector3d&) const;
+
+  // vecteur perpendiculaire de meme phi
+  Vector3d VperpPhi() const;
+
+  // vecteur perpendiculaire de meme theta
+  Vector3d VperpTheta() const;
+
+  Vector3d ETheta() const;
+  Vector3d EPhi() const;
+
+  // rotations d'Euler
+  Vector3d Euler(double, double, double) const;
+
+  // rotation inverse
+  Vector3d InvEuler(double, double, double) const;
+
+  // rotation d'angle phi autour d'un axe omega (regle du tire-bouchon)
+  Vector3d Rotate(const Vector3d& omega,double phi);
+
+  virtual Vector3d& operator=(const Vector3d&);
+  virtual Vector3d& operator+=(const Vector3d&);
+  virtual Vector3d& operator-=(const Vector3d&);
+  virtual Vector3d operator+(const Vector3d&) const;
+  virtual Vector3d operator-(const Vector3d&) const;
+
+  Vector3d& operator+=(double);
+  Vector3d& operator/=(double);
+  Vector3d& operator*=(double);
+
+  Vector3d operator+(double) const;
+  Vector3d operator-(double) const;
+  Vector3d operator*(double) const;
+  Vector3d operator/(double) const;
+  Vector3d operator^(const Vector3d&) const; // produit vectoriel
+  double operator*(const Vector3d&) const; // produit scalaire
+  bool operator==(const Vector3d&);
+  
+  virtual void Print(ostream& os) const;
+
+ protected:
+
+  double _x;
+  double _y;
+  double _z;
+  double _theta;
+  double _phi;
+
+};
+
+inline ostream& operator<<(ostream& s, const Vector3d& v) 
+{  
+  v.Print(s);  
+  return s;  
+}
+
+// fonctions globales
+
+inline Vector3d operator*(double d, const Vector3d& v) 
+{
+  return v*d;
+}
+
+inline Vector3d operator+(double d, const Vector3d& v) 
+{
+  return v+d;
+}
+
+#endif
+
+
