// gen_8d_sphere.c  four to eight dimensional sphere
//       area = d/dr volume         volume = integral area dr  
//              D-1 surface         D volume 
// 8D 8-sphere  1/3 Pi^4 R^7        1/24 Pi^4 R^8 
//
// one definition of sequence of n-spheres 
// a1, a2, a3, a4, a5, a6, a7 are angles, typ: phi, theta, ...
// x1, x2, x3, x4, x5, x6, x7, x8 are orthogonal coordinates, x, y, z, ... 
// x1^2 + x2^2 + x3^2 + x4^2 + x5^2 + x6^2 + x7^2 +x8^2 = R^2 
// Radius R = sqrt(R^2) 
//
// 8D 8-sphere 
// 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) cos(a5) 
// x6 = R sin(a1) sin(a2) sin(a3) sin(a4) sin(a5) cos(a6)  
// x7 = R sin(a1) sin(a2) sin(a3) sin(a4) sin(a5) sin(a6) cos(a7) 
// x8 = R sin(a1) sin(a2) sin(a3) sin(a4) sin(a5) sin(a6) sin(a7) 
//
// R = sqrt(x1^2 + x2^2 + x3^2 + x4^2 + x5^2 + x6^2 + x7^2 + x8^2))
// a1 = acos(x1/R)
// a2 = acos(x2/sqt(x2^2 + x3^2 + x4^2 + x5^2 + x6^2 + x7^2 + x8^2))
// a3 = acos(x3/sqt(x3^2 + x4^2 + x5^2 + x6^2 + x7^2 + x8^2))
// a4 = acos(x4/sqt(x4^2 + x5^2 + x6^2 + x7^2 + x8^2))
// a5 = acos(x5/sqt(x5^2 + x6^2 + x7^2 + x8^2))
// a6 = acos(x6/sqt(x6^2 + x7^2 + x8^2))
// a7 = acos(x7/sqt(x7^2 + x8^2))      if x8>=0
// a7 = 2Pi-acos(x7/sqt(x7^2 + x8^2))  if x8<0
// volume = 1/24 Pi^4 R^9  integral of area wrt R
// area   = 1/3  Pi^4 R^8  derivative of volume wrt R

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define Pi 3.1415926535897932384626433832795028841971
#undef abs
#define abs(a) ((a)<0.0?(-(a)):(a))

static double Rmin = 1.0;
static double Rmax = 1.0;
static int nR = 1;
static double a7min = 0.0;
static double a7max = 2.0*Pi;
static int na7 = 9;          //  2Pi/(na7-1)=Pi/4 steps
static double a6min = 0.0;
static double a6max = Pi;
static int na6 = 5;          // Pi/4 steps
static double a5min = 0.0;
static double a5max = Pi;
static int na5 = 5;
static double a4min = 0.0;
static double a4max = Pi;
static int na4 = 5;
static double a3min = 0.0;
static double a3max = Pi;
static int na3 = 5;
static double a2min = 0.0;
static double a2max = Pi;
static int na2 = 5;
static double a1min = 0.0;
static double a1max = Pi;
static int na1 = 5;
static char fname[] = "sphere8d.dat";
static FILE* fout;
static double pts[65536][9];

static double x8(double R, double a7, double a6, double a5,
                 double a4, double a3, double a2, double a1) 
{
  return R*sin(a1)*sin(a2)*sin(a3)*sin(a4)*sin(a5)*sin(a6)*sin(a7); 
} /* end x8 */

static double x7(double R, double a7, double a6, double a5, 
                 double a4, double a3, double a2, double a1) 
{
  return R*sin(a1)*sin(a2)*sin(a3)*sin(a4)*sin(a5)*sin(a6)*cos(a7);
} /* end x7 */

static double x6(double R, double a7, double a6, double a5, 
                 double a4, double a3, double a2, double a1) 
{
  return R*sin(a1)*sin(a2)*sin(a3)*sin(a4)*sin(a5)*cos(a6);
} /* end x6 */

static double x5(double R, double a7, double a6, double a5, 
                 double a4, double a3, double a2, double a1) 
{
  return R*sin(a1)*sin(a2)*sin(a3)*sin(a4)*cos(a5);
} /* end x5 */

static double x4(double R, double a7, double a6, double a5, 
                 double a4, double a3, double a2, double a1) 
{
  return R*sin(a1)*sin(a2)*sin(a3)*cos(a4);
} /* end x4 */

static double x3(double R, double a7, double a6, double a5, 
                 double a4, double a3, double a2, double a1) 
{
  return R*sin(a1)*sin(a2)*cos(a3);
} /* end x3 */

static double x2(double R, double a7, double a6, double a5,
                 double a4, double a3, double a2, double a1) 
{
  return R*sin(a1)*cos(a2); 
} /* end x2 */

static double x1(double R, double a7, double a6, double a5,
                 double a4, double a3, double a2, double a1) 
{
  return R*cos(a1);
} /* end x1 */


int main(int argc, char* argv[])
{
  double vR, va1, va2, va3, va4, va5, va6, va7, va8;
  double hR, ha1, ha2, ha3, ha4, ha5, ha6, ha7, ha8;
  double da1, da2, da3, da4, da5, da6, da7, da7r, da1c;
  double r, r2;
  double x1min, x2min, x3min, x4min, x5min, x6min, x7min, x8min;
  double x1max, x2max, x3max, x4max, x5max, x6max, x7max, x8max;
  double x1V, x2V, x3V, x4V, x5V, x6V, x7V, x8V;
  int iR, ia1, ia2, ia3, ia4, ia5, ia6, ia7, i;
  int npts = 0;
  int count = 0;
  int delcount = 0;
  int nerr = 0;

  printf("sphere8d.c writing sphere8d.dat \n");
  printf("a1, a2, a3, a4, a5, a6, a7 angles, r spherical radius \n");
  printf("x1, x2, x3, x4, x5, x6, x7, x8 cartesian coordinates \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) cos(a5) \n");
  printf("x6 = R sin(a1) sin(a2) sin(a3) sin(a4) sin(a5) cos(a6) \n");
  printf("x7 = R sin(a1) sin(a2) sin(a3) sin(a4) sin(a5) sin(a6) cos(a7) \n");
  printf("x8 = R sin(a1) sin(a2) sin(a3) sin(a4) sin(a5) sin(a6) sin(a7) \n");
  printf(" \n");
  printf("R = sqrt(x1^2 + x2^2 + x3^2 + x4^2 + x5^2 + x6^2 + x7^2 + x8^2) \n");
  printf("a1 = acos(x1/R) \n");
  printf("a2 = acos(x2/sqrt(x2^2 + x3^2 + x4^2 + x5^2 + x6^2 + x7^2 + x8^2))\n");
  printf("a3 = acos(x3/sqrt(x3^2 + x4^2 + x5^2 + x6^2 + x7^2 + x8^2))\n");
  printf("a4 = acos(x4/sqrt(x4^2 + x5^2 + x6^2 + x7^2 + x8^2))\n");
  printf("a5 = acos(x5/sqrt(x5^2 + x6^2 + x7^2 + x8^2))\n");
  printf("a6 = acos(x6/sqrt(x6^2 + x7^2 + x8^2))\n");
  printf("a7 = acos(x7/sqrt(x7^2 + x8^2))      if x8>=0 \n");
  printf("a7 = 2Pi-acos(x7/sqrt(x7^2 + x8^2))  if x8<0 \n");
  printf(" \n");
  printf("volume = 1/24 Pi^4 R^9 \n");
  printf("area   = 1/3  Pi^4 R^8  derivative of volume wrt R \n");
  printf(" \n");

  if(nR==1) hR = 1.0;
  else hR  = (Rmax-Rmin)/(double)(nR-1);
  
  ha1 = (a1max-a1min)/(double)(na1-1);
  ha2 = (a2max-a2min)/(double)(na2-1);
  ha3 = (a3max-a3min)/(double)(na3-1);
  ha4 = (a4max-a4min)/(double)(na4-1);
  ha5 = (a5max-a5min)/(double)(na5-1);
  ha6 = (a6max-a6min)/(double)(na6-1);
  ha7 = (a7max-a7min)/(double)(na7-1);
  printf("Rmin=%f,  Rmax=%f,  hR=%f,  nR=%d \n",  Rmin,  Rmax,  hR,  nR);
  printf("a1min=%f, a1max=%f, ha1=%f, na1=%d \n", a1min, a1max, ha1, na1);
  printf("a2min=%f, a2max=%f, ha2=%f, na2=%d \n", a2min, a2max, ha2, na2);
  printf("a3min=%f, a3max=%f, ha3=%f, na3=%d \n", a3min, a3max, ha3, na3);
  printf("a4min=%f, a4max=%f, ha4=%f, na4=%d \n", a4min, a4max, ha4, na4);
  printf("a5min=%f, a5max=%f, ha5=%f, na5=%d \n", a5min, a5max, ha5, na5);
  printf("a6min=%f, a6max=%f, ha6=%f, na6=%d \n", a6min, a6max, ha6, na6);
  printf("a7min=%f, a7max=%f, ha7=%f, na7=%d \n", a7min, a7max, ha7, na7);

  fout = fopen(fname, "w");
  fprintf(fout, "%d %d %d %d %d %d %d %d %d \n", 8, nR, na1, na2, na3,
	  na4, na5, na6, na7); /* no lead space */
  for(iR=0; iR<nR; iR++) // actually last
  {
    vR = Rmin+hR*(double)iR; 
    for(ia1=0; ia1<na1; ia1++)
    {
      va1 = a1min+ha1*(double)ia1;
      for(ia2=0; ia2<na2; ia2++)
      {
        va2 = a2min+ha2*(double)ia2;
        for(ia3=0; ia3<na3; ia3++)
        {
          va3 = a3min+ha3*(double)ia3;
          for(ia4=0; ia4<na4; ia4++)
          {
            va4 = a4min+ha4*(double)ia4;
            for(ia5=0; ia5<na5; ia5++)
            {
              va5 = a5min+ha5*(double)ia5;
              for(ia6=0; ia6<na6; ia6++)
              {
                va6 = a6min+ha6*(double)ia6;
                for(ia7=0; ia7<na7; ia7++)
                {
                  va7 = a7min+ha7*(double)ia7;

                  count++;
                  x8V = x8(vR, va7, va6, va5, va4, va3, va2, va1);
	          pts[npts][8] = x8V;
                  x7V = x7(vR, va7, va6, va5, va4, va3, va2, va1);
	          pts[npts][7] = x7V;
                  x6V = x6(vR, va7, va6, va5, va4, va3, va2, va1);
	          pts[npts][6] = x6V;
                  x5V = x5(vR, va7, va6, va5, va4, va3, va2, va1);
	          pts[npts][5] = x5V;
                  x4V = x4(vR, va7, va6, va5, va4, va3, va2, va1);
	          pts[npts][4] = x4V;
                  x3V = x3(vR, va7, va6, va5, va4, va3, va2, va1);
	          pts[npts][3] = x3V;
                  x2V = x2(vR, va7, va6, va5, va4, va3, va2, va1);
	          pts[npts][2] = x2V;
                  x1V = x1(vR, va7, va6, va5, va4, va3, va2, va1);
	          pts[npts][1] = x1V;
	          pts[npts][0] = vR;
	          npts++;

                  // some checking
                  if(ia1==2 && ia2==2 && ia3==2 && ia4==2 && ia5==2 &&
                     ia6==2) printf("ia7=%d, x8V=%f \n", ia7, x8V);

                  r2 = x1V*x1V + x2V*x2V + x3V*x3V + x4V*x4V + x5V*x5V +
                       x6V*x6V + x7V*x7V + x8V*x8V;
                  r = sqrt(r2);
		  if(abs(vR*vR-r2)>0.01) 
                     printf("r2err=%e, at a=%d,%d,%d,%d,%d,%d,%d \n",
                            vR*vR-r2, ia1, ia2, ia3, ia4, ia5, ia6, ia7);


                  da1 = acos(x1V/sqrt(r2));
                  if(abs(va1-da1)>0.01 && nerr<50)
		  {
		    printf("erra1 acos(x81, R) va1=%f, da1=%f, x1=%f \n",
			   va1, da1, x1V); 
                    nerr++;
                  }

                  da1 = x2V*x2V + x3V*x3V + x4V*x4V + x5V*x5V +
                        x6V*x6V + x7V*x7V + x8V*x8V;
                  da1 = sqrt(da1);
                  da1 = atan2(da1, x1V);
		  if(abs(da1-va1)>0.01 && nerr<50) 
		  {
                    printf("a1err a1=%f, da1=%f, x1=%f at %d,%d,%d,%d,%d,%d,%d \n",
			   va1, da1, x1V, ia1, ia2, ia3, ia4, ia5, ia6, ia7);
                    nerr++;
		  }

                  da2 = x2V*x2V + x3V*x3V + x4V*x4V + x5V*x5V +
                        x6V*x6V + x7V*x7V + x8V*x8V;
                  da2 = sqrt(da2);
                  da2 = acos(x2V/da2);
		  if(abs(da2-va2)>0.01 && nerr<50) 
		  {
                    printf("a2err a2=%f, da2=%f, x2=%f, at %d,%d,%d,%d,%d,%d,%d \n",
			   va2, da2, x2V, ia1, ia2, ia3, ia4, ia5, ia6, ia7);
                    nerr++;
		  }

                  da3 = x3V*x3V + x4V*x4V + x5V*x5V + x6V*x6V + x7V*x7V + x8V*x8V;
                  da3 = sqrt(da3);
                  da3 = acos(x3V/da3);
		  if(abs(da3-va3)>0.01 && nerr<50) 
		  {
                    printf("a3err a3=%f, da3=%f, x3=%f, at %d,%d,%d,%d,%d,%d,%d \n",
			   va3, da3, x3V, ia1, ia2, ia3, ia4, ia5, ia6, ia7);
                    nerr++;
		  }

                  da4 = x4V*x4V + x5V*x5V + x6V*x6V + x7V*x7V + x8V*x8V;
                  da4 = sqrt(da4);
                  da4 = acos(x4V/da4);
		  if(abs(da4-va4)>0.01 && nerr<50) 
		  {
                    printf("a4err a4=%f, da4=%f, x4=%f, at %d,%d,%d,%d,%d,%d,%d \n",
			   va4, da4, x4V, ia1, ia2, ia3, ia4, ia5, ia6, ia7);
                    nerr++;
		  }

                  da5 = x5V*x5V + x6V*x6V + x7V*x7V + x8V*x8V;
                  da5 = sqrt(da5);
                  da5 = acos(x5V/da5);
		  if(abs(da5-va5)>0.01 && nerr<50) 
		  {
                    printf("a5err a5=%f, da5=%f, x5=%f, at %d,%d,%d,%d,%d,%d,%d \n",
			   va5, da5, x5V, ia1, ia2, ia3, ia4, ia5, ia6, ia7);
                    nerr++;
		  }

                  da6 = x6V*x6V + x7V*x7V + x8V*x8V;
                  da6 = sqrt(da6);
                  da6 = acos(x6V/da6);
		  if(abs(da6-va6)>0.01 && nerr<50) 
		  {
                    printf("a6err a6=%f, da6=%f, x6=%f, at %d,%d,%d,%d,%d,%d,%d \n",
			   va6, da6, x6V, ia1, ia2, ia3, ia4, ia5, ia6, ia7);
                    nerr++;
		  }

                  da7 = x7V*x7V + x8V*x8V;
                  da7 = sqrt(da7);
                  da7 = acos(x7V/da7);
                  if(x8V<0.0) da7 = 2.0*Pi-da7;
		  if(abs(da7-va7)>0.01 && nerr<50) 
		  {
                    printf("a7err a7=%f, da7=%f, x7=%f, at %d,%d,%d,%d,%d,%d,%d \n",
			   va7, da7, x7V, ia1, ia2, ia3, ia4, ia5, ia6, ia7);
                    nerr++;
		  }

                  da7  = atan2(x8V, x7V);
                  if(da7<0.0) da7 = 2.0*Pi+da7;
                  if(abs(va7-da7)>0.01 && nerr<50 && abs(x8V)+abs(x7V)>0.01)
		  {
		    printf("erra7 atan2(x8V, x7V) va7=%f, da7=%f, x8=%f, x7=%f \n",
			   va7, da7, x8V, x7V); 
                    nerr++;
                  }

                  // eliminate duplicates
	          for(i=0; i<npts-1; i++)
	          {
                    if(abs(pts[npts-1][0] - pts[i][0])<0.01 &&
	               abs(pts[npts-1][1] - pts[i][1])<0.01 &&
	               abs(pts[npts-1][2] - pts[i][2])<0.01 &&
	               abs(pts[npts-1][3] - pts[i][3])<0.01 &&
	               abs(pts[npts-1][4] - pts[i][4])<0.01 &&
	               abs(pts[npts-1][5] - pts[i][5])<0.01 &&
	               abs(pts[npts-1][6] - pts[i][6])<0.01 &&
	               abs(pts[npts-1][7] - pts[i][7])<0.01 &&
                       abs(pts[npts-1][8] - pts[i][8])<0.01 )
                    {
                      npts--;
                      delcount--;
                      break;
	            } // end if
                  } // end i
	        } // end va7
              } // end va6
            } // end va5
          } // end va4
        } // end va3
      } // end va2
    } // end va1
  } // end vR
  for(i=0; i<npts; i++)
  {
    fprintf(fout, "%f %f %f %f %f %f %f %f %f \n",
            pts[i][0], pts[i][1], pts[i][2], pts[i][3], pts[i][4],
            pts[i][5], pts[i][6], pts[i][7], pts[i][8]);
  }
  fclose(fout);
  printf("sphere8d.dat written \n");
  printf("statistics \n");
  x8min = pts[0][8];
  x8max = x8min;
  x7min = pts[0][7];
  x7max = x7min;
  x6min = pts[0][6];
  x6max = x6min;
  x5min = pts[0][5];
  x5max = x5min;
  x4min = pts[0][4];
  x4max = x4min;
  x3min = pts[0][3];
  x3max = x3min;
  x2min = pts[0][2];
  x2max = x2min;
  x1min = pts[0][1];
  x1max = x1min;
  for(i=1; i<npts; i++)
  {
    if(pts[i][8]<x8min) x8min = pts[i][8];
    if(pts[i][8]>x8max) x8max = pts[i][8];
    if(pts[i][7]<x7min) x7min = pts[i][7];
    if(pts[i][7]>x7max) x7max = pts[i][7];
    if(pts[i][6]<x6min) x6min = pts[i][6];
    if(pts[i][6]>x6max) x6max = pts[i][6];
    if(pts[i][5]<x5min) x5min = pts[i][5];
    if(pts[i][5]>x5max) x5max = pts[i][5];
    if(pts[i][4]<x4min) x4min = pts[i][4];
    if(pts[i][4]>x4max) x4max = pts[i][4];
    if(pts[i][3]<x3min) x3min = pts[i][3];
    if(pts[i][3]>x3max) x3max = pts[i][3];
    if(pts[i][2]<x2min) x2min = pts[i][2];
    if(pts[i][2]>x2max) x2max = pts[i][2];
    if(pts[i][1]<x1min) x1min = pts[i][1];
    if(pts[i][1]>x1max) x1max = pts[i][1];
  }
  printf("x1min=%f, x1max=%f \n", x1min, x1max);
  printf("x2min=%f, x2max=%f \n", x2min, x2max);
  printf("x3min=%f, x3max=%f \n", x3min, x3max);
  printf("x4min=%f, x4max=%f \n", x4min, x4max);
  printf("x5min=%f, x5max=%f \n", x5min, x5max);
  printf("x6min=%f, x6max=%f \n", x6min, x6max);
  printf("x7min=%f, x7max=%f \n", x7min, x7max);
  printf("x8min=%f, x8max=%f \n", x8min, x8max);
  printf("count=%d, delcount=%d, npts=%d \n", count, delcount, npts);
  return 0;
} /* end main of gen_8d_sphere.c */