/*******************************************************************
GEODES/METRIC.C   Ambient Space Metric Computation   August 19, 1996

These routines compute the ambient space metric matrix at a set of
coordinates y^r.

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

The C language function prototypes are:
    void ambmet_matrix_init(int m);
    void ambmet_matrix_free();
    void ambmet_euclidean(double *y, double **ambmet_matrix);
    void ambmet_hyperbolic(double *y, double **ambmet_matrix);

Arguments are:
    m ...... Dimension of ambient space y^r.
    y ...... Position vector in ambient rectangular coordinates.
    ambmet_matrix.. Address of address of matrix of ambient
                    space metric values.

The ambmet_matrix_init() routine makes the evaluation routines as
efficient as possible.  The setup routine presets the ambient
space dimension and allocates an ambient metric matrix.  The
ambmet_matrix_free() routine frees this matrix.  The evaluation
routines compute the ambient metric matrix and sets a pointer to
it.  This is more efficient than coping a matrix and encapsulates
the matrix in this ambient metric routine.  This is important
since the ambient space dimension is independent of the manifold
embedded in it.

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

The ambient space dimension is n+1.

Two ambient metrics are provided for: the Euclidean metric and the
non-Euclidean hyperbolic metric.

An ambient metric matrix must be symmetric.

Reference:
    Manfredo do Carmo, Riemannian Geometry, Birkhauser, Boston, 1992

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

This software can be modified or extended by:
    * Including other metrics.
    * Allow ambient space dimension greater than n+1.

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>

/* 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 int m;

static double *ambmet_matrix;

void ambmet_matrix_init(int m_arg) {
int r,s;
m = m_arg;
ambmet_matrix = malloc(sizeof(double)*(m*m));
/* initialize as euclidean metric */
for (r = 0; r < m; r++)
    for (s = 0; s < m; s++)
        if (r == s)
            ambmet_matrix[SUBS2(m,r,s)] = 1.0;
        else
            ambmet_matrix[SUBS2(m,r,s)] = 0.0;
}

void ambmet_matrix_free() {
free(ambmet_matrix);
ambmet_matrix = NULL;
m = 0;
}

void ambmet_euclidean(double *y, double **ambmet_matrix_arg) {
*ambmet_matrix_arg = ambmet_matrix;
}

void ambmet_hyperbolic(double *y, double **ambmet_matrix_arg) {
int r;
double hyperbolic;
hyperbolic = 1.0/(y[0]*y[0]);
for (r = 0; r < m; r++)
    ambmet_matrix[SUBS2(m,r,r)] = hyperbolic;
*ambmet_matrix_arg = ambmet_matrix;
}