LORENE
bin_ns_bh_glob.C
1 /*
2  * Copyright (c) 2005 Philippe Grandclement
3  *
4  * This file is part of LORENE.
5  *
6  * LORENE is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * LORENE is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with LORENE; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 
23 char bin_ns_bh_glob_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bin_ns_bh/bin_ns_bh_glob.C,v 1.8 2014/10/13 08:52:43 j_novak Exp $" ;
24 
25 /*
26  * $Id: bin_ns_bh_glob.C,v 1.8 2014/10/13 08:52:43 j_novak Exp $
27  * $Log: bin_ns_bh_glob.C,v $
28  * Revision 1.8 2014/10/13 08:52:43 j_novak
29  * Lorene classes and functions now belong to the namespace Lorene.
30  *
31  * Revision 1.7 2014/10/06 15:13:01 j_novak
32  * Modified #include directives to use c++ syntax.
33  *
34  * Revision 1.6 2007/04/24 20:13:53 f_limousin
35  * Implementation of Dirichlet and Neumann BC for the lapse
36  *
37  * Revision 1.5 2006/06/23 07:09:24 p_grandclement
38  * Addition of spinning black hole
39  *
40  * Revision 1.4 2006/06/01 12:47:52 p_grandclement
41  * update of the Bin_ns_bh project
42  *
43  * Revision 1.3 2005/12/01 12:59:10 p_grandclement
44  * Files for bin_ns_bh project
45  *
46  * Revision 1.2 2005/11/30 11:09:06 p_grandclement
47  * Changes for the Bin_ns_bh project
48  *
49  * Revision 1.1 2005/08/29 15:10:15 p_grandclement
50  * Addition of things needed :
51  * 1) For BBH with different masses
52  * 2) Provisory files for the mixted binaries (Bh and NS) : THIS IS NOT
53  * WORKING YET !!!
54  *
55  *
56  * $Header: /cvsroot/Lorene/C++/Source/Bin_ns_bh/bin_ns_bh_glob.C,v 1.8 2014/10/13 08:52:43 j_novak Exp $
57  *
58  */
59 
60 
61 
62 //standard
63 #include <cstdlib>
64 #include <cmath>
65 
66 // Lorene
67 #include "nbr_spx.h"
68 #include "tenseur.h"
69 #include "bhole.h"
70 #include "bin_ns_bh.h"
71 #include "proto.h"
72 #include "utilitaires.h"
73 #include "graphique.h"
74 #include "unites.h"
75 #include "metrique.h"
76 
77 namespace Lorene {
78 double Bin_ns_bh::adm_systeme() const {
79  Cmp der_un (hole.psi_auto().dsdr()) ;
80 
81  Map_af auxi_mp (star.get_mp()) ;
82  Cmp der_deux (star.confpsi_auto().dsdr()) ;
83 
84  double masse = hole.mp.integrale_surface_infini(der_un) +
85  auxi_mp.integrale_surface_infini(der_deux) ;
86 
87  masse /= -2*M_PI ;
88  return masse ;
89 }
90 
91 double Bin_ns_bh::adm_systeme_volume() const {
92 
93  using namespace Unites ;
94 
95  Tenseur auxi_bh (flat_scalar_prod(hole.tkij_tot, hole.tkij_auto)) ;
96  Tenseur kk_bh (hole.mp) ;
97  kk_bh = 0 ;
98  Tenseur work_bh(hole.mp) ;
99  work_bh.set_etat_qcq() ;
100  for (int i=0 ; i<3 ; i++) {
101  work_bh.set() = auxi_bh(i, i) ;
102  kk_bh = kk_bh + work_bh ;
103  }
104  Cmp integ_bh (pow(hole.psi_tot(), 5.)*kk_bh()) ;
105  integ_bh.annule(0) ;
106  integ_bh.std_base_scal() ;
107 
108  Cmp integ_hor1 (hole.psi_tot()) ;
109 
110  Cmp tet (hole.mp) ;
111  tet = hole.mp.tet ;
112  Cmp phi (hole.mp) ;
113  phi = hole.mp.phi ;
114  Tenseur rad (hole.mp, 1, COV, hole.mp.get_bvect_cart()) ;
115  rad.set_etat_qcq() ;
116  rad.set(0) = cos(phi)*sin(tet) ;
117  rad.set(1) = sin(phi)*sin(tet) ;
118  rad.set(2) = cos(tet) ;
119 
120  Cmp integ_hor2 (hole.mp) ;
121  integ_hor2.annule_hard() ;
122  integ_hor2.set_dzpuis(2) ;
123  for (int m=0 ; m<3 ; m++)
124  for (int n=0 ; n<3 ; n++)
125  integ_hor2 += rad(m)*rad(n)*hole.tkij_tot(m,n) ;
126  integ_hor2 *= pow(hole.psi_tot(),3.)/4. ;
127  integ_hor2.std_base_scal() ;
128 
129  Tenseur auxi_ns (flat_scalar_prod(star.tkij_tot, star.tkij_auto)) ;
130  Tenseur kk_ns (star.mp) ;
131  kk_ns = 0 ;
132  Tenseur work_ns(star.mp) ;
133  work_ns.set_etat_qcq() ;
134  for (int i=0 ; i<3 ; i++) {
135  work_ns.set() = auxi_ns(i, i) ;
136  kk_ns = kk_ns + work_ns ;
137  }
138  Cmp integ_ns (pow(star.confpsi_comp() + star.confpsi_auto(), 5.)*kk_ns()) ;
139  integ_ns.std_base_scal() ;
140 
141  Cmp integ_matter (pow(star.confpsi_comp()+star.confpsi_auto(), 5.)*star.ener_euler()) ;
142  integ_matter.std_base_scal() ;
143 
144  double masse = (integ_bh.integrale()+integ_ns.integrale())/16/M_PI +
145  hole.mp.integrale_surface(integ_hor1, hole.rayon)/hole.rayon/4/M_PI +
146  hole.mp.integrale_surface(integ_hor2, hole.rayon)/2/M_PI
147  + integ_matter.integrale()*ggrav ;
148 
149  return masse ;
150 }
151 
152 double Bin_ns_bh::komar_systeme() const {
153  Cmp der_un (hole.n_auto().dsdr()) ;
154 
155  Map_af auxi_mp (star.get_mp()) ;
156  Cmp der_deux (star.n_auto().dsdr()) ;
157 
158  double masse = hole.mp.integrale_surface_infini(der_un) +
159  auxi_mp.integrale_surface_infini(der_deux) ;
160 
161  masse /= 4*M_PI ;
162 
163  return masse ;
164 }
165 
166 double Bin_ns_bh::viriel() const {
167  double adm = adm_systeme() ;
168  double komar = komar_systeme() ;
169 
170  return (adm-komar)/adm ;
171 }
172 
173 double Bin_ns_bh::moment_systeme_inf() const {
174 
175  if (omega == 0)
176  return 0 ;
177  else {
178 
179  // On construit une grille et un mapping auxiliaire :
180  int nzones = hole.mp.get_mg()->get_nzone() ;
181  double* bornes = new double [nzones+1] ;
182  double courant = fabs(hole.mp.get_ori_x()-star.mp.get_ori_x())+1 ;
183  for (int i=nzones-1 ; i>0 ; i--) {
184  bornes[i] = courant ;
185  courant /= 2. ;
186  }
187  bornes[0] = 0 ;
188  bornes[nzones] = __infinity ;
189 
190  Map_af mapping (*hole.mp.get_mg(), bornes) ;
191 
192  delete [] bornes ;
193 
194  // On construit k_total dessus :
195  Tenseur_sym k_total (mapping, 2, CON, mapping.get_bvect_cart()) ;
196  k_total.set_etat_qcq() ;
197 
198  Tenseur shift_un (mapping, 1, CON, mapping.get_bvect_cart()) ;
199  shift_un.set_etat_qcq() ;
200 
201  Tenseur shift_deux (mapping, 1, CON, mapping.get_bvect_cart()) ;
202  shift_deux.set_etat_qcq() ;
203 
204  shift_un.set_triad (*hole.shift_auto.get_triad()) ;
205  shift_un.set(0).import(hole.shift_auto(0)) ;
206  shift_un.set(1).import(hole.shift_auto(1)) ;
207  shift_un.set(2).import(hole.shift_auto(2)) ;
208  shift_un.change_triad (mapping.get_bvect_cart()) ;
209 
210  shift_deux.set_triad (*star.shift_auto.get_triad()) ;
211  shift_deux.set(0).import(star.shift_auto(0)) ;
212  shift_deux.set(1).import(star.shift_auto(1)) ;
213  shift_deux.set(2).import(star.shift_auto(2)) ;
214  shift_deux.change_triad(mapping.get_bvect_cart()) ;
215 
216  Tenseur shift_tot (shift_un+shift_deux) ;
217  shift_tot.set_std_base() ;
218  shift_tot.annule(0, nzones-2) ;
219  // On enleve les residus
220  shift_tot.inc2_dzpuis() ;
221  shift_tot.dec2_dzpuis() ;
222 
223  Tenseur grad (shift_tot.gradient()) ;
224  Tenseur trace (grad(0, 0)+grad(1, 1)+grad(2, 2)) ;
225  for (int i=0 ; i<3 ; i++) {
226  k_total.set(i, i) = grad(i, i)-trace()/3. ;
227  for (int j=i+1 ; j<3 ; j++)
228  k_total.set(i, j) = (grad(i, j)+grad(j, i))/2. ;
229  }
230 
231  for (int lig=0 ; lig<3 ; lig++)
232  for (int col=lig ; col<3 ; col++)
233  k_total.set(lig, col).mult_r_zec() ;
234 
235  Tenseur vecteur_un (mapping, 1, CON, mapping.get_bvect_cart()) ;
236  vecteur_un.set_etat_qcq() ;
237  for (int i=0 ; i<3 ; i++)
238  vecteur_un.set(i) = k_total(0, i) ;
239  vecteur_un.change_triad (mapping.get_bvect_spher()) ;
240  Cmp integrant_un (vecteur_un(0)) ;
241 
242  Tenseur vecteur_deux (mapping, 1, CON, mapping.get_bvect_cart()) ;
243  vecteur_deux.set_etat_qcq() ;
244  for (int i=0 ; i<3 ; i++)
245  vecteur_deux.set(i) = k_total(1, i) ;
246  vecteur_deux.change_triad (mapping.get_bvect_spher()) ;
247  Cmp integrant_deux (vecteur_deux(0)) ;
248 
249  // On multiplie par y et x :
250  integrant_un.va = integrant_un.va.mult_st() ;
251  integrant_un.va = integrant_un.va.mult_sp() ;
252 
253  integrant_deux.va = integrant_deux.va.mult_st() ;
254  integrant_deux.va = integrant_deux.va.mult_cp() ;
255 
256  double moment = mapping.integrale_surface_infini (-integrant_un+integrant_deux) ;
257 
258  moment /= 8*M_PI ;
259  return moment ;
260  }
261 }
262 
263 double Bin_ns_bh::moment_systeme_hor() const {
264 
265  using namespace Unites ;
266 
267  if (omega == 0)
268  return 0 ;
269  else {
270  //Contribution du trou noir :
271  Cmp xa_bh (hole.mp) ;
272  xa_bh = hole.mp.xa ;
273  xa_bh.std_base_scal() ;
274 
275  Cmp ya_bh (hole.mp) ;
276  ya_bh = hole.mp.ya ;
277  ya_bh.std_base_scal() ;
278 
279  Tenseur vecteur_bh (hole.mp, 1, CON, hole.mp.get_bvect_cart()) ;
280  vecteur_bh.set_etat_qcq() ;
281  for (int i=0 ; i<3 ; i++)
282  vecteur_bh.set(i) = (-ya_bh*hole.tkij_tot(0, i)+
283  xa_bh * hole.tkij_tot(1, i)) ;
284  vecteur_bh.set_std_base() ;
285  vecteur_bh.annule(hole.mp.get_mg()->get_nzone()-1) ;
286  vecteur_bh.change_triad (hole.mp.get_bvect_spher()) ;
287 
288  Cmp integrant_bh (pow(hole.psi_tot(), 6)*vecteur_bh(0)) ;
289  integrant_bh.std_base_scal() ;
290  double moment_bh = hole.mp.integrale_surface
291  (integrant_bh, hole.rayon)/8/M_PI ;
292 
293  // Contribution de l'étoile :
294  Cmp xa_ns (star.mp) ;
295  xa_ns = star.mp.xa ;
296  xa_ns.std_base_scal() ;
297 
298  Cmp ya_ns (star.mp) ;
299  ya_ns = star.mp.ya ;
300  ya_ns.std_base_scal() ;
301 
302  Cmp integrant_ns (pow(star.confpsi_auto()+star.confpsi_comp(), 10)*(star.ener_euler()+star.press())*
303  (xa_ns*star.u_euler(1) - ya_ns*star.u_euler(0))) ;
304  integrant_ns.std_base_scal() ;
305 
306  double moment_ns = integrant_ns.integrale() * ggrav ;
307  return moment_ns + moment_bh ;
308  }
309 }
310 
311 Tbl Bin_ns_bh::linear_momentum_systeme_inf() const {
312 
313  Tbl res (3) ;
314  res.set_etat_qcq() ;
315 
316  // On construit une grille et un mapping auxiliaire :
317  int nzones = hole.mp.get_mg()->get_nzone() ;
318  double* bornes = new double [nzones+1] ;
319  double courant = fabs(hole.mp.get_ori_x()-star.mp.get_ori_x())+1 ;
320  for (int i=nzones-1 ; i>0 ; i--) {
321  bornes[i] = courant ;
322  courant /= 2. ;
323  }
324  bornes[0] = 0 ;
325  bornes[nzones] = __infinity ;
326 
327  Map_af mapping (*hole.mp.get_mg(), bornes) ;
328 
329  delete [] bornes ;
330 
331  // On construit k_total dessus :
332  Tenseur_sym k_total (mapping, 2, CON, mapping.get_bvect_cart()) ;
333  k_total.set_etat_qcq() ;
334 
335  Tenseur shift_un (mapping, 1, CON, mapping.get_bvect_cart()) ;
336  shift_un.set_etat_qcq() ;
337 
338  Tenseur shift_deux (mapping, 1, CON, mapping.get_bvect_cart()) ;
339  shift_deux.set_etat_qcq() ;
340 
341  shift_un.set_triad (*hole.shift_auto.get_triad()) ;
342  shift_un.set(0).import(hole.shift_auto(0)) ;
343  shift_un.set(1).import(hole.shift_auto(1)) ;
344  shift_un.set(2).import(hole.shift_auto(2)) ;
345  shift_un.change_triad (mapping.get_bvect_cart()) ;
346 
347  shift_deux.set_triad (*star.shift_auto.get_triad()) ;
348  shift_deux.set(0).import(star.shift_auto(0)) ;
349  shift_deux.set(1).import(star.shift_auto(1)) ;
350  shift_deux.set(2).import(star.shift_auto(2)) ;
351  shift_deux.change_triad(mapping.get_bvect_cart()) ;
352 
353  shift_un.set_std_base() ;
354  shift_deux.set_std_base() ;
355 
356  Tenseur shift_tot (shift_un+shift_deux) ;
357  shift_tot.set_std_base() ;
358  shift_tot.annule(0, nzones-2) ;
359 
360  Cmp compy (shift_tot(1)) ;
361  compy.mult_r_zec() ;
362 
363  int nr = mapping.get_mg()->get_nr(nzones-1) ;
364  int nt = mapping.get_mg()->get_nt(nzones-1) ;
365  int np = mapping.get_mg()->get_np(nzones-1) ;
366  Tbl val_inf (nt*np) ;
367  val_inf.set_etat_qcq() ;
368  for (int k=0 ; k<np ; k++)
369  for (int j=0 ; j<nt ; j++)
370  val_inf.set(k*nt + j) = fabs(compy (nzones-1, k, j, nr-1)) ;
371 
372  Tenseur grad (shift_tot.gradient()) ;
373  Tenseur trace (grad(0, 0)+grad(1, 1)+grad(2, 2)) ;
374  for (int i=0 ; i<3 ; i++) {
375  k_total.set(i, i) = grad(i, i)-trace()/3. ;
376  for (int j=i+1 ; j<3 ; j++)
377  k_total.set(i, j) = (grad(i, j)+grad(j, i))/2. ;
378  }
379 
380  for (int comp=0 ; comp<3 ; comp++) {
381  Tenseur vecteur (mapping, 1, CON, mapping.get_bvect_cart()) ;
382  vecteur.set_etat_qcq() ;
383  for (int i=0 ; i<3 ; i++)
384  vecteur.set(i) = k_total(i, comp) ;
385  vecteur.change_triad (mapping.get_bvect_spher()) ;
386  Cmp integrant (vecteur(0)) ;
387 
388  res.set(comp) = mapping.integrale_surface_infini (integrant)/8/M_PI ;
389  }
390  return res ;
391 }
392 
393 double Bin_ns_bh::distance_propre_axe_bh (const int nr) const {
394 
395  double x_bh = hole.mp.get_ori_x() + hole.rayon ;
396 
397  // Les coefficients du changement de variable :
398  double pente = -2./x_bh ;
399  double constante = - 1. ;
400 
401  double ksi ; // variable d'integration.
402  double xabs ; // x reel.
403  double air_un, theta_un, phi_un ; // coordonnee spheriques 1
404  double air_deux, theta_deux, phi_deux ; // coordonnee spheriques 2
405 
406  double* coloc = new double[nr] ;
407  double* coef = new double[nr] ;
408  int* deg = new int[3] ;
409  deg[0] = 1 ; deg[1] = 1 ; deg[2] = nr ;
410 
411  for (int i=0 ; i<nr ; i++) {
412  ksi = -cos (M_PI*i/(nr-1)) ;
413  xabs = (ksi+constante)/pente ;
414 
415  hole.mp.convert_absolute (xabs, 0, 0, air_un, theta_un, phi_un) ;
416  star.mp.convert_absolute (xabs, 0, 0, air_deux, theta_deux, phi_deux) ;
417 
418  coloc[i] = pow (hole.psi_auto().val_point (air_un, theta_un, phi_un) +
419  star.confpsi_auto().val_point (air_deux, theta_deux, phi_deux), 2.) ;
420  }
421 
422  // On prend les coefficients de la fonction
423  cfrcheb(deg, deg, coloc, deg, coef) ;
424 
425  // On integre
426  double* som = new double[nr] ;
427  som[0] = 2 ;
428  for (int i=2 ; i<nr ; i+=2)
429  som[i] = 1./(i+1)-1./(i-1) ;
430  for (int i=1 ; i<nr ; i+=2)
431  som[i] = 0 ;
432 
433  double res = 0 ;
434  for (int i=0 ; i<nr ; i++)
435  res += som[i]*coef[i] ;
436 
437  res /= pente ;
438 
439  delete [] deg ;
440  delete [] coef ;
441  delete [] coloc ;
442  delete [] som ;
443 
444  return res ;
445 }
446 
447 double Bin_ns_bh::distance_propre_axe_ns (const int nr) const {
448 
449  double x_ns = star.mp.get_ori_x() - star.mp.val_r (star.nzet, -1, M_PI/2, M_PI) ;
450 
451  // Les coefficients du changement de variable :
452  double pente = 2./x_ns ;
453  double constante = - 1. ;
454 
455  double ksi ; // variable d'integration.
456  double xabs ; // x reel.
457  double air_un, theta_un, phi_un ; // coordonnee spheriques 1
458  double air_deux, theta_deux, phi_deux ; // coordonnee spheriques 2
459 
460  double* coloc = new double[nr] ;
461  double* coef = new double[nr] ;
462  int* deg = new int[3] ;
463  deg[0] = 1 ; deg[1] = 1 ; deg[2] = nr ;
464 
465  for (int i=0 ; i<nr ; i++) {
466  ksi = -cos (M_PI*i/(nr-1)) ;
467  xabs = (ksi-constante)/pente ;
468 
469  hole.mp.convert_absolute (xabs, 0, 0, air_un, theta_un, phi_un) ;
470  star.mp.convert_absolute (xabs, 0, 0, air_deux, theta_deux, phi_deux) ;
471 
472  coloc[i] = pow (hole.psi_auto().val_point (air_un, theta_un, phi_un) +
473  star.confpsi_auto().val_point (air_deux, theta_deux, phi_deux), 2.) ;
474  }
475 
476  // On prend les coefficients de la fonction
477  cfrcheb(deg, deg, coloc, deg, coef) ;
478 
479  // On integre
480  double* som = new double[nr] ;
481  som[0] = 2 ;
482  for (int i=2 ; i<nr ; i+=2)
483  som[i] = 1./(i+1)-1./(i-1) ;
484  for (int i=1 ; i<nr ; i+=2)
485  som[i] = 0 ;
486 
487  double res = 0 ;
488  for (int i=0 ; i<nr ; i++)
489  res += som[i]*coef[i] ;
490 
491  res /= pente ;
492 
493  delete [] deg ;
494  delete [] coef ;
495  delete [] coloc ;
496  delete [] som ;
497 
498  return res ;
499 }
500 
501 double Bin_ns_bh::smarr() const {
502 
503  using namespace Unites ;
504 
505  // The tests
506  Tenseur psiq_t (pow(star.get_confpsi()(), 4.)) ;
507  psiq_t.set_std_base() ;
508 
509  Tenseur_sym furmet (star.mp, 2, CON, star.mp.get_bvect_cart()) ;
510  furmet.set_etat_qcq() ;
511  for (int i=0 ; i< 3 ; i++) {
512  furmet.set(i,i) = 1/psiq_t() ;
513  for(int j=i+1 ; j<3 ; j++)
514  furmet.set(i,j) = 0 ;
515  }
516  Metrique met (furmet, false) ;
517 
518  Tenseur_sym kij (star.get_tkij_tot()/psiq_t) ;
519  kij.change_triad(star.mp.get_bvect_cart()) ;
520  kij.dec2_dzpuis() ;
521  Tenseur_sym kij_cov (manipule (kij, met)) ;
522  Tenseur shift (star.get_shift()) ;
523  shift.change_triad(star.mp.get_bvect_cart()) ;
524 
525  Tenseur aime (star.mp, 1, CON, star.mp.get_bvect_cart()) ;
526  aime.set_etat_qcq() ;
527  aime.set(0) = -omega*star.mp.ya ;
528  aime.set(1) = omega*star.mp.xa ;
529  aime.set(2) = 0 ;
530  aime.annule(star.mp.get_mg()->get_nzone()-1) ;
531  aime.set_std_base() ;
532  shift = shift - aime ;
533 
534  // La matière :
535  Tenseur u_euler (star.get_u_euler()) ;
536  u_euler.change_triad(star.mp.get_bvect_cart()) ;
537  Tenseur u_i_bas (manipule(u_euler, met)) ;
538  Tenseur mat (qpig*(star.get_nnn()*(star.get_ener_euler() + star.get_s_euler()) - 2*(star.get_ener_euler()+star.get_press())*contract(u_i_bas, 0, shift, 0))) ;
539 
540  // La partie avec la matière :
541  Cmp psiq (pow(star.get_confpsi()(), 4.)) ;
542 
543  Cmp integ_matter (star.get_nnn()()*(star.get_ener_euler()() + star.get_s_euler()())
544  - 2*(star.get_ener_euler()()+star.get_press()())*psiq*flat_scalar_prod(u_euler, shift)()) ;
545  integ_matter = integ_matter * pow(star.get_confpsi()(),6.) ;
546  integ_matter.std_base_scal() ;
547  integ_matter.annule(star.get_nzet(), star.get_mp().get_mg()->get_nzone()-1) ;
548  double matter_term = integ_matter.integrale()*qpig/4/M_PI ;
549 
550  // Integrale sur horizon :
551  Cmp tet (hole.mp) ;
552  tet = hole.mp.tet ;
553  Cmp phi (hole.mp) ;
554  phi = hole.mp.phi ;
555  Tenseur rad (hole.mp, 1, COV, hole.mp.get_bvect_cart()) ;
556  rad.set_etat_qcq() ;
557  rad.set(0) = cos(phi)*sin(tet) ;
558  rad.set(1) = sin(phi)*sin(tet) ;
559  rad.set(2) = cos(tet) ;
560 
561  Cmp temp (hole.mp) ;
562  temp.annule_hard() ;
563  temp.set_dzpuis(2) ;
564  for (int m=0 ; m<3 ; m++)
565  for (int n=0 ; n<3 ; n++)
566  temp += rad(m)*rad(n)*hole.tkij_tot(m,n) ;
567  temp *= pow(hole.psi_tot(),2.) ;
568  temp.std_base_scal() ;
569 
570  Cmp integ_hor ((hole.get_n_tot()().dsdr()-hole.get_n_tot()()*temp)
571  *pow(hole.get_psi_tot()(), 2)) ;
572  integ_hor.std_base_scal() ;
573  integ_hor.raccord(1) ;
574  double hor_term = hole.mp.integrale_surface(integ_hor, hole.get_rayon()) ;
575  hor_term /= 4*M_PI ;
576 
577  double m_test = hor_term + matter_term + 2*omega*moment_systeme_inf() +
578  2*(hole.omega_local-omega)*hole.local_momentum() ;
579 
580  return m_test ;
581  }
582 }
Coord xa
Absolute x coordinate.
Definition: map.h:730
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
const Map & get_mp() const
Returns the mapping.
Definition: etoile.h:659
double integrale_surface_infini(const Cmp &ci) const
Performs the surface integration of ci at infinity.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:783
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
const Tenseur & get_shift() const
Returns the total shift vector .
Definition: etoile.h:730
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:768
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
Tenseur confpsi_comp
Part of the conformal factor $$ generated principaly by the companion star.
Definition: et_bin_nsbh.h:109
const Tenseur & get_psi_tot() const
Returns the total .
Definition: bhole.h:424
int get_nzet() const
Returns the number of domains occupied by the star.
Definition: etoile.h:662
double omega_local
local angular velocity
Definition: bhole.h:276
Tenseur psi_auto
Part of generated by the hole.
Definition: bhole.h:290
const Tenseur & get_n_tot() const
Returns the total N .
Definition: bhole.h:407
Tenseur press
Fluid pressure.
Definition: etoile.h:461
void set(const Map *mp, Mtbl *(*construct)(const Map *))
Semi-constructor from a mapping and a method.
Definition: coord.C:134
Tenseur n_auto
Part of the lapse { N} generated principaly by the star.
Definition: et_bin_nsbh.h:85
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:94
Coord tet
coordinate centered on the grid
Definition: map.h:719
Tenseur shift_auto
Part of generated by the hole.
Definition: bhole.h:297
Tenseur shift_auto
Part of the shift vector generated principaly by the star.
Definition: etoile.h:889
Coord phi
coordinate centered on the grid
Definition: map.h:720
const Tenseur & get_press() const
Returns the fluid pressure.
Definition: etoile.h:682
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: etoile.h:474
Et_bin_nsbh star
The neutron star.
Definition: bin_ns_bh.h:128
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tenseur.h:701
Map_af & mp
Affine mapping.
Definition: bhole.h:273
const Tenseur & get_s_euler() const
Returns the trace of the stress tensor in the Eulerian frame.
Definition: etoile.h:688
Map & mp
Mapping associated with the star.
Definition: etoile.h:429
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
double get_rayon() const
Returns the radius of the horizon.
Definition: bhole.h:347
const Tenseur & get_nnn() const
Returns the total lapse function N.
Definition: etoile.h:727
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition: bin_ns_bh.h:136
const Tenseur & get_confpsi() const
Returns the part of the conformal factor $$.
Definition: et_bin_nsbh.h:257
double rayon
Radius of the horizon in LORENE&#39;s units.
Definition: bhole.h:274
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Tenseur_sym tkij_auto
Auto .
Definition: bhole.h:307
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
const Tenseur & get_ener_euler() const
Returns the total energy density with respect to the Eulerian observer.
Definition: etoile.h:685
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:644
int nzet
Number of domains of *mp occupied by the star.
Definition: etoile.h:432
Tenseur_sym tkij_tot
Total .
Definition: bhole.h:308
Coord ya
Absolute y coordinate.
Definition: map.h:731
Tenseur confpsi_auto
Part of the conformal factor $$ generated principaly by the star.
Definition: et_bin_nsbh.h:104
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
const Tenseur_sym & get_tkij_tot() const
Returns the total extrinsic curvature tensor $K^{ij}$ generated by { shift_auto} and { shift_comp}...
Definition: et_bin_nsbh.h:308
Tenseur psi_tot
Total .
Definition: bhole.h:292
Tenseur_sym tkij_tot
Total extrinsic curvature tensor $K^{ij}$ generated by { shift_auto} and { shift_comp}.
Definition: et_bin_nsbh.h:152
Bhole hole
The black hole.
Definition: bin_ns_bh.h:131
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:69
Tenseur n_auto
Part of N generated by the hole.
Definition: bhole.h:286
Tenseur_sym tkij_auto
Part of the extrinsic curvature tensor $K^{ij}$ generated by { shift_auto}.
Definition: et_bin_nsbh.h:146
const Tenseur & get_u_euler() const
Returns the fluid 3-velocity with respect to the Eulerian observer.
Definition: etoile.h:694
void convert_absolute(double xx, double yy, double zz, double &rr, double &theta, double &pphi) const
Determines the coordinates corresponding to given absolute Cartesian coordinates (X...
Definition: map.C:302
Tenseur ener_euler
Total energy density in the Eulerian frame.
Definition: etoile.h:465
Tenseur manipule(const Tenseur &, const Metrique &, int idx)
Raise or lower the index idx depending on its type, using the given Metrique .