// toro_cube.c  toroidal coordinate system 6 toroides on faces of cube
//  r1, r2, theta, phi  mapped to x,y,z    
//  x = (r1+r2*sin(phi))*cos(theta)
//  y = (r1+r2*sin(phi))*sin(theta)
//  z = r2 * cos(phi)
//  lenctr = r1*2.0*Pi
//  area =   r1*2.0*Pi * r2*2.0*Pi
//  volume = r1*2.0*Pi * r2*r2*Pi
//  0 < theta < 2Pi    -Pi < phi < Pi
//
//        R1 = 3/4   R2 = 1/4
//
// side        center  offset      rotation around axis 
//             xc      yc     zc | rotx  roty  rotz    |
// top    |     0       0  +1.25 |    0     0     0    |
// bottom |     0       0  -1.25 |    0     0     0    |
// right  | +1.25       0      0 |    0   +90     0    |
// left   | -1.25       0      0 |    0   -90     0    |
// front  |     0   +1.25      0 |  +90     0     0    |
// back   |     0   -1.25      0 |  -90     0     0    |
//
// outside cube  xmin=-1, xmax=+1, ymin=-1, ymax=+1, zmin=-1, zmax=+1
//   _______
//  /  T   /|       Z   -Y
// /______/ |       |  / 
// |      |R|       | /
// |  F   | /       |_____ X
// |______|/
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define PI 3.141592653589793238462643
#define nphi 16
#define ntheta 16
#undef abs
#define abs(a) ((a)<0.0?(-(a)):(a))

static double xc, yc, zc, rhox, rhoy, rhoz; // offset and rotation

static double Xpos(double r1, double r2, double theta, double phi)
{
  return (r1 + r2*sin(phi))*cos(theta);
} // end Xpos

static double Ypos(double r1, double r2, double theta, double phi)
{
  return (r1 + r2*sin(phi))*sin(theta);
} // end Ypos

static double Zpos(double r1, double r2, double theta, double phi)
{
  return r2*cos(phi);
} // end Zpos

static void matvec(double A[4][4], double X[], double Y[])
{
  int i, j;
  double val = 0.0;
  for(i=0; i<4; i++)
  {
    val = 0.0;
    for(j=0; j<4; j++)
    {
      val += A[i][j]*X[j];
    } // j
    Y[i] = val;
  } // i
} // end matvec 

void generate_toro(double r1, double r2, FILE * outp)
                   // parameters static  xc, yc, zc, rhox, rhoy, rhoz
{
  double theta, phi, x, y, z;
  double xr, yr, zr;
  double rmat[4][4] = {{1.0,0.0,0.0,0.0},
                       {0.0,1.0,0.0,0.0},
                       {0.0,0.0,1.0,0.0},
                       {0.0,0.0,0.0,1.0}};
  double rvec[4], xyz[4], ang;
  int i, j, ii, jj, nn;
  int p1, p2, p3, p4;

  printf("xc=%f, yc=%f, zc=%f \n", xc, yc, zc);
  printf("rhox=%f, rhoy=%f, rhoz=%f \n", rhox, rhoy, rhoz);

  // write file for gnuplot
  for(i=0; i<=nphi; i++)
  {
    phi = 0.0;
    for(j=0; j<=ntheta; j++)
    {
      if(rhox!=0.0)
      {
        xr = Xpos(r1, r2, theta, phi);
        yr = Ypos(r1, r2, theta, phi);
        zr = Zpos(r1, r2, theta, phi);
        x = xc + xr;
        y = yc + zr;
        z = yr;
      } // end rhox
      else if(rhoy!=0.0)
      {
        xr = Xpos(r1, r2, theta, phi);
        yr = Ypos(r1, r2, theta, phi);
        zr = Zpos(r1, r2, theta, phi);
        x = xc + zr;
        y = yc + yr;
        z = xr;
      } // end rhoy
      else
      {
        x = xc + Xpos(r1, r2, theta, phi);
        y = yc + Ypos(r1, r2, theta, phi);
        z = zc + Zpos(r1, r2, theta, phi);
      }

      fprintf(outp, "%9.5f  %9.5f  %9.5f \n", x, y, z);
      phi = phi + PI/8.0;
    }
    fprintf(outp,"\n");
    theta = theta + PI/8.0;
  }
  fprintf(outp,"\n");
} // end generate_cube

static void check(double r1, double r2, double theta, double phi)
{
  double x, y, z;
  double err = 0.0;
  double errr;
  double errt;
  double errp;

  x = Xpos(r1, r2, theta, phi);
  y = Ypos(r1, r2, theta, phi);
  z = Zpos(r1, r2, theta, phi);
  // inverse and check

  err = abs(errr) + abs(errt) + abs(errp);  
  if(err > 1.0e-5)
  {
    printf("inverse transform error = %g \n",err);
    printf("errr=%g, errt=%g, errp=%g \n", errr, errt, errp);
  }
} // end check


int main(int argc, char * argv[])
{
  FILE * outp;
  double r1, r2, sigma, theta, phi;
  double x, y, z;
  double lenctr, area, volume;
  int i, j, ii, jj, p1, p2, p3, p4;
  int nn;

  printf("toro_cube.c running \n");
 
  outp = fopen("toro_cube.dat","w");
  r1 = 0.75; // 3/4
  r2 = 0.25; // 1/4
  theta = 0.0;
  phi = -PI;  

  printf("r1=%f, r2=%f \n", r1, r2);
  printf("0 < theta < 2Pi,  -Pi < phi < Pi \n");

  lenctr = r1*2.0*PI;
  area=    r1*2.0*PI * r2*2.0*PI;
  volume = r1*2.0*PI * r2*PI*PI;

  printf("lenctr=%f, area=%f, volume=%f \n", lenctr, area, volume);

  xc = 0.0; yc = 0.0; zc =+1.25;
  rhox = 0.0; rhoy = 0.0; rhoz = 0.0;
  generate_toro(r1, r2, outp);

  xc = 0.0; yc = 0.0; zc =-1.25;
  rhox = 0.0; rhoy = 0.0; rhoz = 0.0;
  generate_toro(r1, r2, outp);

  xc = 0.0; yc = 1.25; zc = 0.0;
  rhox = 90.0; rhoy = 0.0; rhoz = 0.0;
  generate_toro(r1, r2, outp);

  xc = 0.0; yc = -1.25; zc = 0.0;
  rhox = 90.0; rhoy = 0.0; rhoz = 0.0;
  generate_toro(r1, r2, outp);

  xc = +1.25; yc = 0.0; zc = 0.0;
  rhox = 0.0; rhoy = 90.0; rhoz = 0.0;
  generate_toro(r1, r2, outp);

  xc = -1.25; yc = 0.0; zc = 0.0;
  rhox = 0.0; rhoy = 90.0; rhoz = 0.0;
  generate_toro(r1, r2, outp);


  fclose(outp);
  printf("toro_cube.dat written for gnuplot \n");

  printf("toro_cube.c finished \n");
  return 0;
} // end toro_cube.c