Warp.hpp

//
// $Id: Warp.hpp,v 1.1 2004/08/12 16:21:33 werner Exp $
//
// $Log: Warp.hpp,v $
// Revision 1.1  2004/08/12 16:21:33  werner
// Integration of numerical geodesics now compiles. Working is not yet satisfying.
//
// Revision 1.2  2004/05/05 15:56:52  werner
// Separation of DOP core routines with dynamic size into the ODE library.
//
// 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 __Kerr_HPP
#define __Kerr_HPP "Created 27.02.2001 21:42:27 by werner"

#include <vecal/Christoffel.hpp>
#include <vecal/Geodesic.hpp>

namespace VecAl
{

struct  Kerr
{
        Scalar_t        m, a, e;

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

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

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

                g.set(0);

        S r2 = point[P4_R]*point[P4_R];
        S tmp = cos(point[P4_H]);

        S sig = r2 + a2*tmp*tmp;
        S del = r2 + a2 + e2 - 2*m*point[P4_R];

                tmp = sin(point[P4_H]);
        S tmp2 = tmp*tmp;

                g[M4_TT] = (a2*tmp2 -  del)/sig;
                g[M4_TP] = g[M4_PT] = -a*tmp2*(r2 + a2 - del)/sig;
                g[M4_RR] = sig/del;
                g[M4_HH] = sig;
                g[M4_PP] = ((r2 + a2)*(r2 + a2) - del*a2*tmp2)/sig*tmp2;
        }

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

                G.set(0);



        S r2 = point[P4_R]*point[P4_R];
        S tmpx = cos(point[P4_H]);

        S sig = r2 + a2*tmpx*tmpx;
        S del = r2 + a2 + e2 - 2*m*point[P4_R];

        S tmp = sin(point[P4_H]);
        S tmp2 = tmp*tmp;

        S t_b = point[P4_R]*point[P4_R] + a2;
        S t_c = t_b*t_b - del*a2*tmp2;
        S t_d = del - a2*tmp2;
        S t_e = t_b - del;
        S t_f = point[P4_R] - m;

        S t_n = t_d*t_c + a2*tmp2*t_e*t_e;



        chr[C4_TTR] =
           chr[C4_TRT] = (t_c*(t_f - point[P4_R]*t_d/sig)
                          + a2*tmp2*t_e*(m - point[P4_R]*t_e/sig))/t_n;
        chr[C4_TTH] =
           chr[C4_THT] = a2*tmp*tmpx/t_n*(-t_c + t_d*t_c/sig
                                          + t_e*t_e*(1 + a2*tmp2/sig));
        chr[C4_TRP] =
           chr[C4_TPR] = a*tmp2/t_n*(m*t_c - t_e*(2*point[P4_R]*t_b
                                                  - a2*tmp2*t_f));
        chr[C4_TPH] =
           chr[C4_THP] = a2*a*del*tmp2*tmp*tmpx*t_e/t_n;

        chr[C4_RTT] = (sig*t_f - point[P4_R]*t_d)/sig/sig/sig*del;
        chr[C4_RTP] =
           chr[C4_RPT] = ((sig*m - point[P4_R]*t_e)/sig/sig/sig
                          *del*a*tmp2);
        chr[C4_RRR] = (del*point[P4_R] - sig*t_f)/del/sig;
        chr[C4_RRH] =
           chr[C4_RHR] = -a2*tmp*tmpx/sig;
        chr[C4_RPP] = (-t_b*2*point[P4_R] + t_f*a2*tmp2
                       + t_c*point[P4_R]/sig)*del*tmp2/sig/sig;
        chr[C4_RHH] = -del*point[P4_R]/sig;
        
        
        chr[C4_PTR] =
           chr[C4_PRT] = -a/t_n*(m*t_d - t_f*t_e);
        chr[C4_PTH] =
           chr[C4_PHT] = -a*tmpx*t_e/t_n*(t_d/tmp + a2*tmp);
        chr[C4_PRP] =
           chr[C4_PPR] = (t_d*(2*point[P4_R]*t_b - a2*tmp2*t_f
                               - point[P4_R]*t_c/sig)
                          +a2*tmp2*t_e*(m - point[P4_R]*t_e/sig))/t_n;
        chr[C4_PPH] =
           chr[C4_PHP] = (t_d*tmpx*(t_c*(1/tmp + a2*tmp/sig)
                                    - a2*del*tmp)
                          + a2*tmp*tmpx*t_e*t_e*(1 + a2*tmp2/sig))/t_n;

        chr[C4_HTT] = -(sig - del + a2*tmp2)/sig/sig/sig*a2*tmp*tmpx;
        chr[C4_HTP] =
           chr[C4_HPT] = tmp*tmpx*t_e/sig/sig*(1 + a2*a*tmp2/sig);
        chr[C4_HRR] = a2*tmp*tmpx/del/sig;
        chr[C4_HRH] =
           chr[C4_HHR] = point[P4_r]/sig;
        chr[C4_HPP] = -tmp*tmpx/sig/sig*(t_b*t_b - 2*del*a2*tmp2
                                         + t_c*a2*tmp2/sig);
        chr[C4_HHH] = -a2*tmp*tmpx/sig;

        }
};

} /* namespace VecAl */ 

#endif /* __Kerr_HPP */