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 
28 char 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 
62 namespace Lorene {
63 void 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 ;
80  rr.std_spectral_base() ;
81  Scalar st(mp) ;
82  st = mp.sint ;
83  st.std_spectral_base() ;
84  Scalar ct(mp) ;
85  ct = mp.cost ;
86  ct.std_spectral_base() ;
87  Scalar sp(mp) ;
88  sp = mp.sinp ;
89  sp.std_spectral_base() ;
90  Scalar cp(mp) ;
91  cp = mp.cosp ;
92  cp.std_spectral_base() ;
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 ;
99  ll.std_spectral_base() ;
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 
409  + taij_quad_tot_bh ;
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 }
Map & mp
Mapping associated with the black hole.
Definition: blackhole.h:80
Sym_tensor taij_auto
Part of the extrinsic curvature tensor generated by the black hole.
Definition: hole_bhns.h:216
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:666
Scalar lapconf_auto_bh
Part of the lapconf function from the analytic background.
Definition: hole_bhns.h:92
virtual void set_etat_qcq()
Sets the logical state of all components to ETATQCQ (ordinary state).
Definition: tensor.C:481
Scalar taij_quad_comp
Part of the scalar from the companion star.
Definition: hole_bhns.h:241
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Sym_tensor taij_tot_rot
Part of the extrinsic curvature tensor from the rotation shift vector.
Definition: hole_bhns.h:195
double get_ori_y() const
Returns the y coordinate of the origin.
Definition: map.h:770
Lorene prototypes.
Definition: app_hor.h:64
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
Standard units of space, time and mass.
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
bool kerrschild
true for a Kerr-Schild background, false for a conformally flat background
Definition: blackhole.h:85
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:768
Sym_tensor taij_comp
Part of the extrinsic curvature tensor generated by the companion star.
Definition: hole_bhns.h:221
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
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:554
Tensor d_shift_comp
Derivative of the shift vector generated by the companion star.
Definition: hole_bhns.h:154
Tensor field of valence 1.
Definition: vector.h:188
Sym_tensor taij_auto_rs
Part of the extrinsic curvature tensor numericalty computed for the black hole.
Definition: hole_bhns.h:211
Scalar lapconf_comp
Lapconf function generated by the companion star.
Definition: hole_bhns.h:98
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
Coord sint
Definition: map.h:721
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric_flat.C:153
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...
Scalar confo_comp
Conformal factor generated by the companion star.
Definition: hole_bhns.h:166
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
Coord sinp
Definition: map.h:723
Vector shift_comp
Shift vector generated by the companion star.
Definition: hole_bhns.h:135
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Sym_tensor taij_tot
Total extrinsic curvature tensor generated by shift_tot , lapse_tot , and confo_tot ...
Definition: hole_bhns.h:206
Scalar lapconf_auto
Lapconf function generated by the black hole.
Definition: hole_bhns.h:95
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
Scalar confo_tot
Total conformal factor.
Definition: hole_bhns.h:169
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric_flat.C:134
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:721
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:791
Scalar taij_quad_tot
Total scalar generated by .
Definition: hole_bhns.h:235
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:654
Metric_flat flat
Flat metric defined on the mapping (Spherical components with respect to the mapping of the black hol...
Definition: blackhole.h:143
Sym_tensor taij_tot_rs
Part of the extrinsic curvature tensor from the numerical computation.
Definition: hole_bhns.h:190
Coord cosp
Definition: map.h:724
Scalar taij_quad_tot_rot
Part of the scalar from the rotation shift vector.
Definition: hole_bhns.h:227
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
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
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 .
const Scalar r_coord(bool neumann, bool first) const
Expresses the areal radial coordinate by that in spatially isotropic coordinates. ...
Scalar taij_quad_tot_bh
Part of the scalar from the analytic background.
Definition: hole_bhns.h:230
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: tensor.C:926
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:223
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: hole_bhns.C:392
Scalar lapconf_tot
Total lapconf function.
Definition: hole_bhns.h:101
Coord r
r coordinate centered on the grid
Definition: map.h:718
Scalar taij_quad_tot_rs
Part of the scalar from the numerical computation.
Definition: hole_bhns.h:224
Coord cost
Definition: map.h:722