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 
29 
30 /*
31  * $Id: hole_bhns_extr_curv.C,v 1.5 2016/12/05 16:17:55 j_novak Exp $
32  * $Log: hole_bhns_extr_curv.C,v $
33  * Revision 1.5 2016/12/05 16:17:55 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.4 2014/10/13 08:53:00 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.3 2014/10/06 15:13:10 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.2 2008/05/15 19:05:49 k_taniguchi
43  * Change of some parameters.
44  *
45  * Revision 1.1 2007/06/22 01:24:56 k_taniguchi
46  * *** empty log message ***
47  *
48  *
49  * $Header: /cvsroot/Lorene/C++/Source/Hole_bhns/hole_bhns_extr_curv.C,v 1.5 2016/12/05 16:17:55 j_novak Exp $
50  *
51  */
52 
53 // C++ headers
54 //#include <>
55 
56 // C headers
57 #include <cmath>
58 
59 // Lorene headers
60 #include "hole_bhns.h"
61 #include "utilitaires.h"
62 #include "unites.h"
63 //#include "graphique.h"
64 
65 namespace Lorene {
66 void Hole_bhns::extr_curv_bhns(double omega_orb, double x_rot, double y_rot) {
67 
68  //----------------------------------
69  // Total extrinsic curvature tensor
70  //----------------------------------
71 
72  // Fundamental constants and units
73  // -------------------------------
74  using namespace Unites ;
75 
76  // Coordinates
77  // -----------
78 
79  double mass = ggrav * mass_bh ;
80 
81  Scalar rr(mp) ;
82  rr = mp.r ;
83  rr.std_spectral_base() ;
84  Scalar st(mp) ;
85  st = mp.sint ;
86  st.std_spectral_base() ;
87  Scalar ct(mp) ;
88  ct = mp.cost ;
89  ct.std_spectral_base() ;
90  Scalar sp(mp) ;
91  sp = mp.sinp ;
92  sp.std_spectral_base() ;
93  Scalar cp(mp) ;
94  cp = mp.cosp ;
95  cp.std_spectral_base() ;
96 
97  Vector ll(mp, CON, mp.get_bvect_cart()) ;
98  ll.set_etat_qcq() ;
99  ll.set(1) = st % cp ;
100  ll.set(2) = st % sp ;
101  ll.set(3) = ct ;
102  ll.std_spectral_base() ;
103 
104  Scalar divshift(mp) ;
105  divshift = shift_auto_rs(1).deriv(1) + shift_auto_rs(2).deriv(2)
106  + shift_auto_rs(3).deriv(3) + d_shift_comp(1,1)
107  + d_shift_comp(2,2) + d_shift_comp(3,3) ;
108  divshift.std_spectral_base() ;
109 
110  if (kerrschild) {
111 
112  Scalar orb_rot_x(mp) ;
113  orb_rot_x = omega_orb * (mp.get_ori_x() - x_rot) ;
114  orb_rot_x.std_spectral_base() ;
115 
116  Scalar orb_rot_y(mp) ;
117  orb_rot_y = omega_orb * (mp.get_ori_y() - y_rot) ;
118  orb_rot_y.std_spectral_base() ;
119 
120  Vector uv(mp, CON, mp.get_bvect_cart()) ; // unit vector
121  uv.set_etat_qcq() ;
122  uv.set(1) = - orb_rot_y ;
123  uv.set(2) = orb_rot_x ;
124  uv.set(3) = 0. ;
125  uv.std_spectral_base() ;
126 
127  // Computation of \tilde{A}^{ij}
128  // -----------------------------
129 
130  Sym_tensor flat_taij(mp, CON, mp.get_bvect_cart()) ;
131  flat_taij.set_etat_qcq() ;
132 
133  for (int i=1; i<=3; i++) {
134  for (int j=1; j<=3; j++) {
135  flat_taij.set(i,j) = shift_auto_rs(j).deriv(i)
136  + shift_auto_rs(i).deriv(j) + d_shift_comp(i,j)
137  + d_shift_comp(j,i)
138  - 2. * divshift * flat.con()(i,j) / 3. ;
139  }
140  }
141 
142  flat_taij.std_spectral_base() ;
143 
144  Sym_tensor curv_taij(mp, CON, mp.get_bvect_cart()) ;
145  curv_taij.set_etat_qcq() ;
146 
147  for (int i=1; i<=3; i++) {
148  for (int j=1; j<=3; j++) {
149  curv_taij.set(i,j) =
150  -2. * lapconf_auto_bh * lapconf_auto_bh * mass
151  * (ll(i) * (shift_auto_rs(j).dsdr()
152  + ll(1)*d_shift_comp(1,j)
153  + ll(2)*d_shift_comp(2,j)
154  + ll(3)*d_shift_comp(3,j))
155  + ll(j) * (shift_auto_rs(i).dsdr()
156  + ll(1)*d_shift_comp(1,i)
157  + ll(2)*d_shift_comp(2,i)
158  + ll(3)*d_shift_comp(3,i))
159  - 2. * ll(i) * ll(j) * divshift / 3.) / rr ;
160  }
161  }
162 
163  curv_taij.std_spectral_base() ;
164 
165  Sym_tensor resi_taij(mp, CON, mp.get_bvect_cart()) ;
166  resi_taij.set_etat_qcq() ;
167 
168  for (int i=1; i<=3; i++) {
169  for (int j=1; j<=3; j++) {
170  resi_taij.set(i,j) =
171  2. * lapconf_auto_bh * lapconf_auto_bh * mass
172  * ( ll(i) * (shift_auto_rs(j) + shift_comp(j))
173  + ll(j) * (shift_auto_rs(i) + shift_comp(i))
174  + ( flat.con()(i,j)
176  *(9.+14.*mass/rr)*ll(i)*ll(j) )
177  * ( ll(1) * (shift_auto_rs(1) + shift_comp(1))
178  + ll(2) * (shift_auto_rs(2) + shift_comp(2))
179  + ll(3) * (shift_auto_rs(3) + shift_comp(3)) ) / 3. )
180  / rr / rr ;
181  }
182  }
183 
184  resi_taij.std_spectral_base() ;
185  resi_taij.inc_dzpuis(2) ;
186 
187  taij_tot_rs = 0.5 * pow(confo_tot, 7.)
188  * (flat_taij + curv_taij + resi_taij) / lapconf_tot ;
189 
192 
194 
195  for (int i=1; i<=3; i++) {
196  for (int j=1; j<=3; j++) {
197  taij_tot_rot.set(i,j) = pow(confo_tot,7.)
199  * ( ll(i) * uv(j) + ll(j) * uv(i)
200  + ( flat.con()(i,j)
202  *(9.+14.*mass/rr)*ll(i)*ll(j) )
203  * (ll(2)*orb_rot_x - ll(1)*orb_rot_y) / 3. )
204  / lapconf_tot / rr / rr ;
205  }
206  }
207 
211 
213 
214  for (int i=1; i<=3; i++) {
215  for (int j=1; j<=3; j++) {
216  taij_tot_bh.set(i,j) = 2. * pow(confo_tot,7.)
217  * pow(lapconf_auto_bh,6.) * mass * (2.+3.*mass/rr)
218  * ( (1.+2.*mass/rr)*flat.con()(i,j)
219  - (3.+2.*mass/rr) * ll(i) * ll(j) )
220  / 3. / lapconf_tot / rr / rr ;
221  }
222  }
223 
227 
229 
232 
233  // Computation of \tilde{A}^{ij} \tilde{A}_{ij}
234  // --------------------------------------------
235 
236  Sym_tensor flat_dshift(mp, COV, mp.get_bvect_cart()) ;
237  flat_dshift.set_etat_qcq() ;
238 
239  for (int i=1; i<=3; i++) {
240  for (int j=1; j<=3; j++) {
241  flat_dshift.set(i,j) =
242  flat.cov()(j,1)*(shift_auto_rs(1).deriv(i)
243  + d_shift_comp(i,1))
244  + flat.cov()(j,2)*(shift_auto_rs(2).deriv(i)
245  + d_shift_comp(i,2))
246  + flat.cov()(j,3)*(shift_auto_rs(3).deriv(i)
247  + d_shift_comp(i,3))
248  + flat.cov()(i,1)*(shift_auto_rs(1).deriv(j)
249  + d_shift_comp(j,1))
250  + flat.cov()(i,2)*(shift_auto_rs(2).deriv(j)
251  + d_shift_comp(j,2))
252  + flat.cov()(i,3)*(shift_auto_rs(3).deriv(j)
253  + d_shift_comp(j,3))
254  - 2. * divshift * flat.cov()(i,j) / 3. ;
255  }
256  }
257 
258  flat_dshift.std_spectral_base() ;
259 
260  Sym_tensor curv_dshift(mp, COV, mp.get_bvect_cart()) ;
261  curv_dshift.set_etat_qcq() ;
262 
263  for (int i=1; i<=3; i++) {
264  for (int j=1; j<=3; j++) {
265  curv_dshift.set(i,j) = 2. * mass
266  *( ll(j) * ( ll(1)*(shift_auto_rs(1).deriv(i)
267  + d_shift_comp(i,1))
268  + ll(2)*(shift_auto_rs(2).deriv(i)
269  + d_shift_comp(i,2))
270  + ll(3)*(shift_auto_rs(3).deriv(i)
271  + d_shift_comp(i,3)))
272  + ll(i) * ( ll(1)*(shift_auto_rs(1).deriv(j)
273  + d_shift_comp(j,1))
274  + ll(2)*(shift_auto_rs(2).deriv(j)
275  + d_shift_comp(j,2))
276  + ll(3)*(shift_auto_rs(3).deriv(j)
277  + d_shift_comp(j,3)))
278  - 2. * divshift * ll(i) * ll(j) / 3. ) / rr ;
279  }
280  }
281 
282  curv_dshift.std_spectral_base() ;
283 
284  Sym_tensor tmp1(mp, COV, mp.get_bvect_cart()) ;
285  tmp1.set_etat_qcq() ;
286 
287  for (int i=1; i<=3; i++) {
288  for (int j=1; j<=3; j++) {
289  tmp1.set(i,j) = 2. * mass
290  *(ll(j)*( (flat.cov()(i,1)+2.*mass*ll(i)*ll(1)/rr)
291  * (shift_auto_rs(1) + shift_comp(1))
292  + (flat.cov()(i,2)+2.*mass*ll(i)*ll(2)/rr)
293  * (shift_auto_rs(2) + shift_comp(2))
294  + (flat.cov()(i,3)+2.*mass*ll(i)*ll(3)/rr)
295  * (shift_auto_rs(3) + shift_comp(3))
296  )
297  + ll(i)*( (flat.cov()(j,1)+2.*mass*ll(j)*ll(1)/rr)
298  * (shift_auto_rs(1) + shift_comp(1))
299  + (flat.cov()(j,2)+2.*mass*ll(j)*ll(2)/rr)
300  * (shift_auto_rs(2) + shift_comp(2))
301  + (flat.cov()(j,3)+2.*mass*ll(j)*ll(3)/rr)
302  * (shift_auto_rs(3) + shift_comp(3)) )
303  ) / rr / rr ;
304  }
305  }
306  tmp1.std_spectral_base() ;
307  tmp1.inc_dzpuis(2) ;
308 
309  Sym_tensor tmp2(mp, COV, mp.get_bvect_cart()) ;
310  tmp2.set_etat_qcq() ;
311 
312  for (int i=1; i<=3; i++) {
313  for (int j=1; j<=3; j++) {
314  tmp2.set(i,j) = 2. * mass * lapconf_auto_bh * lapconf_auto_bh
315  * ( ll(1) * (shift_auto_rs(1) + shift_comp(1))
316  + ll(2) * (shift_auto_rs(2) + shift_comp(2))
317  + ll(3) * (shift_auto_rs(3) + shift_comp(3)) )
318  * (flat.cov()(i,j)
319  - (9.+28.*mass/rr+24.*mass*mass/rr/rr)*ll(i)*ll(j))
320  / 3. / rr / rr ;
321  }
322  }
323  tmp2.std_spectral_base() ;
324  tmp2.inc_dzpuis(2) ;
325 
326  Sym_tensor taij_down_rs(mp, COV, mp.get_bvect_cart()) ;
327  taij_down_rs.set_etat_qcq() ;
328 
329  taij_down_rs = 0.5 * pow(confo_tot,7.)
330  * (flat_dshift + curv_dshift + tmp1 + tmp2) / lapconf_tot ;
331 
332  taij_down_rs.std_spectral_base() ;
333  taij_down_rs.annule_domain(0) ;
334 
335  Sym_tensor taij_down_rot(mp, COV, mp.get_bvect_cart()) ;
336  taij_down_rot.set_etat_qcq() ;
337 
338  for (int i=1; i<=3; i++) {
339  for (int j=1; j<=3; j++) {
340  taij_down_rot.set(i,j) = mass * pow(confo_tot,7.)
341  * ( ll(j)*uv(i) + ll(i)*uv(j)
343  * (ll(2)*orb_rot_x - ll(1)*orb_rot_y)
344  * ( flat.cov()(i,j) - (9.+16.*mass/rr)*ll(i)*ll(j) ) / 3.
345  ) / lapconf_tot / rr / rr ;
346  }
347  }
348  taij_down_rot.std_spectral_base() ;
349  taij_down_rot.annule_domain(0) ;
350  taij_down_rot.inc_dzpuis(2) ;
351 
352  Sym_tensor taij_down_bh(mp, COV, mp.get_bvect_cart()) ;
353  taij_down_bh.set_etat_qcq() ;
354 
355  for (int i=1; i<=3; i++) {
356  for (int j=1; j<=3; j++) {
357  taij_down_bh.set(i,j) = 2. * mass * pow(confo_tot,7.)
358  * pow(lapconf_auto_bh,4.) * (2.+3.*mass/rr)
359  * (flat.cov()(i,j) - (3.+4.*mass/rr) * ll(i) * ll(j))
360  / 3. / lapconf_auto / rr / rr ;
361  }
362  }
363  taij_down_bh.std_spectral_base() ;
364  taij_down_bh.annule_domain(0) ;
365  taij_down_bh.inc_dzpuis(2) ;
366 
367  Scalar taij_quad_rstot(mp) ;
368  taij_quad_rstot = 0. ;
369 
370  for (int i=1; i<=3; i++) {
371  for (int j=1; j<=3; j++) {
372  taij_quad_rstot += taij_down_rs(i,j) * taij_tot(i,j) ;
373  }
374  }
375  taij_quad_rstot.std_spectral_base() ;
376 
377  Scalar taij_quad_rsrotbh(mp) ;
378  taij_quad_rsrotbh = 0. ;
379 
380  for (int i=1; i<=3; i++) {
381  for (int j=1; j<=3; j++) {
382  taij_quad_rsrotbh += taij_tot_rs(i,j)
383  * (taij_down_rot(i,j) + taij_down_bh(i,j)) ;
384  }
385  }
386  taij_quad_rsrotbh.std_spectral_base() ;
387 
388  taij_quad_tot_rs = taij_quad_rstot + taij_quad_rsrotbh ;
390 
391  taij_quad_tot_rot = 8.*mass*mass*pow(confo_tot,14.)
392  * pow(lapconf_auto_bh,6.) * (2.+3.*mass/rr)
393  * (ll(2)*orb_rot_x - ll(1)*orb_rot_y)
394  / 3. / lapconf_tot / lapconf_tot / pow(rr,4.)
395  + 2.*mass*mass*pow(confo_tot,14.)*pow(lapconf_auto_bh,4.)
396  * (3.*(1.+2.*mass/rr)*(orb_rot_x*orb_rot_x+orb_rot_y*orb_rot_y)
397  -2.*(1.+3.*mass/rr)*(ll(2)*orb_rot_x-ll(1)*orb_rot_y)
398  *(ll(2)*orb_rot_x-ll(1)*orb_rot_y))
399  / 3. / lapconf_tot / lapconf_tot / pow(rr,4.) ;
400 
403 
404  taij_quad_tot_bh = 8.*mass*mass*pow(confo_tot,14.)
405  * pow(lapconf_auto_bh,8.) * (2.+3.*mass/rr) * (2.+3.*mass/rr)
406  / 3. / lapconf_tot / lapconf_tot / pow(rr,4.) ;
407 
410 
412  + taij_quad_tot_bh ;
414 
415  }
416  else { // Isotropic coordinates with the maximal slicing
417 
418  // Sets C/M^2 for each case of the lapse boundary condition
419  // --------------------------------------------------------
420  double cc ;
421 
422  if (bc_lapconf_nd) { // Neumann boundary condition
423  if (bc_lapconf_fs) { // First condition
424  // d(\alpha \psi)/dr = 0
425  // ---------------------
426  cc = 2. * (sqrt(13.) - 1.) / 3. ;
427  }
428  else { // Second condition
429  // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
430  // -----------------------------------------
431  cc = 4. / 3. ;
432  }
433  }
434  else { // Dirichlet boundary condition
435  if (bc_lapconf_fs) { // First condition
436  // (\alpha \psi) = 1/2
437  // -------------------
438  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
439  abort() ;
440  }
441  else { // Second condition
442  // (\alpha \psi) = 1/sqrt(2.) \psi_KS
443  // ----------------------------------
444  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
445  abort() ;
446  // cc = 2. * sqrt(2.) ;
447  }
448  }
449 
450  Scalar r_are(mp) ;
452  r_are.std_spectral_base() ;
453 
454  // Computation of \tilde{A}^{ij}
455  // -----------------------------
456 
457  Sym_tensor flat_taij(mp, CON, mp.get_bvect_cart()) ;
458  flat_taij.set_etat_qcq() ;
459 
460  for (int i=1; i<=3; i++) {
461  for (int j=1; j<=3; j++) {
462  flat_taij.set(i,j) = shift_auto_rs(j).deriv(i)
463  + shift_auto_rs(i).deriv(j) + d_shift_comp(i,j)
464  + d_shift_comp(j,i)
465  - 2. * divshift % flat.con()(i,j) / 3. ;
466  }
467  }
468 
469  flat_taij.std_spectral_base() ;
470 
471  taij_tot_rs = 0.5 * pow(confo_tot, 7.) * flat_taij / lapconf_tot ;
472 
475 
476  if (taij_tot_rs(1,2).get_etat() == ETATQCQ) {
477  for (int i=1; i<=3; i++) {
478  for (int j=1; j<=3; j++) {
479  taij_tot_rs.set(i,j).raccord(1) ;
480  }
481  }
482  }
483 
485 
486  for (int i=1; i<=3; i++) {
487  for (int j=1; j<=3; j++) {
488  taij_tot_bh.set(i,j) = pow(confo_tot,7.)*mass*mass*cc
489  * sqrt(1. - 2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
490  * (flat.con()(i,j) - 3.*ll(i)*ll(j)) / lapconf_tot
491  / pow(r_are*rr,3.) ;
492  }
493  }
494 
497 
498  for (int i=1; i<=3; i++) {
499  for (int j=1; j<=3; j++) {
500  taij_tot_bh.set(i,j).raccord(1) ;
501  }
502  }
503 
505 
507 
510 
511  for (int i=1; i<=3; i++) {
512  for (int j=1; j<=3; j++) {
513  taij_tot.set(i,j).raccord(1) ;
514  }
515  }
516 
517  for (int i=1; i<=3; i++) {
518  for (int j=1; j<=3; j++) {
519  taij_tot_rot.set(i,j) = 0. ;
520  }
521  }
523 
524  // Computation of \tilde{A}_{BH}^{ij}
525  // ----------------------------------
526 
527  Scalar divshift_auto(mp) ;
528  divshift_auto = shift_auto_rs(1).deriv(1)
529  + shift_auto_rs(2).deriv(2) + shift_auto_rs(3).deriv(3) ;
530  divshift_auto.std_spectral_base() ;
531 
532  Sym_tensor flat_taij_auto_rs(mp, CON, mp.get_bvect_cart()) ;
533  flat_taij_auto_rs.set_etat_qcq() ;
534 
535  for (int i=1; i<=3; i++) {
536  for (int j=1; j<=3; j++) {
537  flat_taij_auto_rs.set(i,j) = shift_auto_rs(j).deriv(i)
538  + shift_auto_rs(i).deriv(j)
539  - 2. * divshift_auto % flat.con()(i,j) / 3. ;
540  }
541  }
542 
543  flat_taij_auto_rs.std_spectral_base() ;
544 
545  taij_auto_rs = 0.5 * pow(confo_tot, 7.) * flat_taij_auto_rs
546  / lapconf_tot ;
547 
550 
551  if (taij_auto_rs(1,2).get_etat() == ETATQCQ) {
552  for (int i=1; i<=3; i++) {
553  for (int j=1; j<=3; j++) {
554  taij_auto_rs.set(i,j).raccord(1) ;
555  }
556  }
557  }
558 
560 
563 
564  for (int i=1; i<=3; i++) {
565  for (int j=1; j<=3; j++) {
566  taij_auto.set(i,j).raccord(1) ;
567  }
568  }
569 
570  // Computation of \tilde{A}_{NS}^{ij}
571  // ----------------------------------
572 
573  Scalar divshift_comp(mp) ;
574  divshift_comp = d_shift_comp(1,1) + d_shift_comp(2,2)
575  + d_shift_comp(3,3) ;
576  divshift_comp.std_spectral_base() ;
577 
578  Sym_tensor flat_taij_comp(mp, CON, mp.get_bvect_cart()) ;
579  flat_taij_comp.set_etat_qcq() ;
580 
581  for (int i=1; i<=3; i++) {
582  for (int j=1; j<=3; j++) {
583  flat_taij_comp.set(i,j) = d_shift_comp(i,j)
584  + d_shift_comp(j,i)
585  - 2. * divshift_comp % flat.con()(i,j) / 3. ;
586  }
587  }
588 
589  flat_taij_comp.std_spectral_base() ;
590 
591  taij_comp = 0.5 * pow(confo_comp+0.5, 7.) * flat_taij_comp
592  / (lapconf_comp+0.5) ;
593 
596 
597  if (taij_comp(1,2).get_etat() == ETATQCQ) {
598  for (int i=1; i<=3; i++) {
599  for (int j=1; j<=3; j++) {
600  taij_comp.set(i,j).raccord(1) ;
601  }
602  }
603  }
604 
605  // Computation of \tilde{A}^{ij} \tilde{A}_{ij}
606  // --------------------------------------------
607 
608  Sym_tensor flat_dshift(mp, COV, mp.get_bvect_cart()) ;
609  flat_dshift.set_etat_qcq() ;
610 
611  for (int i=1; i<=3; i++) {
612  for (int j=1; j<=3; j++) {
613  flat_dshift.set(i,j) =
614  flat.cov()(j,1)*(shift_auto_rs(1).deriv(i)
615  + d_shift_comp(i,1))
616  + flat.cov()(j,2)*(shift_auto_rs(2).deriv(i)
617  + d_shift_comp(i,2))
618  + flat.cov()(j,3)*(shift_auto_rs(3).deriv(i)
619  + d_shift_comp(i,3))
620  + flat.cov()(i,1)*(shift_auto_rs(1).deriv(j)
621  + d_shift_comp(j,1))
622  + flat.cov()(i,2)*(shift_auto_rs(2).deriv(j)
623  + d_shift_comp(j,2))
624  + flat.cov()(i,3)*(shift_auto_rs(3).deriv(j)
625  + d_shift_comp(j,3))
626  - 2. * divshift * flat.cov()(i,j) / 3. ;
627  }
628  }
629 
630  Sym_tensor taij_down_rs(mp, COV, mp.get_bvect_cart()) ;
631  taij_down_rs.set_etat_qcq() ;
632 
633  taij_down_rs = 0.5 * pow(confo_tot, 7.) * flat_dshift / lapconf_tot ;
634 
635  taij_down_rs.std_spectral_base() ;
636  taij_down_rs.annule_domain(0) ;
637 
638  Sym_tensor taij_down_bh(mp, COV, mp.get_bvect_cart()) ;
639  taij_down_bh.set_etat_qcq() ;
640 
641  for (int i=1; i<=3; i++) {
642  for (int j=1; j<=3; j++) {
643  taij_down_bh.set(i,j) = pow(confo_tot,7.)*mass*mass*cc
644  * sqrt(1. - 2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
645  * (flat.cov()(i,j) - 3.*ll(i)%ll(j)) / lapconf_tot
646  / pow(r_are*rr,3.) ;
647  }
648  }
649  taij_down_bh.std_spectral_base() ;
650  taij_down_bh.annule_domain(0) ;
651 
652  for (int i=1; i<=3; i++) {
653  for (int j=1; j<=3; j++) {
654  taij_down_bh.set(i,j).raccord(1) ;
655  }
656  }
657 
658  taij_down_bh.inc_dzpuis(2) ;
659 
660  Scalar taij_quad_rstot(mp) ;
661  taij_quad_rstot = 0. ;
662 
663  for (int i=1; i<=3; i++) {
664  for (int j=1; j<=3; j++) {
665  taij_quad_rstot += taij_down_rs(i,j) % taij_tot(i,j) ;
666  }
667  }
668  taij_quad_rstot.std_spectral_base() ;
669 
670  Scalar taij_quad_rsbh(mp) ;
671  taij_quad_rsbh = 0. ;
672 
673  for (int i=1; i<=3; i++) {
674  for (int j=1; j<=3; j++) {
675  taij_quad_rsbh += taij_tot_rs(i,j) % taij_down_bh(i,j) ;
676  }
677  }
678  taij_quad_rsbh.std_spectral_base() ;
679 
680  taij_quad_tot_rs = taij_quad_rstot + taij_quad_rsbh ;
682 
683  taij_quad_tot_rot = 0. ;
685 
686  taij_quad_tot_bh = 6.*pow(confo_tot,14.)*pow(mass*mass*cc,2.)
687  * (1. - 2.*mass/r_are/rr + cc*cc*pow(mass/r_are/rr,4.))
688  / lapconf_tot / lapconf_tot / pow(r_are*rr, 6.) ;
692 
694 
699 
700  // -------------------------
701  Scalar taij_quad_auto1(mp) ;
702  taij_quad_auto1 = 0. ;
703  for (int i=1; i<=3; i++) {
704  for (int j=1; j<=3; j++) {
705  taij_quad_auto1 += taij_auto_rs(i,j)
706  * (taij_down_rs(i,j)
707  + pow(confo_tot/(confo_comp+0.5),7.)*(lapconf_comp+0.5)
708  * taij_comp(i,j) / lapconf_tot) ;
709  }
710  }
711  taij_quad_auto1.std_spectral_base() ;
712 
713  Scalar taij_quad_auto2(mp) ;
714  taij_quad_auto2 = 0. ;
715  for (int i=1; i<=3; i++) {
716  for (int j=1; j<=3; j++) {
717  taij_quad_auto2 += taij_tot_bh(i,j) % taij_down_rs(i,j) ;
718  }
719  }
720  taij_quad_auto2.std_spectral_base() ;
721 
722  taij_quad_auto = taij_quad_auto1 + 2.*taij_quad_auto2 ;
725  if (taij_quad_auto.get_etat() == ETATQCQ) {
727  }
728 
729  // Computation of \tilde{A}_{NS}^{ij} \tilde{A}^{NS}_{ij}
730  // ------------------------------------------------------
731 
732  taij_quad_comp = 0. ;
733  for (int i=1; i<=3; i++) {
734  for (int j=1; j<=3; j++) {
735  taij_quad_comp += taij_comp(i,j) % taij_comp(i,j) ;
736  }
737  }
739 
740  }
741 
742  // The derived quantities are obsolete
743  // -----------------------------------
744 
745  del_deriv() ;
746 
747 }
748 }
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:675
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:490
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:223
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:788
Lorene prototypes.
Definition: app_hor.h:67
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:393
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:786
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:790
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
Definition: scalar.h:560
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:739
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric_flat.C:156
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:741
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:351
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:825
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:137
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:724
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:809
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:663
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:742
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:935
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:226
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: hole_bhns.C:395
Scalar lapconf_tot
Total lapconf function.
Definition: hole_bhns.h:101
Coord r
r coordinate centered on the grid
Definition: map.h:736
Scalar taij_quad_tot_rs
Part of the scalar from the numerical computation.
Definition: hole_bhns.h:224
Coord cost
Definition: map.h:740