/* test_faces.c  see faces.out */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#undef abs
#define abs(a) ((a)<0.0?(-(a)):(a))
#undef max
#define max(a,b) ((a)>(b)?(a):(b))

int main(int argc, char *argv[])
{
  double a1, a2, a3, a4, a5; /* angles */
  double a1c, a2c, a3c, a4c, a5c; /* angles recomputed */
  double x1, x2, x3, x4, x5, x6; /* orthoganal coordinates */
  double R = 1.0;
  double sumsq, sq2, sq3, sq4, err, maxerr;
  double Pi = 3.141592653589793238462643;
  double twoPi = 2.0*Pi;
  double halfPi = Pi/2.0;
  double fourthPi = Pi/4.0;
  double eighthPi = Pi/8.0;
  int i1, i2, i3, i4, i5;
  FILE *outp;
  char filename[] = "sphere4d.dat";

  printf("test_faces.c running, R = 1 \n");
  printf("2D circle  1 angle radius 1 \n");
  printf("  x1 = R cos(a1) \n");
  printf("  x2 = R sin(a1) \n");
  maxerr = 0.0;
  a1 = 0.0;
  for(i1=0; i1<8; i1++)
  {
    x1 = R*cos(a1);
    x2 = R*sin(a1);
    sumsq = x1*x1 + x2*x2;
    a1c = acos(x1/sqrt(sumsq));
    if(x2<0.0) a1c = twoPi - a1c;
    maxerr = max(maxerr,abs(R*R-sumsq));
    printf("a1=%f, x1=%f, x2=%f, a1c=%f \n", a1, x1, x2, a1c);
    a1 += fourthPi;
  }
  printf("                              maxerr=%e \n", maxerr);

  printf("3D sphere  2 angles radius 1 \n");
  printf("    x1 = R cos(a1) \n");
  printf("    x2 = R sin(a1)*cos(a2) \n");
  printf("    x3 = R sin(a1)*sin(a2) \n");
  maxerr = 0.0;
  a1 = 0.0;
  for(i1=0; i1<=4; i1++) /* both include poles */
  {
    a2 = 0.0;
    for(i2=0; i2<8; i2++)
    {
      x1 = R*cos(a1);
      x2 = R*sin(a1)*cos(a2);
      x3 = R*sin(a1)*sin(a2);
      sumsq = x1*x1 + x2*x2 + x3*x3;
      sq2 = x2*x2+x3*x3;
      a1c = acos(x1/sqrt(sumsq));
      a2c = acos(x2/sqrt(sq2));
      if(sq2<1.0e-8) a2c = 0.0;
      if(x3<0.0) a2c = twoPi - a2c;
      maxerr = max(maxerr,abs(R*R-sumsq));
      printf("a1 =%f, a2 =%f x1=%f, x2=%f, x3=%f \n", a1, a2, x1, x2, x3);
      printf("a1c=%f, a2c=%f \n", a1c, a2c);
      a2 += fourthPi;
      if(a1==0.0 || abs(a1-Pi)<0.0001) break; /* only one point on pole */
    }
    a1 += fourthPi;
  }
  printf("                              maxerr=%e \n", maxerr);

  printf("4D sphere  3 angels radius 1 \n");
  printf("      x1 = R cos(a1) \n");
  printf("      x2 = R sin(a1)*cos(a2) \n");
  printf("      x3 = R sin(a1)*sin(a2)*cos(a3) \n");
  printf("      x4 = R sin(a1)*sin(a2)*sin(a3) \n");

  maxerr = 0.0;
  a1 = 0.0;
  for(i1=0; i1<4; i1++)
  {
    a2 = 0.0;
    for(i2=0; i2<4; i2++)
    {
      a3 = 0.0;
      for(i3=0; i3<8; i3++)
      {
        x1 = R*cos(a1);
        x2 = R*sin(a1)*cos(a2);
        x3 = R*sin(a1)*sin(a2)*cos(a3);
        x4 = R*sin(a1)*sin(a2)*sin(a3);
        sumsq = x1*x1 + x2*x2 + x3*x3 + x4*x4;
        sq3 = x2*x2+x3*x3+x4*x4;
        sq2 = x3*x3+x4*x4;
        a1c = acos(x1/sqrt(sumsq));
        a2c = acos(x2/sqrt(sq3));
        if(sq3<1.0e-8) a2c = 0.0;
        a3c = acos(x3/sqrt(sq2));
        if(sq2<1.0e-8) a3c = 0.0;
        if(x4<0.0) a3c = twoPi - a3c;

        maxerr = max(maxerr,abs(R*R-sumsq));
        printf("a1 =%f, a2 =%f a3 =%f \nx1=%f, x2=%f, x3=%f x4=%f\n",
               a1, a2, a3, x1, x2, x3, x4);
        printf("a1c=%f, a2c=%f a3c=%f \n\n",
               a1c, a2c, a3c);
        a3 += fourthPi;
        if(a1==0.0 || abs(a1-Pi)<0.0001) break; /* only one point on pole */
        if(a2==0.0 || abs(a2-Pi)<0.0001) break; /* only one point on pole */
      }
      a2 += fourthPi;
      if(a2==0.0 || abs(a2-Pi)<0.0001) break; /* only one point on pole */
    }
    a1 += fourthPi;
  }
  printf("                              maxerr=%e \n", maxerr);


  printf("5D sphere  4 angles radius 1 \n");
  printf("    x1 = R cos(a1) \n");
  printf("    x2 = R sin(a1) cos(a2) \n");
  printf("    x3 = R sin(a1) sin(a2) cos(a3) \n");
  printf("    x4 = R sin(a1) sin(a2) sin(a3) cos(a4) \n");
  printf("    x5 = R sin(a1) sin(a2) sin(a3) sin(a4) \n");
  maxerr = 0.0;
  a1 = 0.0;
  for(i1=0; i1<=4; i1++)
  {
    a2 = 0.0;
    for(i2=0; i2<=4; i2++)
    {
      a3 = 0.0;
      for(i3=0; i3<=4; i3++) /* both include poles */
      {
        a4 = 0.0;
        for(i4=0; i4<8; i4++)
        {
          x1 = R*cos(a1);
          x2 = R*sin(a1)*cos(a2);
          x3 = R*sin(a1)*sin(a2)*cos(a3);
          x4 = R*sin(a1)*sin(a2)*sin(a3)*cos(a4);
          x5 = R*sin(a1)*sin(a2)*sin(a3)*sin(a4);
          sumsq = x1*x1 + x2*x2 + x3*x3 + x4*x4 + x5*x5;
          sq4 = x2*x2 + x3*x3 + x4*x4 + x5*x5;
          sq3 = x3*x3 + x4*x4 + x5*x5;
          sq2 = x4*x4 + x5*x5;
          a1c = acos(x1/sqrt(sumsq));
          a2c = acos(x2/sqrt(sq4));
          if(sq4<1.0e-8) a2c = 0.0;
          a3c = acos(x3/sqrt(sq3));
          if(sq3<1.0e-8) a3c = 0.0;
          a4c = acos(x4/sqrt(sq2));
          if(sq2<1.0e-8) a4c = 0.0;
          if(x5<0.0) a4c = twoPi - a4c;

          maxerr = max(maxerr,abs(R*R-sumsq));
          printf("a1 =%f, a2 =%f a3 =%f a4 =%f\nx1=%f, x2=%f, x3=%f x4=%f x5=%f\n",
                 a1, a2, a3, a4, x1, x2, x3, x4, x5);
          printf("a1c=%f, a2c=%f a3c=%f a4c=%f \n\n",
                 a1c, a2c, a3c, a4c);
          a4 += fourthPi;
          if(a1==0.0 || abs(a1-Pi)<0.0001) break; /* only one point on pole */
          if(a2==0.0 || abs(a2-Pi)<0.0001) break; /* only one point on pole */
          if(a3==0.0 || abs(a3-Pi)<0.0001) break; /* only one point on pole */
        }
        a3 += fourthPi;
        if(a1==0.0 || abs(a1-Pi)<0.0001) break; /* only one point on pole */
        if(a2==0.0 || abs(a2-Pi)<0.0001) break; /* only one point on pole */
      }
      a2 += fourthPi;
      if(a1==0.0 || abs(a1-Pi)<0.0001) break; /* only one point on pole */
    }
    a1 += fourthPi;
  }
  printf("                              maxerr=%e \n", maxerr);

  printf("test_faces.c finished \n");
  return 0;
} /* end test_ faces.c */