/* coordinate.c  sample coordinate conversion code, tailor to your needs */
/* cartesian coordinates are x,y  x,y,z   x,y,z,t                        */
/*                                                                       */
/* polar  radius r,  angle from x axis theta                             */  
/*        x = r * cos(theta)               r = sqrt(x*x+y*y)             */
/*        y = r * sin(theta)               theta = arctan(y/x)  or atan2 */
/*                                                                       */
/*        theta in 0 to 2Pi               need theta over complete range */
/*        r     positive                  any real x,y                   */
/*                                                                       */
/*        area = Pi*r*r                                                  */
/*                                                                       */
/* cylindrical   radius r, angle theta and height z                      */
/*        x = r * cos(theta)               r = sqrt(x*x+y*y)             */
/*        y = r * sin(theta)               theta = arctan(y/x)  or atan2 */
/*        z = z                            z = z                         */
/*                                                                       */
/*        theta in 0 to 2Pi               need theta over complete range */
/*        r     positive                  any real x,y                   */
/*                                                                       */
/*        area = 2*Pi*r*z + 2*Pi*r*r                                     */
/*        volume = Pi*r*r*z                                              */
/*                                                                       */
/* spherical  radius r, angle theta in x,y plane from positive x axis    */
/*            angle phi from positive z axis                             */
/*        x = r * sin(phi) * cos(theta)    r = sqrt(x*x+y*y+z*z)         */
/*        y = r * sin(phi) * sin(theta)    theta = arctan(y/x)           */
/*        z = r * cos(phi)                 phi = arctan(sqrt(x*x+y*y)/z) */
/*                                                                       */
/*        phi   in 0 to Pi                need phi over complete range   */
/*        theta in 0 to 2Pi               need theta over complete range */
/*        r     positive                  any real x,y,z                 */
/*                                                                       */
/*        area = 4*Pi*r*r                                                */
/*        volume = 4/3 Pi*r*r*r                                          */
/*                                                                       */
/* 4D spherical  radius r, angle theta in x,y plane from positive x axis */
/*            angle phi from positive z axis                             */
/*            angle psi from positive t axis                             */
/*        x = r * sin(psi) * sin(phi) * cos(theta)                       */ 
/*                                     r = sqrt(x*x+y*y+z*z+t*t)         */
/*        y = r * sin(psi) * sin(phi) * sin(theta)                       */
/*                                     theta = arctan(y/x)               */
/*        z = r * sin(psi) * cos(phi)                                    */
/*                                     phi = arctan(sqrt(x*x+y*y)/z)     */
/*        t = r * cos(psi)                                               */
/*                                     psi = arctan(sqrt(x*x+y*y+z*z)/t) */
/*                                                                       */
/*        psi   in 0 to Pi             need psi over complete range      */
/*        phi   in 0 to Pi             need phi over complete range      */
/*        theta in 0 to 2Pi            need theta over complete range    */
/*        r     positive               any real x,y,z,t                  */
/*                                                                       */
/*        area = 8*Pi*r*r                                                */
/*        volume = 2 Pi*Pi*r*r*r                                         */
/*        4D volume = 1/2 Pi*Pi*r*r*r*r                                  */
/*                                                                       */
/*                                                                       */
/*  all angles are in radians, convert radians a to degrees d            */
/*                             d = a * 180.0/Pi                          */
/*      but, you may want the range in degrees to be different           */
/*                                                                       */
/*                             convert degrees d to radians a            */
/*                             a = d * Pi/180.0                          */

#include <stdio.h>
#include <math.h>
#define Pi  3.14159265358979323846

void cartesian_to_polar(double x, double y, double *r, double *theta)
{
  *r     = sqrt(x*x+y*y);
  *theta = atan2(y, x);
}

void polar_to_cartesian(double r, double theta, double *x, double *y)
{
  *x = r * cos(theta);
  *y = r * sin(theta);
}

void cartesian_to_cylindrical(double x, double y, double z,
                              double *r, double *theta, double *zs)
{
  *r     = sqrt(x*x+y*y);
  *theta = atan2(y, x);
  *zs    = z;
}

void cylindrical_to_cartesian(double r, double theta, double zs,
                              double *x, double *y, double *z)
{
  *x = r * cos(theta);
  *y = r * sin(theta);
  *z = zs;
}

void cartesian_to_spherical(double x, double y, double z,
                            double *r, double *theta, double *phi)
{
  *r     = sqrt(x*x+y*y+z*z);
  *theta = atan2(y, x);
  *phi   = atan2(sqrt(x*x+y*y),z);
}

void spherical_to_cartesian(double r, double theta, double phi,
                            double *x, double *y, double *z)
{
  *x = r * sin(phi) * cos(theta);
  *y = r * sin(phi) * sin(theta);
  *z = r * cos(phi);
}

void cartesian4D_to_spherical(double x, double y, double z, double t,
			    double *r, double *theta, double *phi, double *psi)
{
  *r     = sqrt(x*x+y*y+z*z+t*t);
  *theta = atan2(y, x);
  *phi   = atan2(sqrt(x*x+y*y),z);
  *psi   = atan2(sqrt(x*x+y*y+z*z),t);
}

void spherical_to_cartesian4D(double r, double theta, double phi, double psi,
			      double *x, double *y, double *z, double *t)
{
  *x = r * sin(psi) * sin(phi) * cos(theta);
  *y = r * sin(psi) * sin(phi) * sin(theta);
  *z = r * sin(psi) * cos(phi);
  *t = r * cos(psi);
}

double to_theta_degrees(double theta)
{ 
  /* latitude   positive  0 to 360 */
  if(theta>0.0)
    return theta * 180.0/Pi;
  else
    return 360.0 + theta * 180.0/Pi; 
}

double to_phi_degrees(double phi)
{
  /* longitude    -90 to +90 degrees */
  if(phi<Pi/2.0)
     return phi * 180.0/Pi;
  else
     return -(180.0 - phi * 180.0/Pi);
}

double to_psi_degrees(double psi)
{
   return psi * 180.0/Pi;
}

int main(int argc, char *argv[])
{
  double x1 = 0.5; /* first quadrant, first octant, first hypertant */
  double y1 = 0.7;
  double z1 = 0.8;
  double t1 = 0.9;
  double x2 = -0.8; /* third quadrant, seventh octant */
  double y2 = -0.5;
  double z2 = -0.7;
  double xz = 0.0;
  double yz = 0.0;
  double zz = 0.0;
  double x, y, z, t;
  double r, theta, phi, psi, zs;

  printf("coordinate.c running \n");
  printf("first quadrant \n");
  printf("convert  x1=%f, y1=%f \n", x1, y1);
  cartesian_to_polar(x1, y1, &r, &theta);
  printf("   to polar  r=%f, theta=%f \n\n", r, theta);
  printf("convert  r=%f, theta=%f \n", r, theta);
  polar_to_cartesian(r, theta, &x, &y);
  printf("   to cartesian  x=%f, y=%f \n", x, y);
  printf("        theta=%f degrees \n\n", to_theta_degrees(theta));

  printf("second quadrant \n");
  printf("convert  x2=%f, y1=%f \n", x2, y1);
  cartesian_to_polar(x2, y1, &r, &theta);
  printf("   to polar  r=%f, theta=%f \n\n", r, theta);
  printf("convert  r=%f, theta=%f \n", r, theta);
  polar_to_cartesian(r, theta, &x, &y);
  printf("   to cartesian  x=%f, y=%f \n", x, y);
  printf("        theta=%f degrees \n\n", to_theta_degrees(theta));

  printf("third quadrant \n");
  printf("convert  x2=%f, y2=%f \n", x2, y2);
  cartesian_to_polar(x2, y2, &r, &theta);
  printf("   to polar  r=%f, theta=%f \n\n", r, theta);
  printf("convert  r=%f, theta=%f \n", r, theta);
  polar_to_cartesian(r, theta, &x, &y);
  printf("   to cartesian  x=%f, y=%f \n", x, y);
  printf("        theta=%f degrees \n\n", to_theta_degrees(theta));

  printf("fourth quadrant \n");
  printf("convert  x1=%f, y2=%f \n", x1, y2);
  cartesian_to_polar(x1, y2, &r, &theta);
  printf("   to polar  r=%f, theta=%f \n\n", r, theta);
  printf("convert  r=%f, theta=%f \n", r, theta);
  polar_to_cartesian(r, theta, &x, &y);
  printf("   to cartesian  x=%f, y=%f \n", x, y);
  printf("        theta=%f degrees \n\n", to_theta_degrees(theta));

  printf("convert  x1=%f, y1=%f, z1=%f \n", x1, y1, z1);
  cartesian_to_cylindrical(x1, y1, z1, &r, &theta, &zs);
  printf("   to cylindrical  r=%f, theta=%f, zs=%f \n\n", r, theta, zs);
  printf("convert  r=%f, theta=%f, zs=%f \n", r, theta, zs);
  cylindrical_to_cartesian(r, theta, zs, &x, &y, &z);
  printf("   to cartesian  x=%f, y=%f, z=%f \n\n", x, y, z);

  printf("convert  x2=%f, y2=%f, z2=%f \n", x2, y2, z2);
  cartesian_to_cylindrical(x2, y2, z2, &r, &theta, &zs);
  printf("   to cylindrical  r=%f, theta=%f, zs=%f \n\n", r, theta, zs);
  printf("convert  r=%f, theta=%f, zs=%f \n", r, theta, zs);
  cylindrical_to_cartesian(r, theta, zs, &x, &y, &z);
  printf("   to cartesian  x=%f, y=%f, z=%f \n\n", x, y, z);

  printf("first octant \n");
  printf("convert  x1=%f, y1=%f, z1=%f \n", x1, y1, z1);
  cartesian_to_spherical(x1, y1, z1, &r, &theta, &phi);
  printf("   to spherical  r=%f, theta=%f, phi=%f \n\n", r, theta, phi);
  printf("convert  r=%f, theta=%f, phi=%f \n", r, theta, phi);
  spherical_to_cartesian(r, theta, phi, &x, &y, &z);
  printf("   to cartesian  x=%f, y=%f, z=%f \n", x, y, z);
  printf("        theta=%f degrees, phi=%f degrees \n\n",
         to_theta_degrees(theta), to_phi_degrees(phi));

  printf("seventh octant \n");
  printf("convert  x2=%f, y2=%f, z2=%f \n", x2, y2, z2);
  cartesian_to_spherical(x2, y2, z2, &r, &theta, &phi);
  printf("   to spherical  r=%f, theta=%f, phi=%f \n\n", r, theta, phi);
  printf("convert  r=%f, theta=%f, phi=%f \n", r, theta, phi);
  spherical_to_cartesian(r, theta, phi, &x, &y, &z);
  printf("   to cartesian  x=%f, y=%f, z=%f \n", x, y, z);
  printf("        theta=%f degrees, phi=%f degrees \n\n",
         to_theta_degrees(theta), to_phi_degrees(phi));

  printf("first hypertant \n");
  printf("convert  x1=%f, y1=%f, z1=%f, t1=%f \n", x1, y1, z1, t1);
  cartesian4D_to_spherical(x1, y1, z1, t1, &r, &theta, &phi, &psi);
  printf("   to spherical  r=%f, theta=%f, phi=%f, psi=%f \n\n",
         r, theta, phi, psi);

  printf("convert  r=%f, theta=%f, phi=%f, psi=%f \n", r, theta, phi, psi);
  spherical_to_cartesian4D(r, theta, phi, psi, &x, &y, &z, &t);
  printf("   to cartesian4D  x=%f, y=%f, z=%f, t=%f \n", x, y, z, t);
  printf("        theta=%f degrees, phi=%f degrees, psi=%f degrees \n\n",
         to_theta_degrees(theta), to_phi_degrees(phi), to_psi_degrees(psi));

  printf("finished coordinate.c \n");
  return 0;
} /* end main of coordinate.c */