LORENE
bin_hor.C
1 /*
2  * Methods of class Bin_hor
3  *
4  * (see file bin_hor.h for documentation)
5  *
6  */
7 
8 /*
9  * Copyright (c) 2004-2005 Francois Limousin
10  * Jose Luis Jaramillo
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License as published by
16  * the Free Software Foundation; either version 2 of the License, or
17  * (at your option) any later version.
18  *
19  * LORENE is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with LORENE; if not, write to the Free Software
26  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27  *
28  */
29 
30 
31 char bin_hor_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bin_hor/bin_hor.C,v 1.12 2014/10/13 08:52:42 j_novak Exp $" ;
32 
33 /*
34  * $Id: bin_hor.C,v 1.12 2014/10/13 08:52:42 j_novak Exp $
35  * $Log: bin_hor.C,v $
36  * Revision 1.12 2014/10/13 08:52:42 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.11 2014/10/06 15:13:00 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.10 2007/04/13 15:28:55 f_limousin
43  * Lots of improvements, generalisation to an arbitrary state of
44  * rotation, implementation of the spatial metric given by Samaya.
45  *
46  * Revision 1.9 2006/08/01 14:37:19 f_limousin
47  * New version
48  *
49  * Revision 1.8 2006/06/29 08:51:00 f_limousin
50  * *** empty log message ***
51  *
52  * Revision 1.7 2006/06/28 13:36:09 f_limousin
53  * Convergence to a given irreductible mass
54  *
55  * Revision 1.6 2006/05/24 16:56:37 f_limousin
56  * Many small modifs.
57  *
58  * Revision 1.5 2005/06/13 15:47:29 jl_jaramillo
59  * Add some quatities in write_global()
60  *
61  * Revision 1.4 2005/06/09 16:12:04 f_limousin
62  * Implementation of amg_mom_adm().
63  *
64  * Revision 1.3 2005/04/29 14:02:44 f_limousin
65  * Important changes : manage the dependances between quantities (for
66  * instance psi and psi4). New function write_global(ost).
67  *
68  * Revision 1.2 2005/03/04 09:38:41 f_limousin
69  * Implement the constructor from a file, operator>>, operator<<
70  * and function sauve.
71  *
72  * Revision 1.1 2004/12/29 16:11:02 f_limousin
73  * First version
74  *
75  *
76  * $Header: /cvsroot/Lorene/C++/Source/Bin_hor/bin_hor.C,v 1.12 2014/10/13 08:52:42 j_novak Exp $
77  *
78  */
79 
80 //standard
81 #include <cstdlib>
82 #include <cmath>
83 
84 // Lorene
85 #include "nbr_spx.h"
86 #include "tenseur.h"
87 #include "tensor.h"
88 #include "isol_hor.h"
89 #include "proto.h"
90 #include "utilitaires.h"
91 //#include "graphique.h"
92 
93 // Standard constructor
94 // --------------------
95 
96 namespace Lorene {
98  hole1(mp1), hole2(mp2), omega(0){
99 
100  holes[0] = &hole1 ;
101  holes[1] = &hole2 ;
102 }
103 
104 // Copy constructor
105 // ----------------
106 
107 Bin_hor::Bin_hor (const Bin_hor& source) :
108  hole1(source.hole1), hole2(source.hole2), omega(source.omega) {
109 
110  holes[0] = &hole1 ;
111  holes[1] = &hole2 ;
112  }
113 
114 // Constructor from a file
115 // -----------------------
116 
117 Bin_hor::Bin_hor(Map_af& mp1, Map_af& mp2, FILE* fich)
118  : hole1(mp1, fich),
119  hole2(mp2, fich),
120  omega(0) {
121 
122  fread_be(&omega, sizeof(double), 1, fich) ;
123  holes[0] = &hole1 ;
124  holes[1] = &hole2 ;
125 
126 }
127 
128  //--------------//
129  // Destructor //
130  //--------------//
131 
133 }
134 
135  //-----------------------//
136  // Mutators / assignment //
137  //-----------------------//
138 
139 void Bin_hor::operator= (const Bin_hor& source) {
140  hole1 = source.hole1 ;
141  hole2 = source.hole2 ;
142 
143  omega = source.omega ;
144 }
145 
146 
147  //--------------------------//
148  // Save in a file //
149  //--------------------------//
150 
151 void Bin_hor::sauve(FILE* fich) const {
152 
153  hole1.sauve(fich) ;
154  hole2.sauve(fich) ;
155  fwrite_be (&omega, sizeof(double), 1, fich) ;
156 
157 }
158 
159 
160 //Initialisation : Sum of two static BH
162  set_omega (0) ;
163  hole1.init_bhole() ;
164  hole2.init_bhole() ;
165 
168 
171 
172  decouple() ;
174 
175 }
176 
177 
178 void Bin_hor::write_global(ostream& ost, double lim_nn, int bound_nn,
179  int bound_psi, int bound_beta, double alpha) const {
180 
181  double distance = hole1.get_mp().get_ori_x() - hole2.get_mp().get_ori_x() ;
182  double mass_adm = adm_mass() ;
183  double mass_komar = komar_mass() ;
184  double mass_area = sqrt(hole1.area_hor()/16/M_PI) +
185  sqrt(hole2.area_hor()/16/M_PI) ;
186  double J_adm = ang_mom_adm() ;
187  double J_hor = ang_mom_hor() ; //hole1.ang_mom_hor() + hole2.ang_mom_hor() ;
188  double j1 = hole1.ang_mom_hor() ;
189  double j2 = hole2.ang_mom_hor() ;
190  double mass_ih1 = hole1.mass_hor() ;
191  double mass_ih2 = hole2.mass_hor() ;
192  double mass_ih = mass_ih1 + mass_ih2 ;
193  double omega1 = hole1.omega_hor() ;
194  double omega2 = hole2.omega_hor() ;
195 
196  // Verification of Smarr :
197  // -----------------------
198 
199  // Les alignemenents pour le signe des CL.
200  double orientation1 = hole1.mp.get_rot_phi() ;
201  assert ((orientation1 == 0) || (orientation1 == M_PI)) ;
202  int aligne1 = (orientation1 == 0) ? 1 : -1 ;
203 
204  Vector angular1 (hole1.mp, CON, hole1.mp.get_bvect_cart()) ;
205  Scalar yya1 (hole1.mp) ;
206  yya1 = hole1.mp.ya ;
207  Scalar xxa1 (hole1.mp) ;
208  xxa1 = hole1.mp.xa ;
209 
210  angular1.set(1) = aligne1 * omega * yya1 ;
211  angular1.set(2) = - aligne1 * omega * xxa1 ;
212  angular1.set(3).annule_hard() ;
213 
214  angular1.set(1).set_spectral_va()
215  .set_base(*(hole1.mp.get_mg()->std_base_vect_cart()[0])) ;
216  angular1.set(2).set_spectral_va()
217  .set_base(*(hole1.mp.get_mg()->std_base_vect_cart()[1])) ;
218  angular1.set(3).set_spectral_va()
219  .set_base(*(hole1.mp.get_mg()->std_base_vect_cart()[2])) ;
220 
221  angular1.change_triad(hole1.mp.get_bvect_spher()) ;
222 
223 
224  double orientation2 = hole2.mp.get_rot_phi() ;
225  assert ((orientation2 == 0) || (orientation2 == M_PI)) ;
226  int aligne2 = (orientation2 == 0) ? 1 : -1 ;
227 
228  Vector angular2 (hole2.mp, CON, hole2.mp.get_bvect_cart()) ;
229  Scalar yya2 (hole2.mp) ;
230  yya2 = hole2.mp.ya ;
231  Scalar xxa2 (hole2.mp) ;
232  xxa2 = hole2.mp.xa ;
233 
234  angular2.set(1) = aligne2 * omega * yya2 ;
235  angular2.set(2) = - aligne2 * omega * xxa2 ;
236  angular2.set(3).annule_hard() ;
237 
238  angular2.set(1).set_spectral_va()
239  .set_base(*(hole2.mp.get_mg()->std_base_vect_cart()[0])) ;
240  angular2.set(2).set_spectral_va()
241  .set_base(*(hole2.mp.get_mg()->std_base_vect_cart()[1])) ;
242  angular2.set(3).set_spectral_va()
243  .set_base(*(hole2.mp.get_mg()->std_base_vect_cart()[2])) ;
244 
245  angular2.change_triad(hole2.mp.get_bvect_spher()) ;
246 
247 
248  Scalar btilde1 (hole1.b_tilde() -
249  contract(angular1, 0, hole1.tgam.radial_vect().up_down(hole1.tgam), 0)) ;
250  Scalar btilde2 (hole2.b_tilde() -
251  contract(angular2, 0, hole2.tgam.radial_vect().up_down(hole2.tgam), 0)) ;
252 
253 
254 
255 
256  Vector integrand_un (hole1.mp, COV, hole1.mp.get_bvect_spher()) ;
257  integrand_un = hole1.nn.derive_cov(hole1.ff)*pow(hole1.psi, 2)
258  - btilde1*contract(hole1.get_k_dd(), 1,
259  hole1.tgam.radial_vect(), 0)*pow(hole1.psi, 2) ;
260  integrand_un.std_spectral_base() ;
261 
262  Vector integrand_deux (hole2.mp, COV, hole2.mp.get_bvect_spher()) ;
263  integrand_deux = hole2.nn.derive_cov(hole2.ff)*pow(hole2.psi, 2)
264  - btilde2*contract(hole2.get_k_dd(), 1,
265  hole2.tgam.radial_vect(), 0)*pow(hole2.psi, 2) ;
266  integrand_deux.std_spectral_base() ;
267 
268  double horizon = hole1.mp.integrale_surface(integrand_un(1),
269  hole1.get_radius())+
270  hole2.mp.integrale_surface(integrand_deux(1), hole2.get_radius()) ;
271 
272  horizon /= 4*M_PI ;
273 
274  double J_smarr = (mass_komar - horizon) / 2. / omega ;
275 
276  ost.precision(8) ;
277  ost << "# Grid : " << hole1.mp.get_mg()->get_nr(1) << "x"
278  << hole1.mp.get_mg()->get_nt(1) << "x"
279  << hole1.mp.get_mg()->get_np(1) << " R_out(l) : " ;
280 
281  for (int ll=0; ll<hole1.mp.get_mg()->get_nzone(); ll++) {
282  ost << " " << hole1.mp.val_r(ll, 1., M_PI/2, 0) ;
283  }
284  ost << endl ;
285  ost << "# bound N, lim N : " << bound_nn << " " << lim_nn
286  << " - bound Psi : " << bound_psi << " - bound shift : " << bound_beta
287  << " alpha = " << alpha << endl ;
288 
289  ost << "# distance omega Mass_ADM Mass_K M_area J_ADM J_hor" << endl ;
290  ost << distance << " " ;
291  ost << omega << " " ;
292  ost << mass_adm << " " ;
293  ost << mass_komar << " " ;
294  ost << mass_area << " " ;
295  ost << J_adm << " " ;
296  ost << J_hor << endl ;
297  ost << "# mass_ih1 mass_ih2 mass_ih j1 J2 omega1 omega2" << endl ;
298  ost << mass_ih1 << " " ;
299  ost << mass_ih2 << " " ;
300  ost << mass_ih << " " ;
301  ost << j1 << " " ;
302  ost << j2 << " " ;
303  ost << omega1 << " " ;
304  ost << omega2 << endl ;
305  ost << "# ADM_mass/M_area J/M_area2 omega*M_area" << endl ;
306  ost << mass_adm / mass_area << " " ;
307  ost << J_adm /mass_area / mass_area << " " ;
308  ost << omega * mass_area << endl ;
309  ost << "# Diff J_hor/J_ADM Diff J_ADM/J_Smarr Diff J_hor/J_smarr"
310  << endl ;
311  ost << fabs(J_adm - J_hor) / J_adm << " " << fabs(J_adm - J_smarr) / J_adm
312  << " " << fabs(J_hor - J_smarr) / J_hor << endl ;
313 
314 
315 }
316 
317 }
Lorene::Bin_hor::~Bin_hor
virtual ~Bin_hor()
Destructor.
Definition: bin_hor.C:132
Lorene::Single_hor::nn
Scalar nn
Lapse function .
Definition: isol_hor.h:921
Lorene::Bin_hor::set_omega
void set_omega(double ome)
Sets the orbital velocity to ome.
Definition: isol_hor.h:1409
Lorene::Bin_hor::ang_mom_adm
double ang_mom_adm() const
Calculates the angular momentum of the black hole.
Definition: binhor_global.C:184
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Lorene::Map::get_ori_x
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:768
Lorene::Single_hor::ff
Metric_flat ff
3 metric flat
Definition: isol_hor.h:980
Lorene::Single_hor::psi
Scalar psi
Conformal factor .
Definition: isol_hor.h:930
Lorene::Bin_hor::decouple
void decouple()
Calculates decouple which is used to obtain tkij_auto and tkij_comp.
Definition: binhor_kij.C:473
Lorene::Map::get_rot_phi
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition: map.h:775
Lorene::Bin_hor::hole1
Single_hor hole1
Black hole one.
Definition: isol_hor.h:1342
Lorene::Single_hor::mass_hor
double mass_hor() const
Mass computed at the horizon
Definition: single_param.C:135
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Map::get_bvect_spher
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::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Single_hor::sauve
virtual void sauve(FILE *fich) const
Total or partial saves in a binary file.
Definition: single_hor.C:244
Lorene::Scalar::derive_cov
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
Lorene::Map_af::integrale_surface
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
Definition: map_af_integ_surf.C:90
Lorene::Tensor::up_down
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
Definition: tensor_calculus.C:305
Lorene::Bin_hor::adm_mass
double adm_mass() const
Calculates the ADM mass of the system.
Definition: binhor_global.C:83
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::Single_hor::get_radius
double get_radius() const
Returns the radius of the horizon.
Definition: isol_hor.h:1046
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Lorene::Bin_hor::holes
Single_hor * holes[2]
Array on the black holes.
Definition: isol_hor.h:1346
Lorene::Map_af::val_r
virtual double val_r(int l, double xi, double theta, double pphi) const
Returns the value of the radial coordinate r for a given in a given domain.
Definition: map_af_radius.C:96
Lorene::Single_hor::omega_hor
double omega_hor() const
Orbital velocity
Definition: single_param.C:160
Lorene::Vector::change_triad
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition: vector_change_triad.C:75
Lorene::Bin_hor::omega
double omega
Angular velocity.
Definition: isol_hor.h:1348
Lorene::Bin_hor::ang_mom_hor
double ang_mom_hor() const
Calculates the angular momentum of the black hole using the formula at the horizon.
Definition: binhor_global.C:118
Lorene::Scalar::set_spectral_va
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:604
Lorene::Single_hor::mp
Map_af & mp
Affine mapping.
Definition: isol_hor.h:900
Lorene::Single_hor::init_bhole
void init_bhole()
Sets the values of the fields to :
Definition: single_hor.C:484
Lorene::Vector::set
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
Lorene::Bin_hor::operator=
void operator=(const Bin_hor &)
Affectation operator.
Definition: bin_hor.C:139
Lorene::Single_hor::b_tilde
const Scalar b_tilde() const
Radial component of the shift with respect to the conformal metric.
Definition: single_param.C:68
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Single_hor::get_mp
const Map_af & get_mp() const
Returns the mapping (readonly).
Definition: isol_hor.h:1038
Lorene::Bin_hor::sauve
void sauve(FILE *fich) const
Total or partial saves in a binary file.
Definition: bin_hor.C:151
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Lorene::fread_be
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition: fread_be.C:69
Lorene::Map_af
Affine radial mapping.
Definition: map.h:2027
Lorene::Valeur::set_base
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition: valeur.C:810
Lorene::Bin_hor
Definition: isol_hor.h:1337
Lorene::Map::get_bvect_cart
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition: map.h:791
Lorene::Bin_hor::init_bin_hor
void init_bin_hor()
Initialisation of the system.
Definition: bin_hor.C:161
Lorene::Mg3d::std_base_vect_cart
Base_val ** std_base_vect_cart() const
Returns the standard spectral bases for the Cartesian components of a vector.
Definition: mg3d_std_base.C:177
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Lorene::Bin_hor::hole2
Single_hor hole2
Black hole two.
Definition: isol_hor.h:1343
Lorene::Vector
Tensor field of valence 1.
Definition: vector.h:188
Lorene::Metric::radial_vect
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity.
Definition: metric.C:362
Lorene::Map::xa
Coord xa
Absolute x coordinate.
Definition: map.h:730
Lorene::Single_hor::psi_comp_import
void psi_comp_import(const Single_hor &comp)
Imports the part of due to the companion hole comp .
Definition: single_hor.C:424
Lorene::Bin_hor::write_global
void write_global(ostream &, double lim_nn, int bound_nn, int bound_psi, int bound_beta, double alpha) const
Write global quantities in a formatted file.
Definition: bin_hor.C:178
Lorene::Single_hor::get_k_dd
const Sym_tensor & get_k_dd() const
k_dd
Definition: single_hor.C:348
Lorene::Scalar::annule_hard
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:380
Lorene::fwrite_be
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition: fwrite_be.C:70
Lorene::Bin_hor::extrinsic_curvature
void extrinsic_curvature()
Calculation of the extrinsic curvature tensor.
Definition: binhor_kij.C:88
Lorene::Bin_hor::Bin_hor
Bin_hor(Map_af &mp1, Map_af &mp2)
Standard constructor.
Definition: bin_hor.C:97
Lorene::Single_hor::n_comp_import
void n_comp_import(const Single_hor &comp)
Imports the part of N due to the companion hole comp .
Definition: single_hor.C:390
Lorene::Single_hor::ang_mom_hor
double ang_mom_hor() const
Angular momentum (modulo)
Definition: single_param.C:107
Lorene::Single_hor::tgam
Metric tgam
3 metric tilde
Definition: isol_hor.h:977
Lorene::Bin_hor::komar_mass
double komar_mass() const
Calculates the Komar mass of the system using : .
Definition: binhor_global.C:101
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition: tenseur_operateur.C:279
Lorene::Single_hor::area_hor
double area_hor() const
Area of the horizon.
Definition: single_param.C:88
Lorene::Map::ya
Coord ya
Absolute y coordinate.
Definition: map.h:731
Lorene::Vector::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: vector.C:316