SchwarzschildGeodesic.hpp

//
// $Id: SchwarzschildGeodesic.hpp,v 1.12 2007/03/06 23:01:16 werner Exp $
//
// $Log: SchwarzschildGeodesic.hpp,v $
// Revision 1.12  2007/03/06 23:01:16  werner
// Integration of the Vector library into the Eagle library, for the sake
// of type registration.
//
// Revision 1.11  2007/02/22 17:08:42  werner
// A major revision of the vector/fixed array interface, unifying three independent
// libraries into one. The vecal (Vector Algebra), fish/vector (Multidimensional
// vector arrays) and the Eagle library (Geometric Algebra) thus share a common
// root now.
//
// Revision 1.10  2005/08/02 18:16:57  werner
// Using new metric interface.
//
// Revision 1.9  2004/09/23 11:22:49  werner
// Little improvements.
//
// Revision 1.8  2004/09/15 21:43:46  werner
// Removed include file.
//
// Revision 1.7  2004/09/15 16:49:41  werner
// Various bugfixes and functionality fixes.
//
// Revision 1.6  2004/09/03 12:26:15  werner
// Adjusted moved header files.
//
// Revision 1.5  2004/08/24 09:09:53  werner
// Towards VC7
//
// Revision 1.4  2004/08/12 16:21:33  werner
// Integration of numerical geodesics now compiles. Working is not yet satisfying.
//
// Revision 1.3  2004/05/14 10:01:24  werner
// Towards position output, temporary commit.
//
// Revision 1.2  2004/05/13 12:19:48  werner
// Adjustments for gcc 3.4.0.
//
// Revision 1.1  2004/05/13 12:07:53  werner
// Schwarzschild geodesics via one-dimensional equatorial solution.
//
#ifndef __SCHWARZSCHILD_GEODESIC_HPP
#define __SCHWARZSCHILD_GEODESIC_HPP "Created 12.05.2004 21:57:45 by werner benger"

#include <ode/Integrate853.hpp>
//#include <vecal/PolarChart.hpp>
//#include <vecal/CartesianChart.hpp>

#include <eagle/STA.hpp>

namespace Traum
{
using namespace Eagle;

#ifdef  DB_LOG
using namespace wb;
#endif

class   EquatorialPlane
{
public: typedef long double real;
        typedef Eagle::STA::sphericalpoint      PolarPt;
        typedef Eagle::STA::point               CartesianPt;
        typedef Eagle::STA::vector              CartesianVec;

        CartesianPt     Center;
        CartesianVec    X, Y;

// Koordinatentransformation Schwarzschild --> kartesisch
        CartesianPt operator()(const PolarPt&P) const
        {
        real    r  = P[P.R],
                phi= P[P.PHI];

                return CartesianPt( Center + r * ( sin(phi) * Y - cos(phi) * X) );
        }

static  real cosine(const CartesianVec&l, const CartesianVec&r)
        {
                return dot( l.spatial(), r.spatial() );
        }

// Koordinatentransformation kartesisch --> Schwarzschild
        PolarPt operator()(const CartesianPt&P,  real&sinphi, real&cosphi) const
        {
        CartesianVec p = P-Center;
//      CartesianVec p = difference(P, Center);
        /*
        Sei p der normalisierte Richtungsvektor zum Punkt,  dann gilt:

                p = sin(phi) * Y - cos(phi) * X

                p*X =           -cos(phi)
                p*Y = sin(phi) 
         */
        PhysicalSpace::vector dist_vec = p.spatial();

        real    r       = norm( dist_vec );
                cosphi  = -cosine(X,p) / r;
                sinphi  =  cosine(Y,p) / r;
                
        real    phi     = atan2(sinphi, cosphi);

                if (phi<0)
                        phi += 2*M_PI; 

        PolarPt R;
                R[R.T] = P[P.T];
                R[R.R] = r;
                R[R.PHI] = phi;
                R[R.THETA] = 0;
                return R;
        }

#if 0

/*
  Tangentialvektor einer Geodäten (u,u') am Punkt phi berechnen
  @returns u^2 dr(phi)/dphi
 */
vector3 CoordinateAxes::direction(double sinphi, double cosphi,
                                  double u, double du) const
{
        return  (-du * sinphi + u * cosphi) * Y +
                ( du * cosphi + u * sinphi) * X;

#if 0
/*
Geradengleichung in radialinversen Polarkoordinaten:

        u (phi) = sin(phi)/d - k/d  cos(phi)
        u'(phi) = cos(phi)/d + k/d  sin(phi)

Daraus: d = (sin(phi)-k*cos(phi))/u = (cos(phi)+k*sin(phi))/u'

Damit:  z = k*n

Mit     */
double  z = sinphi/u  - cosphi/du,
        n = sinphi/du + cosphi/u;
/*
Nun ist
        dir = X + k*Y = X + z/n Y

und daher:
        n*dir = n*X + z*Y;
 */ 
        return  n*X - z*Y; // wo kommt das '-' her??? 
//      return  n*X + z*Y; 
#endif
}

#endif

};


class   SchwarzschildGeodesic : public IntegratePath853<long double>
{
public: typedef long double real;
        real    M;
        EquatorialPlane EQ;

        typedef STA::sphericalpoint     PolarPt;
        typedef STA::sphericalvector    PolarVec;
        typedef STA::point              CartesianPt;
        typedef STA::vector             CartesianVec;

        typedef PhysicalSpace::point    point3;
        typedef PhysicalSpace::vector   vector3;

        SchwarzschildGeodesic(const real&mass=0)
        : M(mass)
        {
                EQ.Center.set(0.0);
        }

/*
  Eine Gerade ist in Polarkoordinaten gegeben durch:

        r(phi) = d / [ sin(phi) - k cos(phi) ]

bzw.
        u(phi) = 1/r(phi) = [sin(phi) - k cos(phi) ] / d

Durch die Wahl der Strahlrichtung als X-Achse wird die Steigung
in der Geraden in diesem Koordinatensystem k=0. Der Parameter d ist
der Minimalabstand der Geraden vom Koordinatenursprung.
 Mit dieser Wahl der X,Y - Achse ergibt sich:

        r(phi) = d / sin(phi)

        u(phi) = sin(phi) / d

        u'(phi) = cos(phi) / d

  
*/
        bool restart(const real&s, const CartesianPt&x0, const CartesianVec&v0)
        {
        enum { N = 1, i=0 }; 


        point3   q0 = x0.spatial();
        vector3 dq0 = v0.spatial();
                dq0 /= norm(dq0);

        const point3 & Center = EQ.Center.spatial(); 

#ifdef  DB_LOG
                db_logf(21, "SchwarzschildGeodesic: EQ Axis: Ray (%lg,%lg,%lg) + l (%lg,%lg,%lg) \n", 
                        double( q0[0]), double( q0[1]), double( q0[2]),
                        double(dq0[0]), double(dq0[1]), double(dq0[2])
                        );
#endif

        // 
        // Nm ist derjenige Punkt auf der Geraden, der dem Massenzentrum 
        // am nächsten liegt. Der Richtungsvektor von diesem Punkt zum 
        // Massenzentrum bildet somit die Y-Achse. Die X-Achse wird 
        // durch den Sichtstrahl selbst gebildet. 
        //
        // Geradenparameter für Punktabstand
        //  l = (P-O)*d / (d*d)
        //
        real    l = dot(Center - q0, dq0);

//      cout << "Closest linear approach at lambda="<< l << endl; 
#ifdef  DB_LOG
        db_logf(21, "SchwarzschildGeodesic: EQ Axis: Closest linear approach is at lambda=%lg\n", double(l) );
#endif

        point3  Nm = q0 + l*dq0; 
//      point3  Nm = q0 + l*dq0;

                EQ.X = tvector4( 0, dq0);

        vector3 y = Nm-Center; 
//      vector3 y = difference(Nm, Center); 

#ifdef  DB_LOG
#define SVECP(V) double(V[V->x]),double(V[V->y]),double(V[V->z])

        db_logf(21,"SchwarzschildGeodesic: CLOSEPOINT: (%lg,%lg,%lg) (lambda=%lg [%lg])\n", 
                                                     SVECP(y),     double(l), double(l/M) );
#endif
        
        // 
        // d ist der Minimale Abstand des Sichtstrahles vom Massenzentrum 
        //   dh. der Impact Parameter 
        //
//      real    d = q - Center; 
        real    d = norm(y);
                //cout << "Impact parameter " << d << endl;

                if (d < 1E-8)
                {
                real *u  = initialize(2*N,s);
                real *du = u + N;
                         u[i] = 0; 
                        du[i] = 0;

#ifdef  DB_LOG 
                        db_logf(21,"SchwarzschildGeodesic: DIRECT VIEW with impact parameter %lg (relative %lg, mass %lg), cannot setup axis!\n",
                                double(d), double(d/M), double(M) );
#endif
                        return false;
                }

                y /= d;
                EQ.Y = tvector4(0, y);

#ifdef  DB_LOG 
                db_logf(21,"SchwarzschildGeodesic: AXIS: (%lg,%lg,%lg) x (%lg,%lg,%lg) \n",
                        SVECP(EQ.X), SVECP(EQ.Y) ); 
#endif


        real sinphi, cosphi;
        PolarPt p0 = EQ(x0, sinphi, cosphi); 
#ifdef  DB_LOG
        db_logf(21,"SchwarzschildGeodesic: Initial point at polar coordinates (r=%lg,theta=%lg,phi=%lg)\n", 
                double(p0[p0->r]),double(p0[p0->th]),double(p0[p0->ph]) );
#endif

        real    phi = p0[p0.PHI]; 
        real    *u  = initialize(2*N, phi);
        real    *du = u + N;
        
                 u[i] = sinphi/d;       // u (phi0) 
                du[i] = cosphi/d;       // u'(phi0) 
#ifdef  DB_LOG 
                db_logf(21,"SchwarzschildGeodesic: sinphi=%lg  cosphi=%lg d=%lg\n",
                        double(sinphi), double(cosphi), double(d) );

                //cout << "u(0)= " << u[i] << "  u'(0)=" << du[0] << endl;
#endif 
                return true;
        }

        void Accel(real, const real *x, const real *v, real *d2x_ds2)
        {
                for(unsigned int i=0; i<nPaths(); i++)
                {
                        d2x_ds2[i] = 3*M*x[i]*x[i]-x[i]; // ddy == y''=3my^2-y
                }
        }

        CartesianPt     position() const
        {
        PolarPt P; 
                P = EQ.Center[P.T], 1/q(0), 0, s1();

                return EQ( P );
        }


        CartesianVec velocity() const
        {
        real    phi = s1(),
                  u = q(0),
                 du = q(1),
             sinphi = sin(phi),
             cosphi = cos(phi);

        double  r_sin_delta = du * sinphi - u * cosphi,
                r_cos_delta =  u * sinphi +du * cosphi; 

                return r_cos_delta * EQ.X + r_sin_delta * EQ.Y; 
        }

        CartesianPt     position(real phi) const
        {
        PolarPt P; 
                P = EQ.Center[P.T],
                    1/q(0, phi),
                    0,
                    phi;

                return EQ( P );
        }


        CartesianVec velocity(real phi) const
        {
        real      u = q(0, phi),
                 du = q(1, phi),
             sinphi = sin(phi),
             cosphi = cos(phi);

        double  r_sin_delta = du * sinphi - u * cosphi,
                r_cos_delta =  u * sinphi +du * cosphi; 

                return r_cos_delta * EQ.X + r_sin_delta * EQ.Y; 
        }


};

} // namespace Vecal

#endif // __SCHWARZSCHILD_GEODESIC_HPP