NumericalGeodesics3D.cpp

Compute geodesics in a 3D tensor field.

Compute geodesics in a 3D tensor field.This example sets up a 10x10x10 metric tensor field with the unit metric describing Euclidean space, but sets an artificial anisotropy at a single point within the data field. Then, a geodesic is computed that will be a straight line far from this singular point, but deviated close to it. For demonstration purposes, the geodesic is computed via Euler integration.

#include <iostream>
 
#include <spacetime/Schwarzschild.hpp>
#include <spacetime/Geodesic853.hpp>
#include <spacetime/Acceleration.hpp>
#include <vecal/CartesianChart.hpp>
#include <spacetime/NumericalSpacetime.hpp>
#include <spacetime/IntegrateGeodesic.hpp>
 
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
 
using namespace VecAl;
using namespace Traum;
using namespace std;
 
void NumericalTestField3D()
{
// symbolic names for coordinates
        enum { x=0, y=1, z=2 };
 
// typed tensorial object:
typedef LowerTriangular<3, double> myMetric;
 
// The data storage of the tensor field is an automatic variable.
myMetric TensorStorage[1000];
 
// This tensor field, consisting of 1000 points, is defined over 
// a three-dimensional base manifold with 10 vertices in each dimension.
MultiArray<3, myMetric, myMetric*> G(TensorStorage, MIndex(10,10,10) );
 
// Initialize all elements through the MultiArray interface
myMetric InitValue;
        InitValue(0,0) = InitValue(1,1) = InitValue(2,2) = 1;
        InitValue(0,1) = InitValue(1,2) = InitValue(2,0) = 0;           
        G.set( InitValue );
 
        // Modify the data values at some point of interest
        G[ MIndex(5,6,7) ] (2,2) = 2;
 
        // Starting point
VVector<3, double> q0, q;
        q = 5.1, 6.8, 5.03;
        q0 = q;
 
        // Velocity
Row<3, double> v0, v;
        v = 0.01,0.02,.2;
        v0 = v;
ChristoffelField<MultiArray<3, myMetric, myMetric*> > Chris(G, q );
 
clock_t start, finish;
        start = clock();
        // Integration
        {for(int is = 0; is<40000; is++)
        {
                // path parameter increment
        double  ds = 0.001;
 
                // coordinate acceleration, \ddot q
        Row<3, double> a;
 
                Chris.getCoordinateAcceleration(a, v);
 
                // Euler step
                q += v*ds;
                v += a*ds;
 
                if (is %1000 == 0)
                        cout << "v="<<v<<"\t q="<<q<<endl;
        }}
        finish = clock();
        printf( "Geodesic computation: %2.1f seconds\n", (double)(finish - start) / CLOCKS_PER_SEC );
}