67 const double& theta_i,
const int& nrk_phi,
68 const int& nrk_theta)
const {
86 cout <<
"Not yet prepared!!!" << endl ;
95 double dp = 2. * M_PI / double(np) ;
106 xi_t.
set(0) = xi_i(0) ;
107 xi_p.
set(0) = xi_i(1) ;
108 xi_l.
set(0) = xi_i(2) ;
115 xi_ini.
set(0) = xi_i(0) ;
116 xi_ini.
set(1) = xi_i(1) ;
117 xi_ini.
set(2) = xi_i(2) ;
119 double pp_0 = phi_i ;
121 for (
int i=1; i<np+1; i++) {
125 xi_t.
set(i) = xi(0) ;
126 xi_p.
set(i) = xi(1) ;
127 xi_l.
set(i) = xi(2) ;
136 cout <<
"Hole_bhns::killing_vect :" << endl ;
138 cout <<
" xi_p(0) : " << xi_p(0) << endl ;
139 cout <<
" xi_p(np) : " << xi_p(np) << endl ;
150 for (
int k=0; k<np; k++) {
172 double hh = 2. * M_PI / double(nn) ;
176 double t1, t2, t3, t4, t5 ;
184 for (
int i=0; i<mm; i++) {
186 t1 = hh * double(4*i) ;
187 t2 = hh * double(4*i+1) ;
188 t3 = hh * double(4*i+2) ;
189 t4 = hh * double(4*i+3) ;
190 t5 = hh * double(4*i+4) ;
192 integ += (hh/45.) * (14.*source_phi.
val_point(rah,M_PI/2.,t1)
193 + 64.*source_phi.
val_point(rah,M_PI/2.,t2)
194 + 24.*source_phi.
val_point(rah,M_PI/2.,t3)
195 + 64.*source_phi.
val_point(rah,M_PI/2.,t4)
196 + 14.*source_phi.
val_point(rah,M_PI/2.,t5)
201 cout <<
"Hole_bhns:: t_f = " << integ << endl ;
202 double ratio = 0.5 * integ / M_PI ;
204 cout <<
"Hole_bhns:: t_f / 2M_PI = " << ratio << endl ;
206 for (
int k=0; k<np; k++) {
214 double dt = 0.5 * M_PI / double(nt-1) ;
225 for (
int i=0; i<np; i++) {
227 xi_th.
set(nt-1, i) = xi_t(i) ;
228 xi_ph.
set(nt-1, i) = xi_p(i) ;
229 xi_ll.
set(nt-1, i) = xi_l(i) ;
233 for (
int i=0; i<np; i++) {
236 xi_ini.
set(0) = xi_t(i) ;
237 xi_ini.
set(1) = xi_p(i) ;
238 xi_ini.
set(2) = xi_l(i) ;
240 double pp = double(i) * dp ;
241 double tt_0 = theta_i ;
243 for (
int j=1; j<nt; j++) {
247 xi_th.
set(nt-1-j, i) = xi(0) ;
248 xi_ph.
set(nt-1-j, i) = xi(1) ;
249 xi_ll.
set(nt-1-j, i) = xi(2) ;
260 cout <<
"Hole_bhns::killing_vect :" << endl ;
262 cout <<
" xi_ph(nt-1,0) : " << xi_ph(nt-1,0) << endl ;
263 cout <<
" xi_ph(0,0) : " << xi_ph(0,0) << endl ;
264 cout <<
" xi_ph(0,np/4) : " << xi_ph(0,np/4) << endl ;
265 cout <<
" xi_ph(0,np/2) : " << xi_ph(0,np/2) << endl ;
266 cout <<
" xi_ph(0,3np/4) : " << xi_ph(0,3*np/4) << endl ;
274 killing.
set(1) = 0.5 ;
276 killing.
set(2) = 0.5 ;
277 killing.
set(3) = 0.5 ;
279 for (
int l=0; l<2; l++) {
280 for (
int i=0; i<nr; i++) {
281 for (
int j=0; j<nt; j++) {
282 for (
int k=0; k<np; k++) {
283 (killing.
set(1)).set_grid_point(l, k, j, i) = xi_ll(j, k) ;
284 (killing.
set(2)).set_grid_point(l, k, j, i) = xi_th(j, k) ;
285 (killing.
set(3)).set_grid_point(l, k, j, i) = xi_ph(j, k) ;
295 double check_norm1 = 0. ;
296 double check_norm2 = 0. ;
297 source_phi =
pow(
confo_tot, 2.) * rr * st / killing(3) ;
300 for (
int i=0; i<mm; i++) {
302 t1 = hh * double(4*i) ;
303 t2 = hh * double(4*i+1) ;
304 t3 = hh * double(4*i+2) ;
305 t4 = hh * double(4*i+3) ;
306 t5 = hh * double(4*i+4) ;
308 check_norm1 += (hh/45.) *
309 ( 14.*source_phi.
val_point(rah,M_PI/4.,t1)
310 + 64.*source_phi.
val_point(rah,M_PI/4.,t2)
311 + 24.*source_phi.
val_point(rah,M_PI/4.,t3)
312 + 64.*source_phi.
val_point(rah,M_PI/4.,t4)
313 + 14.*source_phi.
val_point(rah,M_PI/4.,t5) ) ;
315 check_norm2 += (hh/45.) *
316 ( 14.*source_phi.
val_point(rah,M_PI/8.,t1)
317 + 64.*source_phi.
val_point(rah,M_PI/8.,t2)
318 + 24.*source_phi.
val_point(rah,M_PI/8.,t3)
319 + 64.*source_phi.
val_point(rah,M_PI/8.,t4)
320 + 14.*source_phi.
val_point(rah,M_PI/8.,t5) ) ;
325 cout <<
"Hole_bhns:: t_f for M_PI/4 = " << check_norm1 / M_PI
327 cout <<
"Hole_bhns:: t_f for M_PI/8 = " << check_norm2 / M_PI
335 dldt = killing(1).dsdt() ;
338 dldp = killing(1).stdsdp() ;
341 xidl = killing(2) * dldt + killing(3) * dldp ;
344 double xidl_error1 = 0. ;
345 double xidl_error2 = 0. ;
346 double xidl_norm = 0. ;
348 for (
int j=0; j<nt; j++) {
349 for (
int k=0; k<np/2; k++) {
354 for (
int j=0; j<nt; j++) {
355 for (
int k=np/2; k<np; k++) {
360 for (
int j=0; j<nt; j++) {
361 for (
int k=0; k<np; k++) {
367 cout <<
"Relative error in the 1st condition : "
368 << xidl_error1 / xidl_norm / double(nt) / double(np) * 2.
370 << xidl_error2 / xidl_norm / double(nt) / double(np) * 2.
372 << (xidl_error1+xidl_error2)/xidl_norm/
double(nt)/double(np)
381 dxitstdt = xitst.
dsdt() ;
391 dxi = dxitstdt + st * dxipdp ;
394 double dxi_error1 = 0. ;
395 double dxi_error2 = 0. ;
396 double dxi_norm = 0. ;
398 for (
int j=0; j<nt; j++) {
399 for (
int k=0; k<np/2; k++) {
404 for (
int j=0; j<nt; j++) {
405 for (
int k=np/2; k<np; k++) {
410 for (
int j=0; j<nt; j++) {
411 for (
int k=0; k<np; k++) {
416 cout <<
"Relative error in the 2nd condition : "
417 << dxi_error1 / dxi_norm / double(nt) / double(np) * 2.
419 << dxi_error2 / dxi_norm / double(nt) / double(np) * 2.
421 << (dxi_error1+dxi_error2)/dxi_norm/
double(nt)/double(np)
430 dxipstdt = xipst.
dsdt() ;
440 dxi2l = dxipstdt - st * dxitdp
444 double dxi2l_error1 = 0. ;
445 double dxi2l_error2 = 0. ;
446 double dxi2l_norm = 0. ;
448 for (
int j=0; j<nt; j++) {
449 for (
int k=0; k<np/2; k++) {
454 for (
int j=0; j<nt; j++) {
455 for (
int k=np/2; k<np; k++) {
460 for (
int j=0; j<nt; j++) {
461 for (
int k=0; k<np; k++) {
466 cout <<
"Relative error in the 3rd condition : "
467 << dxi2l_error1 / dxi2l_norm / double(nt) / double(np) * 2.
469 << dxi2l_error2 / dxi2l_norm / double(nt) / double(np) * 2.
471 << (dxi2l_error1+dxi2l_error2)/dxi2l_norm/
double(nt)/double(np)