69 const int& nrk_theta)
const {
81 cout <<
"Not yet prepared!!!" << endl ;
88 double xi_t0 = xi_i(0) ;
89 double xi_p0 = xi_i(1) ;
90 double xi_l0 = xi_i(2) ;
91 double theta0 = theta_i ;
93 double dt = - 0.5 * M_PI / double(nt-1) / double(nrk_theta) ;
110 double xi_t1, xi_t2, xi_t3, xi_t4, xi_tf ;
111 double xi_p1, xi_p2, xi_p3, xi_p4, xi_pf ;
112 double xi_l1, xi_l2, xi_l3, xi_l4, xi_lf ;
113 double f1, f2, f3, f4 ;
114 double g1, g2, g3, g4 ;
115 double h1, h2, h3, h4 ;
121 for (
int i=0; i<nrk_theta; i++) {
124 f1 = -2. * xi_p0 * dlnconfo.
val_point(rah, theta0, phi) ;
125 g1 = xi_l0 * rah * confo2.
val_point(rah, theta0, phi)
126 + 2. * xi_t0 * dlnconfo.
val_point(rah, theta0, phi) ;
127 h1 = - (1. - 2.*laplnconfo.
val_point(rah, theta0, phi)) * xi_p0
128 / rah / confo2.
val_point(rah, theta0, phi) ;
135 f2 = -2. * (xi_p0+0.5*xi_p1)
136 * dlnconfo.
val_point(rah, theta0+0.5*dt, phi) ;
137 g2 = (xi_l0+0.5*xi_l1) * rah
138 * confo2.
val_point(rah, theta0+0.5*dt, phi)
139 + 2. * (xi_t0+0.5*xi_t1)
140 * dlnconfo.
val_point(rah, theta0+0.5*dt, phi) ;
141 h2 = - (1. - 2.*laplnconfo.
val_point(rah, theta0+0.5*dt, phi))
142 * (xi_p0+0.5*xi_p1) / rah
143 / confo2.
val_point(rah, theta0+0.5*dt, phi) ;
150 f3 = -2. * (xi_p0+0.5*xi_p2)
151 * dlnconfo.
val_point(rah, theta0+0.5*dt, phi) ;
152 g3 = (xi_l0+0.5*xi_l2) * rah
153 * confo2.
val_point(rah, theta0+0.5*dt, phi)
154 + 2. * (xi_t0+0.5*xi_t2)
155 * dlnconfo.
val_point(rah, theta0+0.5*dt, phi) ;
156 h3 = - (1. - 2.*laplnconfo.
val_point(rah, theta0+0.5*dt, phi))
157 * (xi_p0+0.5*xi_p2) / rah
158 / confo2.
val_point(rah, theta0+0.5*dt, phi) ;
165 f4 = -2. * (xi_p0+xi_p3) * dlnconfo.
val_point(rah, theta0+dt, phi) ;
166 g4 = (xi_l0+xi_l3) * rah * confo2.
val_point(rah, theta0+dt, phi)
167 + 2. * (xi_t0+xi_t3) * dlnconfo.
val_point(rah, theta0+dt, phi) ;
168 h4 = - (1. - 2.*laplnconfo.
val_point(rah, theta0+dt, phi))
169 * (xi_p0+xi_p3) / rah / confo2.
val_point(rah, theta0+dt, phi) ;
177 xi_tf = xi_t0 + (xi_t1 + 2.*xi_t2 + 2.*xi_t3 + xi_t4) / 6. ;
178 xi_pf = xi_p0 + (xi_p1 + 2.*xi_p2 + 2.*xi_p3 + xi_p4) / 6. ;
179 xi_lf = xi_l0 + (xi_l1 + 2.*xi_l2 + 2.*xi_l3 + xi_l4) / 6. ;
190 xi_f.
set(0) = xi_tf ;
191 xi_f.
set(1) = xi_pf ;
192 xi_f.
set(2) = xi_lf ;
Tbl runge_kutta_theta(const Tbl &xi_i, const double &theta_i, const double &phi, const int &nrk) const
Compute a forth-order Runge-Kutta integration to the theta direction for the solution of the Killing ...
double val_point(double r, double theta, double phi) const
Computes the value of the field at an arbitrary point , by means of the spectral expansion.