/******************************************************************* 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 #include #include #include /* 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); }