// toroidal_coord.c  toroidal coordinate system      
//  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                  
//                                                    
//  r2, x, y, z mapped to r1, theta, phi              
//  theta = arctan(y/x)                               
//  phi = arccos(z/r2)                                
//  r1 = x/cos(theta) - r2*sin(phi)   or              
//  r1 = y/sin(theta) - r2*sin(phi) no divide by zero 
//                                                    
//  r1, x, y, z mapped to r2, theta, phi              
//  theta = arctan(y/x)                               
//  x1 = r1*cos(theta)                                
//  y1 = r1*sin(theta)
//  phi = arctan(sqrt((x-x1)^2+(y-y1)^2)/z)
//  r2 = sqrt((x-x1)^2+(y-y1)^2 + z^2)



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

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

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

static double Zc(double r1, double r2, double theta, double phi)
{
  return r2 * cos(phi);
}


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 = Xc(r1, r2, theta, phi);
  y = Yc(r1, r2, theta, phi);
  z = Zc(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);
  }
}

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;

  printf("toroidal1_coord.c running \n");
 
  outp = fopen("toro1.dat","w");
  r1 = 1.5;
  r2 = 1.0;
  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);

  for(i=0; i<33; i++)
  {
    phi = 0.0;
    for(j=0; j<17; j++)
    {
      x = Xc(r1, r2, theta, phi);
      y = Yc(r1, r2, theta, phi);
      z = Zc(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;
  }
  fclose(outp);
  printf("toro1.dat written for gnuplot \n");
  printf("toroidal1_coord.c finished \n");
  return 0;
}