NumericalGeodesics4D.cpp

Compute geodesics in a 4D tensor field.

Compute geodesics in a 4D tensor field.This example sets up a 10x10x10x10 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.

This example is similar to NumericalGeodesics3D.cpp, but just extended to four dimensions. In the same manner, geodesics can be computed in arbitrary dimensions. However, the computational increases enormously.

#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;
 
// Compute geodesics in a 4D tensor field
void NumericalTestField4D()
{
// symbolic names for coordinates
        enum { x=0, y=1, z=2, t=3 };
 
// typed tensorial object:
typedef LowerTriangular<4, double> Tensor;
 
// The data storage of the tensor field is an automatic variable.
Tensor TensorStorage[10000];
 
// This tensor field, consisting of 1000 points, is defined over 
// a three-dimensional base manifold with 10 vertices in each dimension.
MultiArray<4, Tensor, Tensor*> G(TensorStorage, MIndex(10,10,10,10) );
 
// Initialize all elements through the MultiArray interface
Tensor  InitValue;
        InitValue.set(0);
        InitValue(0,0) = InitValue(1,1) = InitValue(2,2) = InitValue(3,3) = 1;
        G.set( InitValue );
 
        // Modify the data values at some point of interest
        G[ MIndex(5,6,7,3) ] (2,2) = 2;
 
        // Starting point
VVector<4, double> q0, q;
        q = 5.1, 6.8, 5.03, 3;
        q0 = q;
 
        // Velocity
Row<4, double> v0, v;
        v = 0.01,0.02,.2,0;
        v0 = v;
ChristoffelField<MultiArray<4, Tensor, Tensor*> > Chris(G,q);
 
clock_t start, finish;
        start = clock();
        {for(int is = 0; is<40000; is++)
        {
                // path parameter increment
        double  ds = 0.001;
 
                // coordinate acceleration, \ddot q
        Row<4, 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( "4D Geodesic computation: %2.1f seconds\n", (double)(finish - start) / CLOCKS_PER_SEC );
}