LORENE
hole_bhns_extr_curv.C
1/*
2 * Method of class Hole_bhns to compute the extrinsic curvature tensor
3 *
4 * (see file hole_bhns.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2005-2007 Keisuke Taniguchi
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License version 2
15 * as published by the Free Software Foundation.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28char hole_bhns_extr_curv_C[] = "$Header: /cvsroot/Lorene/C++/Source/Hole_bhns/hole_bhns_extr_curv.C,v 1.4 2014/10/13 08:53:00 j_novak Exp $" ;
29
30/*
31 * $Id: hole_bhns_extr_curv.C,v 1.4 2014/10/13 08:53:00 j_novak Exp $
32 * $Log: hole_bhns_extr_curv.C,v $
33 * Revision 1.4 2014/10/13 08:53:00 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.3 2014/10/06 15:13:10 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.2 2008/05/15 19:05:49 k_taniguchi
40 * Change of some parameters.
41 *
42 * Revision 1.1 2007/06/22 01:24:56 k_taniguchi
43 * *** empty log message ***
44 *
45 *
46 * $Header: /cvsroot/Lorene/C++/Source/Hole_bhns/hole_bhns_extr_curv.C,v 1.4 2014/10/13 08:53:00 j_novak Exp $
47 *
48 */
49
50// C++ headers
51//#include <>
52
53// C headers
54#include <cmath>
55
56// Lorene headers
57#include "hole_bhns.h"
58#include "utilitaires.h"
59#include "unites.h"
60//#include "graphique.h"
61
62namespace Lorene {
63void Hole_bhns::extr_curv_bhns(double omega_orb, double x_rot, double y_rot) {
64
65 //----------------------------------
66 // Total extrinsic curvature tensor
67 //----------------------------------
68
69 // Fundamental constants and units
70 // -------------------------------
71 using namespace Unites ;
72
73 // Coordinates
74 // -----------
75
76 double mass = ggrav * mass_bh ;
77
78 Scalar rr(mp) ;
79 rr = mp.r ;
81 Scalar st(mp) ;
82 st = mp.sint ;
84 Scalar ct(mp) ;
85 ct = mp.cost ;
87 Scalar sp(mp) ;
88 sp = mp.sinp ;
90 Scalar cp(mp) ;
91 cp = mp.cosp ;
93
94 Vector ll(mp, CON, mp.get_bvect_cart()) ;
95 ll.set_etat_qcq() ;
96 ll.set(1) = st % cp ;
97 ll.set(2) = st % sp ;
98 ll.set(3) = ct ;
100
101 Scalar divshift(mp) ;
102 divshift = shift_auto_rs(1).deriv(1) + shift_auto_rs(2).deriv(2)
103 + shift_auto_rs(3).deriv(3) + d_shift_comp(1,1)
104 + d_shift_comp(2,2) + d_shift_comp(3,3) ;
105 divshift.std_spectral_base() ;
106
107 if (kerrschild) {
108
109 Scalar orb_rot_x(mp) ;
110 orb_rot_x = omega_orb * (mp.get_ori_x() - x_rot) ;
111 orb_rot_x.std_spectral_base() ;
112
113 Scalar orb_rot_y(mp) ;
114 orb_rot_y = omega_orb * (mp.get_ori_y() - y_rot) ;
115 orb_rot_y.std_spectral_base() ;
116
117 Vector uv(mp, CON, mp.get_bvect_cart()) ; // unit vector
118 uv.set_etat_qcq() ;
119 uv.set(1) = - orb_rot_y ;
120 uv.set(2) = orb_rot_x ;
121 uv.set(3) = 0. ;
122 uv.std_spectral_base() ;
123
124 // Computation of \tilde{A}^{ij}
125 // -----------------------------
126
127 Sym_tensor flat_taij(mp, CON, mp.get_bvect_cart()) ;
128 flat_taij.set_etat_qcq() ;
129
130 for (int i=1; i<=3; i++) {
131 for (int j=1; j<=3; j++) {
132 flat_taij.set(i,j) = shift_auto_rs(j).deriv(i)
133 + shift_auto_rs(i).deriv(j) + d_shift_comp(i,j)
134 + d_shift_comp(j,i)
135 - 2. * divshift * flat.con()(i,j) / 3. ;
136 }
137 }
138
139 flat_taij.std_spectral_base() ;
140
141 Sym_tensor curv_taij(mp, CON, mp.get_bvect_cart()) ;
142 curv_taij.set_etat_qcq() ;
143
144 for (int i=1; i<=3; i++) {
145 for (int j=1; j<=3; j++) {
146 curv_taij.set(i,j) =
147 -2. * lapconf_auto_bh * lapconf_auto_bh * mass
148 * (ll(i) * (shift_auto_rs(j).dsdr()
149 + ll(1)*d_shift_comp(1,j)
150 + ll(2)*d_shift_comp(2,j)
151 + ll(3)*d_shift_comp(3,j))
152 + ll(j) * (shift_auto_rs(i).dsdr()
153 + ll(1)*d_shift_comp(1,i)
154 + ll(2)*d_shift_comp(2,i)
155 + ll(3)*d_shift_comp(3,i))
156 - 2. * ll(i) * ll(j) * divshift / 3.) / rr ;
157 }
158 }
159
160 curv_taij.std_spectral_base() ;
161
162 Sym_tensor resi_taij(mp, CON, mp.get_bvect_cart()) ;
163 resi_taij.set_etat_qcq() ;
164
165 for (int i=1; i<=3; i++) {
166 for (int j=1; j<=3; j++) {
167 resi_taij.set(i,j) =
168 2. * lapconf_auto_bh * lapconf_auto_bh * mass
169 * ( ll(i) * (shift_auto_rs(j) + shift_comp(j))
170 + ll(j) * (shift_auto_rs(i) + shift_comp(i))
171 + ( flat.con()(i,j)
173 *(9.+14.*mass/rr)*ll(i)*ll(j) )
174 * ( ll(1) * (shift_auto_rs(1) + shift_comp(1))
175 + ll(2) * (shift_auto_rs(2) + shift_comp(2))
176 + ll(3) * (shift_auto_rs(3) + shift_comp(3)) ) / 3. )
177 / rr / rr ;
178 }
179 }
180
181 resi_taij.std_spectral_base() ;
182 resi_taij.inc_dzpuis(2) ;
183
184 taij_tot_rs = 0.5 * pow(confo_tot, 7.)
185 * (flat_taij + curv_taij + resi_taij) / lapconf_tot ;
186
189
191
192 for (int i=1; i<=3; i++) {
193 for (int j=1; j<=3; j++) {
194 taij_tot_rot.set(i,j) = pow(confo_tot,7.)
196 * ( ll(i) * uv(j) + ll(j) * uv(i)
197 + ( flat.con()(i,j)
199 *(9.+14.*mass/rr)*ll(i)*ll(j) )
200 * (ll(2)*orb_rot_x - ll(1)*orb_rot_y) / 3. )
201 / lapconf_tot / rr / rr ;
202 }
203 }
204
208
210
211 for (int i=1; i<=3; i++) {
212 for (int j=1; j<=3; j++) {
213 taij_tot_bh.set(i,j) = 2. * pow(confo_tot,7.)
214 * pow(lapconf_auto_bh,6.) * mass * (2.+3.*mass/rr)
215 * ( (1.+2.*mass/rr)*flat.con()(i,j)
216 - (3.+2.*mass/rr) * ll(i) * ll(j) )
217 / 3. / lapconf_tot / rr / rr ;
218 }
219 }
220
224
226
229
230 // Computation of \tilde{A}^{ij} \tilde{A}_{ij}
231 // --------------------------------------------
232
233 Sym_tensor flat_dshift(mp, COV, mp.get_bvect_cart()) ;
234 flat_dshift.set_etat_qcq() ;
235
236 for (int i=1; i<=3; i++) {
237 for (int j=1; j<=3; j++) {
238 flat_dshift.set(i,j) =
239 flat.cov()(j,1)*(shift_auto_rs(1).deriv(i)
240 + d_shift_comp(i,1))
241 + flat.cov()(j,2)*(shift_auto_rs(2).deriv(i)
242 + d_shift_comp(i,2))
243 + flat.cov()(j,3)*(shift_auto_rs(3).deriv(i)
244 + d_shift_comp(i,3))
245 + flat.cov()(i,1)*(shift_auto_rs(1).deriv(j)
246 + d_shift_comp(j,1))
247 + flat.cov()(i,2)*(shift_auto_rs(2).deriv(j)
248 + d_shift_comp(j,2))
249 + flat.cov()(i,3)*(shift_auto_rs(3).deriv(j)
250 + d_shift_comp(j,3))
251 - 2. * divshift * flat.cov()(i,j) / 3. ;
252 }
253 }
254
255 flat_dshift.std_spectral_base() ;
256
257 Sym_tensor curv_dshift(mp, COV, mp.get_bvect_cart()) ;
258 curv_dshift.set_etat_qcq() ;
259
260 for (int i=1; i<=3; i++) {
261 for (int j=1; j<=3; j++) {
262 curv_dshift.set(i,j) = 2. * mass
263 *( ll(j) * ( ll(1)*(shift_auto_rs(1).deriv(i)
264 + d_shift_comp(i,1))
265 + ll(2)*(shift_auto_rs(2).deriv(i)
266 + d_shift_comp(i,2))
267 + ll(3)*(shift_auto_rs(3).deriv(i)
268 + d_shift_comp(i,3)))
269 + ll(i) * ( ll(1)*(shift_auto_rs(1).deriv(j)
270 + d_shift_comp(j,1))
271 + ll(2)*(shift_auto_rs(2).deriv(j)
272 + d_shift_comp(j,2))
273 + ll(3)*(shift_auto_rs(3).deriv(j)
274 + d_shift_comp(j,3)))
275 - 2. * divshift * ll(i) * ll(j) / 3. ) / rr ;
276 }
277 }
278
279 curv_dshift.std_spectral_base() ;
280
281 Sym_tensor tmp1(mp, COV, mp.get_bvect_cart()) ;
282 tmp1.set_etat_qcq() ;
283
284 for (int i=1; i<=3; i++) {
285 for (int j=1; j<=3; j++) {
286 tmp1.set(i,j) = 2. * mass
287 *(ll(j)*( (flat.cov()(i,1)+2.*mass*ll(i)*ll(1)/rr)
288 * (shift_auto_rs(1) + shift_comp(1))
289 + (flat.cov()(i,2)+2.*mass*ll(i)*ll(2)/rr)
290 * (shift_auto_rs(2) + shift_comp(2))
291 + (flat.cov()(i,3)+2.*mass*ll(i)*ll(3)/rr)
292 * (shift_auto_rs(3) + shift_comp(3))
293 )
294 + ll(i)*( (flat.cov()(j,1)+2.*mass*ll(j)*ll(1)/rr)
295 * (shift_auto_rs(1) + shift_comp(1))
296 + (flat.cov()(j,2)+2.*mass*ll(j)*ll(2)/rr)
297 * (shift_auto_rs(2) + shift_comp(2))
298 + (flat.cov()(j,3)+2.*mass*ll(j)*ll(3)/rr)
299 * (shift_auto_rs(3) + shift_comp(3)) )
300 ) / rr / rr ;
301 }
302 }
303 tmp1.std_spectral_base() ;
304 tmp1.inc_dzpuis(2) ;
305
306 Sym_tensor tmp2(mp, COV, mp.get_bvect_cart()) ;
307 tmp2.set_etat_qcq() ;
308
309 for (int i=1; i<=3; i++) {
310 for (int j=1; j<=3; j++) {
311 tmp2.set(i,j) = 2. * mass * lapconf_auto_bh * lapconf_auto_bh
312 * ( ll(1) * (shift_auto_rs(1) + shift_comp(1))
313 + ll(2) * (shift_auto_rs(2) + shift_comp(2))
314 + ll(3) * (shift_auto_rs(3) + shift_comp(3)) )
315 * (flat.cov()(i,j)
316 - (9.+28.*mass/rr+24.*mass*mass/rr/rr)*ll(i)*ll(j))
317 / 3. / rr / rr ;
318 }
319 }
320 tmp2.std_spectral_base() ;
321 tmp2.inc_dzpuis(2) ;
322
323 Sym_tensor taij_down_rs(mp, COV, mp.get_bvect_cart()) ;
324 taij_down_rs.set_etat_qcq() ;
325
326 taij_down_rs = 0.5 * pow(confo_tot,7.)
327 * (flat_dshift + curv_dshift + tmp1 + tmp2) / lapconf_tot ;
328
329 taij_down_rs.std_spectral_base() ;
330 taij_down_rs.annule_domain(0) ;
331
332 Sym_tensor taij_down_rot(mp, COV, mp.get_bvect_cart()) ;
333 taij_down_rot.set_etat_qcq() ;
334
335 for (int i=1; i<=3; i++) {
336 for (int j=1; j<=3; j++) {
337 taij_down_rot.set(i,j) = mass * pow(confo_tot,7.)
338 * ( ll(j)*uv(i) + ll(i)*uv(j)
340 * (ll(2)*orb_rot_x - ll(1)*orb_rot_y)
341 * ( flat.cov()(i,j) - (9.+16.*mass/rr)*ll(i)*ll(j) ) / 3.
342 ) / lapconf_tot / rr / rr ;
343 }
344 }
345 taij_down_rot.std_spectral_base() ;
346 taij_down_rot.annule_domain(0) ;
347 taij_down_rot.inc_dzpuis(2) ;
348
349 Sym_tensor taij_down_bh(mp, COV, mp.get_bvect_cart()) ;
350 taij_down_bh.set_etat_qcq() ;
351
352 for (int i=1; i<=3; i++) {
353 for (int j=1; j<=3; j++) {
354 taij_down_bh.set(i,j) = 2. * mass * pow(confo_tot,7.)
355 * pow(lapconf_auto_bh,4.) * (2.+3.*mass/rr)
356 * (flat.cov()(i,j) - (3.+4.*mass/rr) * ll(i) * ll(j))
357 / 3. / lapconf_auto / rr / rr ;
358 }
359 }
360 taij_down_bh.std_spectral_base() ;
361 taij_down_bh.annule_domain(0) ;
362 taij_down_bh.inc_dzpuis(2) ;
363
364 Scalar taij_quad_rstot(mp) ;
365 taij_quad_rstot = 0. ;
366
367 for (int i=1; i<=3; i++) {
368 for (int j=1; j<=3; j++) {
369 taij_quad_rstot += taij_down_rs(i,j) * taij_tot(i,j) ;
370 }
371 }
372 taij_quad_rstot.std_spectral_base() ;
373
374 Scalar taij_quad_rsrotbh(mp) ;
375 taij_quad_rsrotbh = 0. ;
376
377 for (int i=1; i<=3; i++) {
378 for (int j=1; j<=3; j++) {
379 taij_quad_rsrotbh += taij_tot_rs(i,j)
380 * (taij_down_rot(i,j) + taij_down_bh(i,j)) ;
381 }
382 }
383 taij_quad_rsrotbh.std_spectral_base() ;
384
385 taij_quad_tot_rs = taij_quad_rstot + taij_quad_rsrotbh ;
387
388 taij_quad_tot_rot = 8.*mass*mass*pow(confo_tot,14.)
389 * pow(lapconf_auto_bh,6.) * (2.+3.*mass/rr)
390 * (ll(2)*orb_rot_x - ll(1)*orb_rot_y)
391 / 3. / lapconf_tot / lapconf_tot / pow(rr,4.)
392 + 2.*mass*mass*pow(confo_tot,14.)*pow(lapconf_auto_bh,4.)
393 * (3.*(1.+2.*mass/rr)*(orb_rot_x*orb_rot_x+orb_rot_y*orb_rot_y)
394 -2.*(1.+3.*mass/rr)*(ll(2)*orb_rot_x-ll(1)*orb_rot_y)
395 *(ll(2)*orb_rot_x-ll(1)*orb_rot_y))
396 / 3. / lapconf_tot / lapconf_tot / pow(rr,4.) ;
397
400
401 taij_quad_tot_bh = 8.*mass*mass*pow(confo_tot,14.)
402 * pow(lapconf_auto_bh,8.) * (2.+3.*mass/rr) * (2.+3.*mass/rr)
403 / 3. / lapconf_tot / lapconf_tot / pow(rr,4.) ;
404
407
411
412 }
413 else { // Isotropic coordinates with the maximal slicing
414
415 // Sets C/M^2 for each case of the lapse boundary condition
416 // --------------------------------------------------------
417 double cc ;
418
419 if (bc_lapconf_nd) { // Neumann boundary condition
420 if (bc_lapconf_fs) { // First condition
421 // d(\alpha \psi)/dr = 0
422 // ---------------------
423 cc = 2. * (sqrt(13.) - 1.) / 3. ;
424 }
425 else { // Second condition
426 // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
427 // -----------------------------------------
428 cc = 4. / 3. ;
429 }
430 }
431 else { // Dirichlet boundary condition
432 if (bc_lapconf_fs) { // First condition
433 // (\alpha \psi) = 1/2
434 // -------------------
435 cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
436 abort() ;
437 }
438 else { // Second condition
439 // (\alpha \psi) = 1/sqrt(2.) \psi_KS
440 // ----------------------------------
441 cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
442 abort() ;
443 // cc = 2. * sqrt(2.) ;
444 }
445 }
446
447 Scalar r_are(mp) ;
449 r_are.std_spectral_base() ;
450
451 // Computation of \tilde{A}^{ij}
452 // -----------------------------
453
454 Sym_tensor flat_taij(mp, CON, mp.get_bvect_cart()) ;
455 flat_taij.set_etat_qcq() ;
456
457 for (int i=1; i<=3; i++) {
458 for (int j=1; j<=3; j++) {
459 flat_taij.set(i,j) = shift_auto_rs(j).deriv(i)
460 + shift_auto_rs(i).deriv(j) + d_shift_comp(i,j)
461 + d_shift_comp(j,i)
462 - 2. * divshift % flat.con()(i,j) / 3. ;
463 }
464 }
465
466 flat_taij.std_spectral_base() ;
467
468 taij_tot_rs = 0.5 * pow(confo_tot, 7.) * flat_taij / lapconf_tot ;
469
472
473 if (taij_tot_rs(1,2).get_etat() == ETATQCQ) {
474 for (int i=1; i<=3; i++) {
475 for (int j=1; j<=3; j++) {
476 taij_tot_rs.set(i,j).raccord(1) ;
477 }
478 }
479 }
480
482
483 for (int i=1; i<=3; i++) {
484 for (int j=1; j<=3; j++) {
485 taij_tot_bh.set(i,j) = pow(confo_tot,7.)*mass*mass*cc
486 * sqrt(1. - 2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
487 * (flat.con()(i,j) - 3.*ll(i)*ll(j)) / lapconf_tot
488 / pow(r_are*rr,3.) ;
489 }
490 }
491
494
495 for (int i=1; i<=3; i++) {
496 for (int j=1; j<=3; j++) {
497 taij_tot_bh.set(i,j).raccord(1) ;
498 }
499 }
500
502
504
507
508 for (int i=1; i<=3; i++) {
509 for (int j=1; j<=3; j++) {
510 taij_tot.set(i,j).raccord(1) ;
511 }
512 }
513
514 for (int i=1; i<=3; i++) {
515 for (int j=1; j<=3; j++) {
516 taij_tot_rot.set(i,j) = 0. ;
517 }
518 }
520
521 // Computation of \tilde{A}_{BH}^{ij}
522 // ----------------------------------
523
524 Scalar divshift_auto(mp) ;
525 divshift_auto = shift_auto_rs(1).deriv(1)
526 + shift_auto_rs(2).deriv(2) + shift_auto_rs(3).deriv(3) ;
527 divshift_auto.std_spectral_base() ;
528
529 Sym_tensor flat_taij_auto_rs(mp, CON, mp.get_bvect_cart()) ;
530 flat_taij_auto_rs.set_etat_qcq() ;
531
532 for (int i=1; i<=3; i++) {
533 for (int j=1; j<=3; j++) {
534 flat_taij_auto_rs.set(i,j) = shift_auto_rs(j).deriv(i)
535 + shift_auto_rs(i).deriv(j)
536 - 2. * divshift_auto % flat.con()(i,j) / 3. ;
537 }
538 }
539
540 flat_taij_auto_rs.std_spectral_base() ;
541
542 taij_auto_rs = 0.5 * pow(confo_tot, 7.) * flat_taij_auto_rs
543 / lapconf_tot ;
544
547
548 if (taij_auto_rs(1,2).get_etat() == ETATQCQ) {
549 for (int i=1; i<=3; i++) {
550 for (int j=1; j<=3; j++) {
551 taij_auto_rs.set(i,j).raccord(1) ;
552 }
553 }
554 }
555
557
560
561 for (int i=1; i<=3; i++) {
562 for (int j=1; j<=3; j++) {
563 taij_auto.set(i,j).raccord(1) ;
564 }
565 }
566
567 // Computation of \tilde{A}_{NS}^{ij}
568 // ----------------------------------
569
570 Scalar divshift_comp(mp) ;
571 divshift_comp = d_shift_comp(1,1) + d_shift_comp(2,2)
572 + d_shift_comp(3,3) ;
573 divshift_comp.std_spectral_base() ;
574
575 Sym_tensor flat_taij_comp(mp, CON, mp.get_bvect_cart()) ;
576 flat_taij_comp.set_etat_qcq() ;
577
578 for (int i=1; i<=3; i++) {
579 for (int j=1; j<=3; j++) {
580 flat_taij_comp.set(i,j) = d_shift_comp(i,j)
581 + d_shift_comp(j,i)
582 - 2. * divshift_comp % flat.con()(i,j) / 3. ;
583 }
584 }
585
586 flat_taij_comp.std_spectral_base() ;
587
588 taij_comp = 0.5 * pow(confo_comp+0.5, 7.) * flat_taij_comp
589 / (lapconf_comp+0.5) ;
590
593
594 if (taij_comp(1,2).get_etat() == ETATQCQ) {
595 for (int i=1; i<=3; i++) {
596 for (int j=1; j<=3; j++) {
597 taij_comp.set(i,j).raccord(1) ;
598 }
599 }
600 }
601
602 // Computation of \tilde{A}^{ij} \tilde{A}_{ij}
603 // --------------------------------------------
604
605 Sym_tensor flat_dshift(mp, COV, mp.get_bvect_cart()) ;
606 flat_dshift.set_etat_qcq() ;
607
608 for (int i=1; i<=3; i++) {
609 for (int j=1; j<=3; j++) {
610 flat_dshift.set(i,j) =
611 flat.cov()(j,1)*(shift_auto_rs(1).deriv(i)
612 + d_shift_comp(i,1))
613 + flat.cov()(j,2)*(shift_auto_rs(2).deriv(i)
614 + d_shift_comp(i,2))
615 + flat.cov()(j,3)*(shift_auto_rs(3).deriv(i)
616 + d_shift_comp(i,3))
617 + flat.cov()(i,1)*(shift_auto_rs(1).deriv(j)
618 + d_shift_comp(j,1))
619 + flat.cov()(i,2)*(shift_auto_rs(2).deriv(j)
620 + d_shift_comp(j,2))
621 + flat.cov()(i,3)*(shift_auto_rs(3).deriv(j)
622 + d_shift_comp(j,3))
623 - 2. * divshift * flat.cov()(i,j) / 3. ;
624 }
625 }
626
627 Sym_tensor taij_down_rs(mp, COV, mp.get_bvect_cart()) ;
628 taij_down_rs.set_etat_qcq() ;
629
630 taij_down_rs = 0.5 * pow(confo_tot, 7.) * flat_dshift / lapconf_tot ;
631
632 taij_down_rs.std_spectral_base() ;
633 taij_down_rs.annule_domain(0) ;
634
635 Sym_tensor taij_down_bh(mp, COV, mp.get_bvect_cart()) ;
636 taij_down_bh.set_etat_qcq() ;
637
638 for (int i=1; i<=3; i++) {
639 for (int j=1; j<=3; j++) {
640 taij_down_bh.set(i,j) = pow(confo_tot,7.)*mass*mass*cc
641 * sqrt(1. - 2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
642 * (flat.cov()(i,j) - 3.*ll(i)%ll(j)) / lapconf_tot
643 / pow(r_are*rr,3.) ;
644 }
645 }
646 taij_down_bh.std_spectral_base() ;
647 taij_down_bh.annule_domain(0) ;
648
649 for (int i=1; i<=3; i++) {
650 for (int j=1; j<=3; j++) {
651 taij_down_bh.set(i,j).raccord(1) ;
652 }
653 }
654
655 taij_down_bh.inc_dzpuis(2) ;
656
657 Scalar taij_quad_rstot(mp) ;
658 taij_quad_rstot = 0. ;
659
660 for (int i=1; i<=3; i++) {
661 for (int j=1; j<=3; j++) {
662 taij_quad_rstot += taij_down_rs(i,j) % taij_tot(i,j) ;
663 }
664 }
665 taij_quad_rstot.std_spectral_base() ;
666
667 Scalar taij_quad_rsbh(mp) ;
668 taij_quad_rsbh = 0. ;
669
670 for (int i=1; i<=3; i++) {
671 for (int j=1; j<=3; j++) {
672 taij_quad_rsbh += taij_tot_rs(i,j) % taij_down_bh(i,j) ;
673 }
674 }
675 taij_quad_rsbh.std_spectral_base() ;
676
677 taij_quad_tot_rs = taij_quad_rstot + taij_quad_rsbh ;
679
680 taij_quad_tot_rot = 0. ;
682
683 taij_quad_tot_bh = 6.*pow(confo_tot,14.)*pow(mass*mass*cc,2.)
684 * (1. - 2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
685 / lapconf_tot / lapconf_tot / pow(r_are*rr, 6.) ;
689
691
696
697 // -------------------------
698 Scalar taij_quad_auto1(mp) ;
699 taij_quad_auto1 = 0. ;
700 for (int i=1; i<=3; i++) {
701 for (int j=1; j<=3; j++) {
702 taij_quad_auto1 += taij_auto_rs(i,j)
703 * (taij_down_rs(i,j)
704 + pow(confo_tot/(confo_comp+0.5),7.)*(lapconf_comp+0.5)
705 * taij_comp(i,j) / lapconf_tot) ;
706 }
707 }
708 taij_quad_auto1.std_spectral_base() ;
709
710 Scalar taij_quad_auto2(mp) ;
711 taij_quad_auto2 = 0. ;
712 for (int i=1; i<=3; i++) {
713 for (int j=1; j<=3; j++) {
714 taij_quad_auto2 += taij_tot_bh(i,j) % taij_down_rs(i,j) ;
715 }
716 }
717 taij_quad_auto2.std_spectral_base() ;
718
719 taij_quad_auto = taij_quad_auto1 + 2.*taij_quad_auto2 ;
722 if (taij_quad_auto.get_etat() == ETATQCQ) {
724 }
725
726 // Computation of \tilde{A}_{NS}^{ij} \tilde{A}^{NS}_{ij}
727 // ------------------------------------------------------
728
729 taij_quad_comp = 0. ;
730 for (int i=1; i<=3; i++) {
731 for (int j=1; j<=3; j++) {
732 taij_quad_comp += taij_comp(i,j) % taij_comp(i,j) ;
733 }
734 }
736
737 }
738
739 // The derived quantities are obsolete
740 // -----------------------------------
741
742 del_deriv() ;
743
744}
745}
const Scalar r_coord(bool neumann, bool first) const
Expresses the areal radial coordinate by that in spatially isotropic coordinates.
Map & mp
Mapping associated with the black hole.
Definition blackhole.h:80
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the black hol...
Definition blackhole.h:143
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
Definition blackhole.h:85
double mass_bh
Gravitational mass of BH.
Definition blackhole.h:88
Sym_tensor taij_tot_bh
Part of the extrinsic curvature tensor from the analytic background.
Definition hole_bhns.h:200
Scalar lapconf_auto
Lapconf function generated by the black hole.
Definition hole_bhns.h:95
Sym_tensor taij_comp
Part of the extrinsic curvature tensor generated by the companion star.
Definition hole_bhns.h:221
Sym_tensor taij_auto
Part of the extrinsic curvature tensor generated by the black hole.
Definition hole_bhns.h:216
Sym_tensor taij_tot_rs
Part of the extrinsic curvature tensor from the numerical computation.
Definition hole_bhns.h:190
Scalar taij_quad_auto
Part of the scalar from the black hole.
Definition hole_bhns.h:238
Vector shift_auto_rs
Part of the shift vector from the numerical computation.
Definition hole_bhns.h:126
virtual void del_deriv() const
Deletes all the derived quantities.
Definition hole_bhns.C:392
bool bc_lapconf_fs
true for the first type BC for the lapconf function, false for the second type BH
Definition hole_bhns.h:78
Scalar taij_quad_tot_rot
Part of the scalar from the rotation shift vector.
Definition hole_bhns.h:227
Scalar taij_quad_tot
Total scalar generated by .
Definition hole_bhns.h:235
Sym_tensor taij_tot_rot
Part of the extrinsic curvature tensor from the rotation shift vector.
Definition hole_bhns.h:195
Sym_tensor taij_tot
Total extrinsic curvature tensor generated by shift_tot , lapse_tot , and confo_tot .
Definition hole_bhns.h:206
Scalar confo_comp
Conformal factor generated by the companion star.
Definition hole_bhns.h:166
Scalar lapconf_comp
Lapconf function generated by the companion star.
Definition hole_bhns.h:98
Scalar lapconf_auto_bh
Part of the lapconf function from the analytic background.
Definition hole_bhns.h:92
Scalar taij_quad_comp
Part of the scalar from the companion star.
Definition hole_bhns.h:241
void extr_curv_bhns(double omega_orb, double x_rot, double y_rot)
Computes taij_tot , taij_quad_tot from shift_tot , lapse_tot , confo_tot .
bool bc_lapconf_nd
true for the Neumann type BC for the lapconf function, false for the Dirichlet type BH
Definition hole_bhns.h:73
Scalar lapconf_tot
Total lapconf function.
Definition hole_bhns.h:101
Tensor d_shift_comp
Derivative of the shift vector generated by the companion star.
Definition hole_bhns.h:154
Vector shift_comp
Shift vector generated by the companion star.
Definition hole_bhns.h:135
Scalar confo_tot
Total conformal factor.
Definition hole_bhns.h:169
Sym_tensor taij_auto_rs
Part of the extrinsic curvature tensor numericalty computed for the black hole.
Definition hole_bhns.h:211
Scalar taij_quad_tot_rs
Part of the scalar from the numerical computation.
Definition hole_bhns.h:224
Scalar taij_quad_tot_bh
Part of the scalar from the analytic background.
Definition hole_bhns.h:230
Coord cosp
Definition map.h:724
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition map.h:791
Coord sint
Definition map.h:721
Coord r
r coordinate centered on the grid
Definition map.h:718
double get_ori_y() const
Returns the y coordinate of the origin.
Definition map.h:770
Coord sinp
Definition map.h:723
double get_ori_x() const
Returns the x coordinate of the origin.
Definition map.h:768
Coord cost
Definition map.h:722
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition scalar.C:784
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition scalar.h:554
Class intended to describe valence-2 symmetric tensors.
Definition sym_tensor.h:223
Tensor field of valence 1.
Definition vector.h:188
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:316
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition tensor.C:816
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition tensor.C:666
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition tensor.C:481
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition tensor.C:926
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.