/* toroidal_coord.c  toroidal coordinate system      */
/*  sigma, theta, phi  mapped to x,y,z               */
/*  denom = cosh(theta)-cos(sigma)                   */
/*  x = a * sinh(theta) * cos(phi) / denom           */
/*  y = a * sinh(theta) * sin(phi) / denom           */
/*  z = a * sin(sigma)             / denom           */
/*  -Pi < sigma < Pi   theta > 0    0 < phi < 2Pi    */
/*                                                   */
/*  x,y,z mapped to sigma, theta, phi                */
/*  phi = arctan(y/x)                                */
/*       need temporaries r1 = sqrt(x^2 + y^2)       */
/*                        d1 = sqrt((r1+a)^2 + z^2)  */
/*                        d2 = sqrt((r1-a)^2 + z^2)  */
/* theta = ln(d1/d2)                                 */
/* sigma = arccos((d1^2+d2^2-4*a^2)/(2*d1*d2))       */
/*                                                   */
/* incremental volume                                */
/* dV = a^3 * sinh(theta)/denom^3 dsigma dtheta dphi */
/* Laplacian nabla^2 Psi = many terms                */
/* solutions involve e^i*terms ...                   */
/* fyi sinh(x) = (e^x - e^-x)/2                      */
/*     cosh(x) = (e^x + e^-x)/2                      */
/* thus, try a=1, theta=1, sigma and phi over range  */

#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 a, double sigma, double theta, double phi)
{
  double denom = cosh(theta)-cos(sigma);
  return a * sinh(theta) * cos(phi) / denom;
}

static double Yc(double a, double sigma, double theta, double phi)
{
  double denom = cosh(theta)-cos(sigma);
  return a * sinh(theta) * sin(phi) / denom;
}

static double Zc(double a, double sigma, double theta, double phi)
{
  double denom = cosh(theta)-cos(sigma);
  return a * sin(sigma) / denom;    
}

static double sigmac(double a, double x, double y, double z)
{
  double r1 = sqrt(x*x + y*y);
  double d1 = sqrt((r1+a)*(r1+a) + z*z);
  double d2 = sqrt((r1-a)*(r1-a) + z*z);
  double ang = acos((d1*d1 + d2*d2 - 4.0*a*a)/(2.0*d1*d2));
  if(z<0.0) return -abs(ang);
  else      return  abs(ang);
}

static double thetac(double a, double x, double y, double z)
{
  double r1 = sqrt(x*x + y*y);
  double d1 = sqrt((r1+a)*(r1+a) + z*z);
  double d2 = sqrt((r1-a)*(r1-a) + z*z);
  return log(d1/d2);
}

static double phic(double a, double x, double y, double z)
{
  double ang = atan2(y,x);
  if(ang<0.0) ang = 2.0*PI+ang; /* 0 to 2PI */ 
  return ang; 
}

static void pr(double a, double sigma, double theta, double phi)
{
  double x, y, z;
  x = Xc(a, sigma, theta, phi);
  y = Yc(a, sigma, theta, phi);
  z = Zc(a, sigma, theta, phi);
  printf("theta=%4.2f, sigma=%5.2f, phi=%5.2f, ", theta, sigma, phi);
  printf("Xc=%7.3f, Yc=%7.3f, Zc=%7.3f \n", x, y, z);
}

static void check(double a, double sigma, double theta, double phi)
{
  double x, y, z;
  double err = 0.0;
  double errs;
  double errt;
  double errp;

  x = Xc(a, sigma, theta, phi);
  y = Yc(a, sigma, theta, phi);
  z = Zc(a, sigma, theta, phi);
  errs = sigma - sigmac(a, x, y, z);
  errt = theta - thetac(a, x, y, z);
  errp = phi   -   phic(a, x, y, z);

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

int main(int argc, char * argv[])
{
  FILE * outp;
  double a, sigma, theta, phi;
  double x, y, z;
  int i, j;

  printf("toroidal_coord.c running \n");
  a = 1.0;
  sigma = 0.0;
  theta = 1.0;
  phi = 0.0;
  
  printf("a=%4.2f \n", a);

  printf("\nchanging sigma\n");
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  sigma = 0.2;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  sigma = 0.5;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  sigma = PI/2.0;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  sigma = PI;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  sigma = -PI/2.0;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  sigma = -0.5;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  sigma = -0.2;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);

  printf("\nchanging phi\n");
  sigma = 0.0;
  phi = 0.0;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  phi = 0.2;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  phi = 0.5;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  phi = PI/2.0;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  phi = PI;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  phi = 3.0*PI/2.0;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  phi = 2.0*PI-0.5;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  phi = 2.0*PI-0.2;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);

  printf("\nchanging theta\n");
  sigma = 0.0;
  phi = 0.0;
  theta = 4.0;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  theta = 3.0;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  theta = 2.0;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  theta = 1.0;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  theta = 0.9;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  theta = 0.7;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  theta = 0.5;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);

  printf("\nchanging theta and phi, sigma just changes z\n");
  sigma = 0.0;
  phi = 0.1;
  theta = 1.1;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  phi = 0.2;
  theta = 1.2;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);
  phi = 0.3;
  theta = 1.3;
  pr(a, sigma, theta, phi);
  check(a, sigma, theta, phi);

  outp = fopen("toro.dat","w");
  a = 1.0;
  theta = 1.0;
  sigma = -PI;  
  for(i=0; i<33; i++)
  {
    phi = 0.0;
    for(j=0; j<17; j++)
    {
      x = Xc(a, sigma, theta, phi);
      y = Yc(a, sigma, theta, phi);
      z = Zc(a, sigma, theta, phi);
      fprintf(outp, "%9.5f  %9.5f,  %9.5f \n", x, y, z);
      phi = phi + PI/8.0;
    }
    fprintf(outp,"\n");
    sigma = sigma + PI/16.0;
  }
  fclose(outp);
  printf("toro.dat written for gnuplot \n");
  printf("toroidal_coord.c finished \n");
  return 0;
}