SphericalSymmetric.hpp

//
// $Id: SphericalSymmetric.hpp,v 1.2 2004/05/11 16:53:47 werner Exp $
//
// $Log: SphericalSymmetric.hpp,v $
// Revision 1.2  2004/05/11 16:53:47  werner
// Introduction of coordinate systems for improved type-safety.
// Will support easy coordinate transformations soon.
//
// 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 <vecal/Christoffel.hpp>
#include <vecal/Geodesic.hpp>

namespace VecAl
{

struct  SphericalSymmetric
{
        enum { t, r, h, p,
               theta=h, phi = p
        };

        Scalar_t A(Scalar_t r, Scalar_t t);
        Scalar_t B(Scalar_t r, Scalar_t t);


        void metric(Metric<Scalar_t,4>&g, const Point_t&P) const
        {
        Scalar_t sinTheta = sin(P[h]);

                g.set(0);
                g(t,t) =  A(P[r], P[t]);
                g(r,r) = -B(P[r], P[t]);
                g(h,h) = -P[r]*P[r];
                g(p,p) = -P[r]*P[r]*sinTheta*sinTheta;
        }

        void cometric(Metric<Scalar_t,4>&g, const Point_t&P) const
        {
        Scalar_t sinTheta = sin(P[h]);

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


        void getChristoffel(Christoffel<Scalar_t,4>&G, const Point_t&P) const
        {
        Scalar_t sinTheta = sin(P[h]);

                G.set(0);
                G(t,t,t) = dotA(P[r], P[t]) / (2*A(P[r], P[t]) );
                /*
                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_ph = cos(theta)/sin(theta)
                G(p,h,p) = G(p,p,h) = 1/tan( P[h] );
                */
        }


        void Ricci(Metric<Scalar_t,4>&g, const Point_t&P) const
        {}
};

} /* namespace VecAl */ 

#endif /* __Schwarzschild_HPP */