===================================================================
RCS file: /home/cvs/OpenXM/src/hgm/so3/src/so3_nc.c,v
retrieving revision 1.1
retrieving revision 1.3
diff -u -p -r1.1 -r1.3
--- OpenXM/src/hgm/so3/src/so3_nc.c	2013/02/07 07:38:23	1.1
+++ OpenXM/src/hgm/so3/src/so3_nc.c	2013/02/12 23:51:25	1.3
@@ -17,8 +17,10 @@ void so3_evalByS(int deg,double a,double b,double c,do
 
 #define MDEG 20
 
-int SO3_Quiet = 0;
-int SO3_Deg = 10;
+#define SO3_QUIET_DEFAULT 0
+#define SO3_Deg_DEFAULT 10
+int SO3_Quiet = SO3_QUIET_DEFAULT;
+int SO3_Deg = SO3_Deg_DEFAULT;
 #ifdef STANDALONE
 main(int argc,char *argv[]) {
   double a[3];
@@ -53,6 +55,8 @@ void so3_main(double *in1,double *in2,double *in3,doub
   double y[4];
   double t0;
   int i,j;
+  SO3_Quiet = SO3_QUIET_DEFAULT;
+  SO3_Deg = SO3_Deg_DEFAULT;
   if (*quiet) SO3_Quiet = 1;
   if (*deg) {
     SO3_Deg = *deg;
@@ -69,7 +73,7 @@ void so3_main(double *in1,double *in2,double *in3,doub
     *out = 0.0; return;
   }
   a[0] = *in1; a[1] = *in2; a[2] = *in3;
-  fprintf(stderr,"DEBUG: %lf,%lf,%lf,%lf\n",t0,a[0],a[1],a[2]);
+  //  fprintf(stderr,"DEBUG: %lf,%lf,%lf,%lf\n",t0,a[0],a[1],a[2]);
   so3_nc(a,t0,y);
 // printf("%lg, %lg, %lg, %lg\n",y[0],y[1],y[2],y[3]);
   if (!SO3_Quiet) printf("%.15e, %.15e, %.15e, %.15e\n",y[0],y[1],y[2],y[3]);
@@ -96,8 +100,17 @@ void so3_nc(double a[3],double t0,double y[4]) {
   double r;
   double y0[4];
   double myerr,myerr2;
+  double aa;
   deg = SO3_Deg;
-  for (i=0; i<4; i++) SO3_A[i]=a[i];
+  for (i=0; i<3; i++) SO3_A[i]=a[i];
+
+  /* When the argument is small, eval it only by series */
+  aa = 0.0; for (i=0; i<3; i++) aa += a[i]*a[i];
+  if (aa < 0.01) {
+    so3_evalByS(deg,a[0],a[1],a[2], 1.0, y);
+    return;
+  }
+
   SO3_R = 0.0;
   r = a[0]-a[1]-a[2]; if (r > SO3_R) SO3_R = r;
   r = -a[0]+a[1]-a[2]; if (r > SO3_R) SO3_R = r;