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
23char 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
77namespace Lorene {
78double 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
91double 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
152double 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
166double Bin_ns_bh::viriel() const {
167 double adm = adm_systeme() ;
168 double komar = komar_systeme() ;
169
170 return (adm-komar)/adm ;
171}
172
173double 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
263double 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
311Tbl 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
393double 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
447double 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
501double 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}
Lorene::Bin_ns_bh::star
Et_bin_nsbh star
The neutron star.
Definition bin_ns_bh.h:128
Lorene::Bin_ns_bh::hole
Bhole hole
The black hole.
Definition bin_ns_bh.h:131
Lorene::Bin_ns_bh::omega
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition bin_ns_bh.h:136
Lorene::Cmp::std_base_scal
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition cmp.C:644
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:69
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
Lorene::flat_scalar_prod
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...
Definition tenseur_operateur.C:653
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition tenseur_operateur.C:279
Lorene::manipule
Tenseur manipule(const Tenseur &, const Metrique &, int idx)
Raise or lower the index idx depending on its type, using the given Metrique .
Definition tenseur_operateur.C:509
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::phi
Coord phi
coordinate centered on the grid
Definition map.h:720
Lorene::tet
Coord tet
coordinate centered on the grid
Definition map.h:719