*** empty log message ***

Geo/gmshSurface.cpp

-// $Id: gmshSurface.cpp,v 1.7 2007-03-02 14:36:38 geuzaine Exp $
+// $Id: gmshSurface.cpp,v 1.8 2007-03-14 10:21:12 remacle Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -70,22 +70,14 @@ gmshSurface * gmshPolarSphere::NewPolarSphere(int iSphere, double x, double y, d
   allGmshSurfaces[iSphere] = sph;
   return sph;
 }
-
 SPoint3  gmshPolarSphere::point (double parA, double parB) const
 {
-	/*double ra,phi;
-	ra=hypot(parA,parB);
-	phi=2*atan((parB/ra)/(1+parA/ra));
-	double par1=-phi;	
-	double par2=M_PI-ra;
-
-  //par2 += M_PI*.5;
-  const double x = xc + r * sin(par2) * cos(par1);
-  const double y = yc + r * sin(par2) * sin(par1);
-  const double z = zc - r * cos(par2); 
-  // printf("%g %g - %g %g %g\n",par1,par2,x,y,z);*/
-	double f=2*r/(parA*parA+parB*parB+4*r*r);
-  return SPoint3(f*2*parA*r, f*2*parB*r, f*(parA*parA+parB*parB));
+	//stereographic projection from the south pole, origin of the axis at the center of the sphere
+	//parA=2rx/(r+z)
+	//parB=2ry/(r+z)
+	double rp2=parA*parA+parB*parB;
+	double z=r*(4*r*r-rp2)/(4*r*r+rp2);
+  return SPoint3((r+z)*parA/(2*r),(r+z)*parB/(2*r),z);
 }