Fix precision problem with asin()

Common/Numeric.cpp

-// $Id: Numeric.cpp,v 1.3 2001-12-03 08:41:43 geuzaine Exp $
+// $Id: Numeric.cpp,v 1.4 2002-01-24 17:48:28 geuzaine Exp $
 
 #include "Gmsh.h"
 #include "Numeric.h"
@@ -20,6 +20,12 @@ double myacos (double a){
   return acos (a);
 }
 
+double myasin(double a){
+  if(a<=-1.) return -Pi/2.;
+  else if(a>=1.) return Pi/2.;
+  else return asin(a);
+}
+
 void prodve (double a[3], double b[3], double c[3]){
   c[2] = a[0] * b[1] - a[1] * b[0];
   c[1] = -a[0] * b[2] + a[2] * b[0];

Common/Numeric.h

@@ -36,6 +36,7 @@
 
 double myatan2 (double a, double b);
 double myacos (double a);
+double myasin (double a);
 void prodve (double a[3], double b[3], double c[3]);
 void prosca (double a[3], double b[3], double *c);
 void norme (double a[3]);

Mesh/Create.cpp

-// $Id: Create.cpp,v 1.33 2002-01-23 16:28:00 geuzaine Exp $
+// $Id: Create.cpp,v 1.34 2002-01-24 17:48:28 geuzaine Exp $
 
 #include "Gmsh.h"
 #include "Numeric.h"
@@ -318,17 +318,20 @@ void End_Curve (Curve * c){
 	A1 = A3 = 0.;
 	f1 = f2 = R ;
       }
-      else{
+      else{ 
 	f1 = sqrt(1./sol[0]);
 	f2 = sqrt(1./sol[1]);
+	// myasin() permet de contourner les problemes de precision
+	// sur y1/f2 ou y3/f2, qui peuvent legerement etre hors de
+	// [-1,1]
 	if(x1 < 0) 
-	  A1 = -asin(y1/f2) + A4 + Pi; 
+	  A1 = -myasin(y1/f2) + A4 + Pi; 
 	else
-	  A1 = asin(y1/f2) + A4; 
+	  A1 = myasin(y1/f2) + A4; 
 	if(x3 < 0) 
-	  A3 = -asin(y3/f2) + A4 + Pi; 
+	  A3 = -myasin(y3/f2) + A4 + Pi; 
 	else
-	  A3 = asin(y3/f2) + A4; 
+	  A3 = myasin(y3/f2) + A4; 
       }
     }
     else{