/*******************************************************************
GEODES/GEODTEST.C   Geodesic Computation Examples    August 19, 1996

This routine includes example geodesic problems. It calls the
geodes routine that approximates a minimum-length geodesic.

********************************************************************

A geodesic problem is specified by:
    * manifold representation in rectangular coordinates.
    * ambient space metric (symmetric) that manifold inherits.
    * geodesic endpoints.

This routine calls the geodesic computing routine in geodes.c.  It
passes two callback functions that specify the manifold and metric.
The manifold function computes position vector, Jacobian matrix, and
Hessian matrix.  The metric function computes the ambient space
metric.  Both Euclidean metric and the non-Euclidean hyperbolic
metric are in the examples.

Some geodes routine arguments control performance. See geodes
routine for a description of all arguments.

********************************************************************

This software can be modified or extended by:
    * Allowing indefinite metric with null geodesic(s).
    * Optimizing grid points dynamically. Using recursion.
    * Optimizing performance control values dynamically.
    * Finding families of geodesics.
    * Add localized salients (bumps) to the manifold.
    * Plotting geodesic in both parameter and ambient space.

Refer comments and questions to:
    William L. Anderson
    Elements Research
    7501 Windyrush Road
    Charlotte, NC 28226
    704-543-9180
    elements@ix.netcom.com
    http://www.netcom.com/~elements/

********************************************************************

Copyright (c) 1994-1996 Elements Research, William L. Anderson.

Permission to use, copy, modify, and distribute this software and
its documentation for any purpose without fee is hereby granted,
provided that the above copyright notice appear in all copies and
that both that copyright notice and this permission notice appear
in supporting documentation.

This software is provided "as is" without express or implied
warranty.

*******************************************************************/

#include <stdio.h>
#include <math.h>
#include <float.h>
#include <malloc.h>

/* C language function prototypes */
void ntorus_init(int n_arg, double *a_arg, double *b_arg);
void ntorus_free();
void ntorus(double *x, double *y, double *jacobian, double *hessian);
void ncylinder(double *x, double *y, double *jacobian, double 
*hessian);
void ambmet_matrix_init(int m_arg);
void ambmet_matrix_free();
void ambmet_euclidean(double *y, double **ambient_metric_matrix_arg);
void ambmet_hyperbolic(double *y, double **ambient_metric_matrix_arg);
void geodes(int n, int m,
    void yeval(double *,double *,double *,double *),
    void ambmet(double *,double **),
    int pmax, double accel, double eps, int itmax,
    double *xend, double *xgeod);

/* macros reduce multidimensional array to C language array */
#define SUBS2(n,i,j) ((i)*(n)+(j))
#define SUBS3(n,m,i,j,k) (((i)*(n)+(j))*(m)+(k))

static double random_value;

static void frandoi(double seed) {
random_value = seed;
random_value = fabs(random_value);
random_value = log(random_value)*1.0e+9;
}

static double frandom() {
static const double two_power_31_minus_one = 2147483647.0;
static const double two_power_31 = 2147483648.0;
random_value = fmod(16807.0 * random_value,two_power_31_minus_one);
return (random_value/two_power_31);
}

static void printgeod(int n, double *xgeod,
int pmax) {
int p,i;
for (p = 0; p < pmax; p++) {
    for (i = 0; i < n; i++)
        printf("%g\t",xgeod[SUBS2(n,p,i)]);
    printf("\n");
    }
}

static void nsphere_geod(int n,void ambmet(double *,double **)) {
double *a;
double *b;
double *xend;
double *xgeod;
int pmax;
double accel;
double eps;
int itmax;
int i,r;
printf("Geodesic on unit %d-sphere\n",n);
a = malloc(sizeof(double)*(n+1));
b = malloc(sizeof(double)*n);
for (r = 0; r < n+1; r++)
    a[r] = 1.0;
for (i = 0; i < n; i++)
    b[i] = 0.0;
pmax = 100;
accel = 0.5;
eps = 1.0e-4;
itmax = 50;
xend = malloc(sizeof(double)*(2*n));
xend[SUBS2(n,0,0)] = 0.5;
xend[SUBS2(n,0,1)] = 1.0;
xend[SUBS2(n,1,0)] = 1.0;
xend[SUBS2(n,1,1)] = 0.5;
for (i = 2; i < n; i++) {
    xend[SUBS2(n,0,i)] = frandom();
    xend[SUBS2(n,1,i)] = frandom();
    }
xgeod = malloc(sizeof(double)*(pmax+2)*n);
ambmet_matrix_init(n+1);
ntorus_init(n,a,b);
geodes(n,n+1,ntorus,ambmet,pmax,accel,eps,itmax,xend,xgeod);
ntorus_free();
ambmet_matrix_free();
printgeod(n,xgeod,pmax+2);
free(a);
free(b);
free(xend);
free(xgeod);
}

static void ncylinder_geod(int n,void ambmet(double *,double **)) {
/* geodesic on cylinder (2D cylinder is isometric with plane) */
double *a;
double *b;
double *xend;
double *xgeod;
int pmax;
double accel;
double eps;
int itmax;
int i,r;
printf("Geodesic on unit %d-cylinder\n",n);
a = malloc(sizeof(double)*(n+1));
b = malloc(sizeof(double)*n);
for (r = 0; r < n+1; r++)
    a[r] = 1.0;
for (i = 0; i < n; i++)
    b[i] = 0.0;
pmax = 100;
accel = 0.5;
eps = 1.0e-4;
itmax = 50;
xend = malloc(sizeof(double)*(2*n));
xend[SUBS2(n,0,0)] = 0.5;
xend[SUBS2(n,0,1)] = 1.0;
xend[SUBS2(n,1,0)] = 1.0;
xend[SUBS2(n,1,1)] = 0.5;
for (i = 2; i < n; i++) {
    xend[SUBS2(n,0,i)] = frandom();
    xend[SUBS2(n,1,i)] = frandom();
    }
xgeod = malloc(sizeof(double)*(pmax+2)*n);
ambmet_matrix_init(n+1);
ntorus_init(n,a,b);
geodes(n,n+1,ncylinder,ambmet,pmax,accel,eps,itmax,xend,xgeod);
ntorus_free();
ambmet_matrix_free();
printgeod(n,xgeod,pmax+2);
free(a);
free(b);
free(xend);
free(xgeod);
}

static void nellipsoid_geod(int n,void ambmet(double *,double **)) {
double *a;
double *b;
double *xend;
double *xgeod;
int pmax;
double accel;
double eps;
int itmax;
int i,r;
printf("Geodesic on %d-Ellipsoid\n",n);
a = malloc(sizeof(double)*(n+1));
b = malloc(sizeof(double)*n);
for (r = 0; r < n+1; r++)
    a[r] = (double)r + 1.0;
for (i = 0; i < n; i++)
    b[i] = 0.0;
pmax = 100;
accel = 0.5;
eps = 1.0e-4;
itmax = 50;
xend = malloc(sizeof(double)*(2*n));
xend[SUBS2(n,0,0)] = 0.5;
xend[SUBS2(n,0,1)] = 1.0;
xend[SUBS2(n,1,0)] = 1.0;
xend[SUBS2(n,1,1)] = 0.5;
for (i = 2; i < n; i++) {
    xend[SUBS2(n,0,i)] = frandom();
    xend[SUBS2(n,1,i)] = frandom();
    }
xgeod = malloc(sizeof(double)*(pmax+2)*n);
ambmet_matrix_init(n+1);
ntorus_init(n,a,b);
geodes(n,n+1,ntorus,ambmet,pmax,accel,eps,itmax,xend,xgeod);
ntorus_free();
ambmet_matrix_free();
printgeod(n,xgeod,pmax+2);
free(a);
free(b);
free(xend);
free(xgeod);
}

static void ntorus_geod(int n,void ambmet(double *,double **)) {
double *a;
double *b;
double *xend;
double *xgeod;
int pmax;
double accel;
double eps;
int itmax;
int i,r;
printf("Geodesic on unit %d-torus\n",n);
a = malloc(sizeof(double)*(n+1));
b = malloc(sizeof(double)*n);
for (r = 0; r < n+1; r++)
    a[r] = 1.0;
for (i = 0; i < n; i++)
    b[i] = 1.0;
pmax = 100;
accel = 0.5;
eps = 1.0e-4;
itmax = 50;
xend = malloc(sizeof(double)*(2*n));
xend[SUBS2(n,0,0)] = 0.5;
xend[SUBS2(n,0,1)] = 1.0;
xend[SUBS2(n,1,0)] = 1.0;
xend[SUBS2(n,1,1)] = 0.5;
for (i = 2; i < n; i++) {
    xend[SUBS2(n,0,i)] = frandom();
    xend[SUBS2(n,1,i)] = frandom();
    }
xgeod = malloc(sizeof(double)*(pmax+2)*n);
ambmet_matrix_init(n+1);
ntorus_init(n,a,b);
geodes(n,n+1,ntorus,ambmet,pmax,accel,eps,itmax,xend,xgeod);
ntorus_free();
ambmet_matrix_free();
printgeod(n,xgeod,pmax+2);
free(a);
free(b);
free(xend);
free(xgeod);
}

static void ntorus_random_geod(int n,void ambmet(double *,double **)) {
double *a;
double *b;
double *xend;
double *xgeod;
int pmax;
double accel;
double eps;
int itmax;
int i,r;
printf("Geodesic on random radius %d-torus\n",n);
a = malloc(sizeof(double)*(n+1));
b = malloc(sizeof(double)*n);
for (r = 0; r < n+1; r++)
    a[r] = frandom();
for (i = 0; i < n; i++)
    b[i] = frandom() + 1.0;
pmax = 100;
accel = 0.5;
eps = 1.0e-4;
itmax = 50;
xend = malloc(sizeof(double)*(2*n));
xend[SUBS2(n,0,0)] = 0.5;
xend[SUBS2(n,0,1)] = 1.0;
xend[SUBS2(n,1,0)] = 1.0;
xend[SUBS2(n,1,1)] = 0.5;
for (i = 2; i < n; i++) {
    xend[SUBS2(n,0,i)] = frandom();
    xend[SUBS2(n,1,i)] = frandom();
    }
xgeod = malloc(sizeof(double)*(pmax+2)*n);
ambmet_matrix_init(n+1);
ntorus_init(n,a,b);
geodes(n,n+1,ntorus,ambmet,pmax,accel,eps,itmax,xend,xgeod);
ntorus_free();
ambmet_matrix_free();
printgeod(n,xgeod,pmax+2);
free(a);
free(b);
free(xend);
free(xgeod);
}

static void test_suite(void ambmet(double *,double **)) {
ncylinder_geod(2,ambmet);
ncylinder_geod(3,ambmet);
nsphere_geod(2,ambmet);
nsphere_geod(3,ambmet);
nsphere_geod(6,ambmet);
nellipsoid_geod(2,ambmet);
nellipsoid_geod(3,ambmet);
nellipsoid_geod(6,ambmet);
ntorus_geod(2,ambmet);
ntorus_geod(3,ambmet);
ntorus_geod(6,ambmet);
ntorus_geod(10,ambmet);
ntorus_random_geod(2,ambmet);
ntorus_random_geod(3,ambmet);
ntorus_random_geod(6,ambmet);
ntorus_random_geod(10,ambmet);
}

void main() {
frandoi(acos(-1.0));
printf("Ambient metric is Euclidean\n");
test_suite(ambmet_euclidean);
printf("Ambient metric is Hyperbolic\n");
test_suite(ambmet_hyperbolic);
}