Schwarzschild.hpp

//
// $Id: Schwarzschild.hpp,v 1.9 2007/02/22 23:11:52 werner Exp $
//
// $Log: Schwarzschild.hpp,v $
// Revision 1.9  2007/02/22 23:11:52  werner
// Using Eagle library instead of Vecal
//
// Revision 1.8  2004/09/03 12:26:15  werner
// Adjusted moved header files.
//
// Revision 1.7  2004/08/12 16:21:33  werner
// Integration of numerical geodesics now compiles. Working is not yet satisfying.
//
// Revision 1.6  2004/05/13 16:44:23  werner
// Formulate Christoffel symbols via direct vectorization.
//
// Revision 1.5  2004/05/11 16:53:47  werner
// Introduction of coordinate systems for improved type-safety.
// Will support easy coordinate transformations soon.
//
// Revision 1.4  2004/05/06 22:42:16  werner
// Towards a specification of a spacetime via the Acceleration structure.
//
// Revision 1.3  2004/05/03 13:33:33  werner
// integration improved
//
// Revision 1.2  2004/03/29 11:51:02  werner
// Common interface among simple integrators, DiffMe and Vecal, and preliminiary work on integrating dop853.
//
// Revision 1.1  2004/03/22 11:55:02  werner
// Schwarzschild geodesic integration.
//
// Revision 1.1  2004/02/13 16:36:21  werner
// Initial preliminiary version of the Vector Algebra Library.
//
#ifndef __Schwarzschild_HPP
#define __Schwarzschild_HPP "Created 27.02.2001 21:42:27 by werner"

#include "Geodesic.hpp"
#include <eagle/STA.hpp>

namespace Traum
{
        using namespace Eagle;

struct  Schwarzschild : Chart<STA::SphericalChart4D>
{
        typedef Chart<STA::SphericalChart4D> Polar4Dd;

        using Polar4Dd::Point_t;
        using Polar4Dd::Scalar_t;
        using Polar4Dd::Vector_t;
        using Polar4Dd::Metric;
        using Polar4Dd::Christoffel;

        enum
        {
                t = Polar4Dd::T,
                h = Polar4Dd::THETA,
                p = Polar4Dd::PHI,
                r = Polar4Dd::R
        };

        Scalar_t        m;

        Schwarzschild(const Scalar_t&mass)
        : m(mass)
        {}

        void getMetric(Metric&g, const Point_t&P) const
        {
        Scalar_t sinTheta = sin(P[h]);

                g.set(0);
                g(t,t) = 1 - 2*m/P[r];
                g(r,r) = -1/g(t,t);
                g(h,h) = -P[r]*P[r];
                g(p,p) =  g(h,h)*sinTheta*sinTheta ;
        }

        void getChristoffel(Christoffel_t&G, const Point_t&P) const
        {
        Scalar_t sinTheta = sin(P[h]); 
                G.set(0.0); 
                G(t,t,r) = G(t,r,r)  = m/P[r]/(P[r] - 2*m);
                G(r,t,t) = m*(1 - 2*m/P[r]) / (P[r] * P[r] );
                G(r,r,r) = -G(t,t,r);
                G(r,h,h) = 2*m - P[r];
                G(r,p,p) = G(r,h,h)*sinTheta*sinTheta;
                G(h,r,h) = G(h,h,r)= 1/P[r];
                G(h,p,p) = -sinTheta*cos( P[h] );
                G(p,r,p) = G(p,p,r) = G(h,r,h);
                G(p,h,p) = G(p,p,h) = 1/tan( P[h] );
        }
};

} /* namespace VecAl */ 

#endif /* __Schwarzschild_HPP */