Geodesics.hpp

00001 #ifndef __GEODESCICS_HPP_
00002 #define __GEODESCICS_HPP_
00003 
00004 #include <fish/fiber/vector/CubicIpol.hpp>
00005 #include <fish/fiber/vector/LinearIpol.hpp>
00006 #include <eagle/QuadraticMatrix.hpp> 
00007 #include <eagle/Coordinates.hpp>
00008 #include <fish/fiber/baseop/LocalFromWorldPoint.hpp>
00009 #include <ode/dop853.hpp>
00010 
00011 //#define VERBOSE
00012 
00013 namespace Fiber
00014 {
00015 
00016 template <>
00017 struct  LinearIpolZeroDerivativeTrait<metric33>
00018 {
00019 static  metric33 zero()
00020         {
00021         metric33 m;
00022                 m.set(0.0);
00023                 return m;
00024         }
00025 };
00026 
00027 template <>
00028 struct  LinearIpolZeroDerivativeTrait<metric44>
00029 {
00030 static  metric44 zero()
00031         {
00032         metric44 m;
00033                 m.set(0.0);
00034                 return m;
00035         }
00036 };
00037 
00038 //template <>
00039 // struct  CubicIpolZeroDerivativeTrait<metric33>
00040 // {
00041 // static       metric33 zero()
00042 //         {
00043 //         metric33 m;
00044 //              m.set(0.0);
00045 //                 return m;
00046 //         }
00047 // };
00048 
00049 // template <>
00050 // struct  CubicIpolZeroDerivativeTrait<metric44>
00051 // {
00052 // static       metric44 zero()
00053 //         {
00054 //         metric44 m;
00055 //              m.set(0.0);
00056 //                 return m;
00057 //         }
00058 // };
00059 
00060 }
00061 
00062 
00063 namespace FieldLines
00064 {
00065 
00066 struct Geodesic{};
00067 
00068 struct GeodesicSpatial{};
00069 
00070 template<class T>
00071 double getDet(const Eagle::Quadratic<T::Dims,double>&q)
00072 {
00073         return  Eagle::Determinantor<T::Dims>::compute(q);
00074 }
00075 
00076 template<class T>
00077 class ChristoffelXYZ : public MemCore::ReferenceBase< ChristoffelXYZ<T> >
00078 {
00079         typedef typename T::Chart_t Chart_t;
00080         typedef typename Eagle::Coordinates<Chart_t>::point point;
00081         // store metric tensor
00082         T g;
00083         
00084         Eagle::Quadratic<T::Dims, double> g_inv;
00085  
00086         // store precomputed finite differences of g // array size 3
00087         T g_diff[T::Dims];
00088 
00089 
00090         bool ok;
00091         // consts for access
00092 //      enum { MU = 0, NU, LAMBDA  };
00093 //      enum { X = 0, Y, Z};
00094 
00095 public:
00096         ChristoffelXYZ( const T&gP, const T&g_xP, const T&g_yP, const T&g_zP,
00097                         const double dxP, const double dyP, const double dzP  )
00098         : MemCore::ReferenceBase< ChristoffelXYZ >(this)
00099         , g(gP)
00100         , ok(true)
00101         {
00102                 std::cout << "aber jetzt nit im ersnst!!!";
00103                 if(T::Dims == 3)
00104                 {
00105                         enum { X = 0, Y, Z};
00106 
00107                 if(dxP != 0.0)
00108                         for(unsigned i = 0; i < T::SIZE; i++)
00109                                 g_diff[X][i] = ( g_xP[i] - g[i] ) / dxP; 
00110 
00111                 if(dyP != 0.0)
00112                         for(unsigned i = 0; i < T::SIZE; i++)
00113                                 g_diff[Y][i] = ( g_yP[i] - g[i] ) / dyP; 
00114 
00115                 if(dzP != 0.0)
00116                         for(unsigned i = 0; i < T::SIZE; i++)
00117                                 g_diff[Z][i] = ( g_zP[i] - g[i] ) / dzP;
00118 
00119                 } 
00120 
00121                 if(T::Dims == 4)
00122                 {
00123                         enum {  t = Eagle::STA::CartesianChart4D<>::T, 
00124                                 X = Eagle::STA::CartesianChart4D<>::X,
00125                                 Y = Eagle::STA::CartesianChart4D<>::Y,
00126                                 Z = Eagle::STA::CartesianChart4D<>::Z }; 
00127 
00128                         if(dxP != 0.0)
00129                                 for(unsigned i = 0; i < T::SIZE; i++)
00130                                         g_diff[X][i] = ( g_xP[i] - g[i] ) / dxP; 
00131 
00132                         if(dyP != 0.0)
00133                                 for(unsigned i = 0; i < T::SIZE; i++)
00134                                         g_diff[Y][i] = ( g_yP[i] - g[i] ) / dyP; 
00135 
00136                         if(dzP != 0.0)
00137                                 for(unsigned i = 0; i < T::SIZE; i++)
00138                                         g_diff[Z][i] = ( g_zP[i] - g[i] ) / dzP;
00139                         
00140                 }
00141         }
00142 
00143         ChristoffelXYZ( const T&gP, const T&g_dxP, const T&g_dyP, const T&g_dzP )
00144         : MemCore::ReferenceBase< ChristoffelXYZ<T> >(this)
00145         , g(gP)
00146         , ok(true)      
00147         {
00148                 //std::cout << g_dxP << endl << g_dyP << endl << g_dzP << endl;
00149                 //std::cout << "------";
00150                 if(T::Dims == 3)
00151                 {
00152                         enum { X = 0, Y, Z};
00153                         g_diff[X] = g_dxP; 
00154                         g_diff[Y] = g_dyP; 
00155                         g_diff[Z] = g_dzP; 
00156                 } 
00157                 if(T::Dims == 4)
00158                 {
00159                         enum {  t = Eagle::STA::CartesianChart4D<>::T, 
00160                                 X = Eagle::STA::CartesianChart4D<>::X,
00161                                 Y = Eagle::STA::CartesianChart4D<>::Y,
00162                                 Z = Eagle::STA::CartesianChart4D<>::Z }; 
00163 
00164                         g_diff[t].set(0.0); 
00165                         g_diff[X] = g_dxP; 
00166                         g_diff[Y] = g_dyP; 
00167                         g_diff[Z] = g_dzP; 
00168                 }
00169                         try
00170                         {
00171                                 Eagle::Quadratic<T::Dims, double> q (g); 
00172 
00173                                 Eagle::ComputeInverse(g_inv, q); 
00174 
00175                                 Eagle::Matrix<T::Dims,T::Dims, double> &mg = q ; 
00176                                 Eagle::Matrix<T::Dims,T::Dims, double> &mg_i = g_inv; 
00177                                 Eagle::Matrix<T::Dims,T::Dims, double> I = mg * mg_i; 
00178                                 //std::cout << "controll inversion: " << I << std::endl; 
00179                         } 
00180                         catch(const Eagle::DegeneratedMatrix&)
00181                         {
00182                                 puts("ChristoffelXYZ<>::ChristoffelXYZ ERROR Degenerated Matrix in inversion."); fflush(stdout); 
00183                                 ok = false;
00184                         } 
00185         }
00186 
00187         bool isOk(){ return ok;}
00188 
00189         void  getMetricfromQuadratic(const Eagle::Quadratic<T::Dims,double>&q, T&g)
00190                 {
00191                         for(unsigned i = 0; T::SIZE; i++)
00192                                 g[i] = q[i];
00193                 }
00194 
00195 
00196 
00197         // christoffel symbols function for on demand symbol computing, later there may be the need to 
00198         // store them and reuse them when possible.
00199         const double operator()(const unsigned mu, const unsigned nu, const unsigned lambda) const 
00200         {
00201         double gamma = 0.0; 
00202 
00203                 if(mu > T::Dims-1 || nu > T::Dims-1 || lambda > T::Dims-1)
00204                         return 0.0; 
00205 
00206                 for(unsigned i = 0; i < T::Dims; i++)
00207                         gamma += g_inv(lambda, i) * ( g_diff[nu](mu,i) + g_diff[mu](nu,i) - g_diff[i](mu, nu) ); 
00208 
00209                 gamma /= 2;
00210 
00211                 return gamma;
00212         }
00213         
00214         const double gammaDebug(const unsigned mu, const unsigned nu, const unsigned lambda) const 
00215         {
00216         double gamma = 0.0; 
00217 
00218                 if(mu > T::Dims-1 || nu > T::Dims-1 || lambda > T::Dims-1)
00219                         return 0.0; 
00220                         printf( "mu nu lam %d %d %d", mu, nu, lambda);
00221                 for(unsigned i = 0; i < T::Dims; i++)
00222                 {
00223                         
00224                 double tmp = g_inv(lambda, i) * ( g_diff[nu](mu,i) + g_diff[mu](nu,i) - g_diff[i](mu, nu) );
00225                     gamma += tmp; 
00226                 printf( "g_inv: %f, g1: %f, g2:% f, g3: %f, term: %f\n", g_inv(lambda, i), g_diff[nu](mu,i), g_diff[mu](nu,i), g_diff[i](mu, nu), tmp);
00227                 }
00228                 gamma /= 2; 
00229                 printf( "gamma: %f\n", gamma);
00230                 return gamma;
00231         }
00232 
00233 
00234         const T& getG() { return g; }
00235         
00236         void speak() const
00237        {
00238                
00239                std::cout << "Christoffel content:\n" << "g " << endl << g << "g_inv: " << endl << g_inv; 
00240 
00241                if(T::Dims == 3)
00242                        std::cout << "g,x" << std::endl <<g_diff[0] 
00243                                  << "g,y" << std::endl <<g_diff[1] 
00244                                  << "g,z" << std::endl <<g_diff[2]; 
00245                                
00246                if(T::Dims == 4)
00247                        std::cout << "g,t" << std::endl <<g_diff[0] 
00248                                  << "g,x" << std::endl <<g_diff[1] 
00249                                  << "g,y" << std::endl <<g_diff[2] 
00250                                  << "g,z" << std::endl <<g_diff[3]; 
00251 
00252                cout << "Gamma (T)|X|Y|Z|:" << endl;
00253 
00254                for( unsigned lam = 0; lam <  T::Dims; lam++)
00255                 {
00256                         for( unsigned mu = 0; mu <  T::Dims; mu++)
00257                         {
00258                                 for( unsigned nu = 0; nu <  T::Dims; nu++)
00259                                 {
00260                                         printf("%f ", (*this)(mu,nu,lam));
00261                                 } 
00262                                 printf("\n");
00263                         } 
00264                         printf("-----------------------------\n");
00265                 }
00266 
00267 
00268                 // for( unsigned lam = 0; lam <  T::Dims; lam++)
00269                 // {
00270                 //      for( unsigned mu = 0; mu <  T::Dims; mu++)
00271                 //      {
00272                 //              for( unsigned nu = 0; nu <  T::Dims; nu++)
00273                 //              {
00274                 //                      gammaDebug(mu, nu, lam);        
00275                 //              } 
00276                 //              printf("\n");
00277                 //      } 
00278                 //      printf("-----------------------------\n");
00279                 // }
00280        }
00281         
00282 };
00283 
00284 
00287 template<class T>
00288 typename Eagle::Coordinates<typename T::Chart_t >::vector getQddot( RefPtr< ChristoffelXYZ<T> >& Gamma, 
00289                                                                     typename Eagle::Coordinates<typename T::Chart_t >::vector&q_dot )
00290 {
00291 typename Eagle::Coordinates<typename T::Chart_t >::vector q_ddot; 
00292 
00293         q_ddot.set(0.0);
00294 
00295         // q-ddot^lamba = - Gamma^lambda_mu_nu * q-dot^mu * q-dot^nu 
00296 
00297         //cout << "getQddot  qdot:" << q_dot << endl;
00298 
00299         for(index_t lambda = 0; lambda < T::Dims; lambda++)
00300                 for( index_t mu = 0; mu < T::Dims; mu++)
00301                         for(index_t nu = 0; nu < T::Dims; nu++)
00302                                 q_ddot[lambda] -= (*Gamma)( mu, nu, lambda ) * q_dot[mu] * q_dot[nu]; 
00303 
00304         //cout << "getQddot  qddot:" << q_ddot << endl;
00305 
00306         return q_ddot;
00307 }
00308 
00309 
00310 template<class TT>
00311 typename Eagle::Coordinates<typename TT::Chart_t >::vector getQdot( const typename Eagle::Coordinates<typename TT::Chart_t >::vector& q_dot, 
00312                                                                     const TT&g, const bool big = true)
00313 {
00314         if (TT::Dims==3)
00315                 return q_dot; 
00316 
00317         enum {  T = STA::CartesianChart4D<>::T, 
00318                 X = STA::CartesianChart4D<>::X,
00319                 Y = STA::CartesianChart4D<>::Y,
00320                 Z = STA::CartesianChart4D<>::Z }; 
00321 
00322 const double eps = 1.0e-4; 
00323 
00324         //std::cout << "qdot before:" << q_dot << std::endl;
00325 
00326 typename Eagle::Coordinates<typename TT::Chart_t >::vector v = q_dot; 
00327 
00328         if( !isnormal(v[X]) ) v[X] = 0.0; 
00329         if( !isnormal(v[Y]) ) v[Y] = 0.0; 
00330         if( !isnormal(v[Z]) ) v[Z] = 0.0; 
00331                 
00332 
00333         if(v[X] < eps && v[X] > -eps &&
00334            v[Y] < eps && v[Y] > -eps &&
00335            v[Z] < eps && v[Z] > -eps)
00336         {
00337                 //puts("yep");
00338                 v[T] = 1.0; 
00339                 return v;
00340         } 
00341 
00342         //puts("hm");
00343 
00344 double A = 2 * ( g(X,T)*v[X]    + g(Y,T)*v[Y]      + g(Z,T)*v[Z] ); 
00345 double B = g(X,X)*v[X]*v[X]     + g(Y,Y)*v[Y]*v[Y] + g(Z,Z)*v[Z]*v[Z]+
00346            2*( g(X,Y)*v[X]*v[Y] + g(X,Z)*v[X]*v[Z] + g(Y,Z)*v[Y]*v[Z] ); 
00347 
00348         if(g(T,T) == 0.0)
00349         {
00350                 puts("getQdot::Error g(T,T) == 0! Division by 0."); 
00351                 return v;
00352         }
00353 
00354 double p = A/g(T,T); 
00355 double q= B/g(T,T); 
00356 
00357 double w = p*p/4 - q; 
00358 
00359         if(w < 0.0)
00360         {
00361                 puts("getQdot::Error direction initialization srqt < 0!!!!"); 
00362                 return v;
00363         }
00364 
00365         double sw = sqrt( w ); 
00366 
00367         double v1 = -p/2 + sw; 
00368         double v2 = -p/2 - sw; 
00369 
00370                 if(big)
00371                         v[T] = (v1 > v2)? v1 : v2; // use bigger value of the solution of the quadratic equation 
00372                 else
00373                         v[T] = (v1 > v2)? v2 : v1; 
00374 
00375         return v;
00376 } 
00377 
00378 
00381 template<class T>
00382 bool getChristoffelXYZ( const typename Eagle::Coordinates<typename T::Chart_t >::point& CurrentPoint, 
00383                         const MemCore::RefPtr<LocalFromWorldPoint>&LocalPointFinder, 
00384                         const FieldSelector&FieldSelection, const double step_size, 
00385                         RefPtr< ChristoffelXYZ<T> >& Christoffels,
00386                         T& CurrentMetric,
00387                         pair<Eagle::PhysicalSpace::point, string>& localPoint ) 
00388         {
00389         typedef typename T::Chart_t Chart_t;
00390         typedef typename Coordinates<Chart_t>::point point; 
00391 
00392                 Eagle::PhysicalSpace::point current (CurrentPoint[0+T::Dims-3], CurrentPoint[1+T::Dims-3], CurrentPoint[2+T::Dims-3]); 
00393 
00394         bool test = LocalPointFinder->get( current, localPoint ); 
00395 
00396                 if(!test)
00397                 {
00398                         puts("getChristoffelXYZ()1 error finding local point at:"); 
00399                         cout << current << endl;
00400                         return false;
00401                 } 
00402                 
00403         const string&FragName = localPoint.second; 
00404         const Eagle::PhysicalSpace::point &float_index = localPoint.first; 
00405 
00406 
00407 
00408         RefPtr<Field> Field; 
00409                 Field = FieldSelection.getField(); 
00410 
00411 #ifdef VERBOSE
00412                 cout << "FieldLines::getChristoffelXYZ local point: " << float_index; 
00413                 cout << "in Fragment:" << FragName << endl;
00414                 cout << "getting Field: " << FieldSelection.getFieldName() << endl;
00415 #endif          
00416                 if(!Field) 
00417                 { 
00418                         puts("AtomicIntegrator<metric33, Geodesic>doit2(); ERROR: No Field0 found in Time iteration"); 
00419                         return false; 
00420                 } 
00421 
00422 
00423          RefPtr<FragmentID> FragID;
00424 
00425                 if(FragName.compare("") != 0 )
00426                         FragID   = new FragmentID( FragName ); 
00427 
00428                 RefPtr<CreativeArrayBase> cBase = Field->getCreator( FragID ); 
00429                 assert(cBase); 
00430 
00431                 RefPtr<MemBase> mBase = cBase->create();                
00432                 assert(mBase); 
00433 
00434 
00435 
00436                 RefPtr<MemArray<3, T> > ToIntArr0 = mBase; 
00437                 ToIntArr0   = Field->getCreator( FragID )->create(); 
00438 
00439                 if(!ToIntArr0)
00440                 {
00441                         cout << "getChristoffelXYZ<T> dims:" << T::Dims << endl; 
00442                         mBase.speak("getChristoffelXYZ<T> ------>");
00443 
00444                         puts("getChristoffelXYZ()<T>:ERROR retrieving matching field. Maybe Dims of instantiated metric and Dims of field data does not match. 33 44?"); 
00445                         return false;   
00446                 } 
00447                 
00448                 assert( ToIntArr0 );
00449 
00450 
00451 
00452                         // NOTE: Here optionally use another interpolator, eventually cubic, or 
00453                         //       self-aligning for eigenvector fields.
00454                         // 
00455 
00456 //              cout << "before interpol" << endl; 
00457 
00458         T g_x, g_y, g_z; 
00459         T g_xl, g_xp, g_xxp, g_yp, g_zp, g_yzp, g_yyzp; 
00460 
00461         Interpolate<3, T, CubicIpol<T> > myField0  ( *ToIntArr0,   float_index); 
00462         Interpolate<3, T, CubicIpol<T>, double, NoDelimiter<T>, 0  > myField0_x( *ToIntArr0, float_index); 
00463         Interpolate<3, T, CubicIpol<T>, double, NoDelimiter<T>, 1  > myField0_y( *ToIntArr0, float_index); 
00464         Interpolate<3, T, CubicIpol<T>, double, NoDelimiter<T>, 2  > myField0_z( *ToIntArr0, float_index); 
00465         // Interpolate<3, T, LinearIpol<T> > myField0  ( *ToIntArr0,   float_index); 
00466 //      Interpolate<3, T, LinearIpol<T>, double, NoDelimiter<T>, 0  > myField0_xl( *ToIntArr0, float_index); 
00467         // Interpolate<3, T, LinearIpol<T>, double, NoDelimiter<T>, 1  > myField0_y( *ToIntArr0, float_index); 
00468         // Interpolate<3, T, LinearIpol<T>, double, NoDelimiter<T>, 2  > myField0_z( *ToIntArr0, float_index); 
00469         
00470                   CurrentMetric = myField0.eval(); 
00471                   g_x = myField0_x.eval(); 
00472                   g_y = myField0_y.eval(); 
00473                   g_z = myField0_z.eval(); 
00474 
00475 //                g_xl= myField0_xl.eval();
00476 
00477           // double eps = -1; 
00478 
00479           // Eagle::point3 xp =  float_index + Eagle::tvector3(eps,0,0 ); 
00480           // Eagle::point3 xxp =  float_index + Eagle::tvector3(-eps,0,0 ); 
00481 
00482           // Eagle::point3 yp =  float_index + Eagle::tvector3(0,eps,0 ); 
00483           // Eagle::point3 zp =  float_index + Eagle::tvector3(0,0,eps ); 
00484           // Eagle::point3 yzp = float_index + Eagle::tvector3(0,0,eps ) + Eagle::tvector3(0,eps,0 ); 
00485           // Eagle::point3 yyzp = float_index + Eagle::tvector3(0,0,eps ) - Eagle::tvector3(0,eps,0 ); 
00486 
00487           //Interpolate<3, T, LinearIpol<T> >  myField0_xp( *ToIntArr0, xp ); 
00488           //Interpolate<3, T, LinearIpol<T> >  myField0_xxp( *ToIntArr0, xxp ); 
00489           //Interpolate<3, T, LinearIpol<T> >  myField0_yp( *ToIntArr0, yp ); 
00490           //Interpolate<3, T, LinearIpol<T> >  myField0_zp( *ToIntArr0, zp ); 
00491           // Interpolate<3, T, LinearIpol<T> >  myField0_yzp( *ToIntArr0, yzp ); 
00492           // Interpolate<3, T, LinearIpol<T> >  myField0_yyzp( *ToIntArr0, yyzp ); 
00493 
00494           //         g_xp  = myField0_xp.eval(); 
00495           //        g_xxp  = myField0_xxp.eval(); 
00496 
00497           //      g_yp  = myField0_yp.eval(); 
00498           //      g_zp  = myField0_zp.eval(); 
00499                   // g_yzp = myField0_yzp.eval(); 
00500                   // g_yyzp = myField0_yyzp.eval(); 
00501 
00502                   RefPtr<Fiber::Field> field = FieldSelection.getCartesianPositions(); 
00503                   RefPtr<Fiber::DirectProductMemArray<Eagle::point3,3> > FData = field->getData(); 
00504                   Eagle::point3 delta = FData->Delta(); 
00505                   
00506 
00507         Eagle::Quadratic<T::Dims, double> q (CurrentMetric); 
00508         double det = getDet<T>(q); 
00509 
00510                 
00511 
00512                 if(det < 0.0 && T::Dims == 3) // check if det gets negative 
00513                 {
00514                         cout << "g:" << q << endl;
00515                         cout << "det: " << det << endl;
00516                         puts("getChristoffelXYZ()<T>:warning, Det of G < 0."); 
00517                         return false;
00518                 } 
00519 
00520                 g_x  /= delta[0]; // normalize index derivative by interval length  
00521                 g_xl /= delta[0]; 
00522                 g_y  /= delta[1]; 
00523                 g_z  /= delta[2]; 
00524 
00525 #ifdef VERBOSE 
00526 /*
00527         cout << "det(g): " << det << endl; 
00528         cout << "delta: " << delta << endl;
00529         cout << "aus dem interpolator:" << endl << CurrentMetric << g_x << g_y << g_z; 
00530 
00531                 cout << endl <<"g:" << endl 
00532                      << float_index << endl << CurrentMetric << "variiert xyz:" << endl 
00533                      << xp  << endl << g_xp  
00534                      << xxp  << endl << g_xxp  
00535                         
00536                      << yp  << endl << g_yp
00537                      << zp  << endl << g_zp 
00538 
00539                      << "g,x lin" << endl << g_xl;
00540                      // << yzp << endl << g_yzp  
00541                      // << yyzp << endl << g_yyzp  ; 
00542 */
00543 #endif          
00544                   Christoffels = new ChristoffelXYZ<T>( CurrentMetric, g_x, g_y, g_z ); 
00545 
00546 #ifdef VERBOSE
00547                   std::cout << "Hilfe Christoffel wo bist du?" << std::endl; 
00548                   cout << "current point: " << CurrentPoint;
00549                   Christoffels->speak();
00550 #endif  
00551                   if( Christoffels->isOk() == false)
00552                   {
00553                           puts("getChristoffelXYZ()<T>:warning, christoffels probably degenerated matrix."); 
00554                           return false;
00555                   }
00556 
00557                   return true;
00558         }
00559 
00560 template<class T>
00561 struct  GeodesicDifferentialEquation : public Traum::VirtualSecondOrderDiffEquation <double>
00562 {
00563 typedef typename T::Chart_t Chart_t;
00564 typedef typename Eagle::Coordinates<Chart_t>::point point;
00565 typedef typename Eagle::Coordinates<Chart_t>::vector tvector;
00566 
00567 RefPtr<LocalFromWorldPoint> LocalPointFinder;   
00568 FieldSelector FieldSelection;
00569 
00570 double step_size; // depricated, define stepsize for finite differential derivation of the metric tensors used for ChristoffelsXYZ 
00571         
00572 double scale;
00573 
00574 bool valid_point;
00575 
00576 pair<Eagle::PhysicalSpace::point, std::string>localPoint;
00577 
00578 T CurrentMetric;
00579 tvector DataDir;
00580 RefPtr<ChristoffelXYZ<T> > Gamma;
00581 
00582         enum { N = T::Dims * 2 };
00583         typedef double real;
00584         
00585         GeodesicDifferentialEquation()
00586         {}
00587         
00588         virtual ~GeodesicDifferentialEquation() {}
00589         
00590         struct RealArray_t : FixedArray<real, N>               
00591         {
00592                 void init(int) {}
00593         };
00594         
00595         struct DOPVarsArray_t : FixedArray<Traum::dop_vars<real>, N >
00596         {
00597                 void init(int) {}
00598         };
00599         
00600         void setData(const tvector&DataDirP, const RefPtr<LocalFromWorldPoint>&LocalPointFinderP, const FieldSelector&FieldSelectionP,
00601                      const double step_sizeP, const double scaleP)
00602         {
00603                 DataDir = DataDirP; 
00604                 LocalPointFinder = LocalPointFinderP; 
00605                 FieldSelection = FieldSelectionP; 
00606                 step_size = step_sizeP; 
00607                 scale = scaleP;
00608         }
00609 
00610         void Accel(real s, const real *x, const real *v, real *d2x_ds2)
00611         {
00612                 point Q; 
00613                 tvector q_dot; 
00614 
00615                         for(int k=0; k < point::Dims; k++)
00616                         {
00617                                 Q[k] = x[k]; 
00618                                 q_dot[k] = v[k]; 
00619                                 DataDir[k] = v[k];
00620                         } 
00621 //              printf("Q: %f, %f, %f\n", Q[0], Q[1], Q[2] ); 
00622 //              printf("V: %f, %f, %f\n", v[0], v[1], v[2] ); 
00623 
00624         bool test = FieldLines::getChristoffelXYZ<T>( Q, LocalPointFinder, FieldSelection, step_size, 
00625                                                      Gamma,
00626                                                      CurrentMetric,
00627                                                      localPoint ); 
00628                 if(!test)
00629                 {
00630                         puts("GeodesicDifferentialEquation::DifferentialEquation::Accel() No valid point found!");fflush(stdout); 
00631                         return;
00632                 } 
00633 
00634         tvector qdd = getQddot<T>( Gamma, q_dot); 
00635 
00636                 for(int k=0; k < point::Dims; k++)
00637                         d2x_ds2[k] = qdd[k];
00638         }
00639 
00640 };
00641 
00642 // double linearInterpolate(const double t0, const double d0, const double t1, const double d1, const double t)
00643 // {
00644 //      return d0 + (d1 - d0) / (t1 - t0) * (t - t0);
00645 // }
00646 
00647 // metric33 linearInterpolate(const double t0, const metric33&m0, const double t1, const metric33&m1, const double t )
00648 // {
00649 // metric33 tmp; 
00650 //      for(unsigned i = 0; i < metric33::SIZE; i++)
00651 //              tmp[i] = linearInterpolate(t0, m0[i], t1, m1[i], t); 
00652 
00653 //      return tmp;
00654 // }
00655 
00656 
00657 }
00658 #endif // __GEODESCICS_HPP_