LORENE
binary_global_xcts.C
1/*
2 * Methods of class Binary_xcts to compute global quantities
3 * (see file binary_xcts.h for documentation)
4 */
5
6/*
7 * Copyright (c) 2010 Michal Bejger
8 *
9 * This file is part of LORENE.
10 *
11 * LORENE is free software; you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License version 2
13 * as published by the Free Software Foundation.
14 *
15 * LORENE is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with LORENE; if not, write to the Free Software
22 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 *
24 */
25
26char binary_global_xcts_C[] = "$Header: /cvsroot/Lorene/C++/Source/Binary_xcts/binary_global_xcts.C,v 1.11 2014/10/13 08:52:45 j_novak Exp $" ;
27
28/*
29 * $Id: binary_global_xcts.C,v 1.11 2014/10/13 08:52:45 j_novak Exp $
30 * $Log: binary_global_xcts.C,v $
31 * Revision 1.11 2014/10/13 08:52:45 j_novak
32 * Lorene classes and functions now belong to the namespace Lorene.
33 *
34 * Revision 1.10 2010/12/21 11:15:46 m_bejger
35 * Linear momentum properly defined
36 *
37 * Revision 1.9 2010/12/20 15:45:40 m_bejger
38 * Spectral basis in lin_mom() was not defined
39 *
40 * Revision 1.8 2010/12/20 09:54:09 m_bejger
41 * Angular momentum correction, stub for linear momentum added
42 *
43 * Revision 1.7 2010/12/09 10:39:41 m_bejger
44 * Further corrections to integral quantities
45 *
46 * Revision 1.6 2010/10/26 19:16:26 m_bejger
47 * Cleanup of some diagnostic messages
48 *
49 * Revision 1.5 2010/10/25 15:02:08 m_bejger
50 * mass_kom_vol() corrected
51 *
52 * Revision 1.4 2010/10/24 21:45:24 m_bejger
53 * mass_adm() corrected
54 *
55 * Revision 1.3 2010/06/17 14:48:14 m_bejger
56 * Minor corrections
57 *
58 * Revision 1.2 2010/06/04 19:54:19 m_bejger
59 * Minor corrections, mass volume integrals need to be checked out
60 *
61 * Revision 1.1 2010/05/04 07:35:54 m_bejger
62 * Initial version
63 *
64 * $Header: /cvsroot/Lorene/C++/Source/Binary_xcts/binary_global_xcts.C,v 1.11 2014/10/13 08:52:45 j_novak Exp $
65 *
66 */
67
68// Headers C
69#include "math.h"
70
71// Headers Lorene
72#include "nbr_spx.h"
73#include "binary_xcts.h"
74#include "unites.h"
75#include "metric.h"
76
77 //---------------------------------//
78 // ADM mass //
79 //---------------------------------//
80
81namespace Lorene {
82double Binary_xcts::mass_adm() const {
83
84 using namespace Unites ;
85
86 if (p_mass_adm == 0x0) { // a new computation is required
87
88 p_mass_adm = new double ;
89 *p_mass_adm = 0 ;
90
91 const Map_af map0 (et[0]->get_mp()) ;
92 const Metric& flat = (et[0]->get_flat()) ;
93
94 Vector dpsi((et[0]->get_Psi()).derive_cov(flat)) ;
95
96 dpsi.change_triad(map0.get_bvect_spher()) ;
97
98 Scalar integrand ( dpsi(1) ) ;
99
100 *p_mass_adm = - 2.* map0.integrale_surface_infini(integrand)/ qpig ;
101
102 }
103
104 return *p_mass_adm ;
105
106}
107
108 //---------------------------------//
109 // ADM mass //
110 // (volume integral) //
111 //---------------------------------//
112
114
115 using namespace Unites ;
116
117 double massadm = 0. ;
118
119 for (int i=0; i<=1; i++) { // loop on the stars
120
121 // Declaration of all fields
122
123 const Scalar& psi(et[i]->get_Psi()) ;
124 Scalar psi5 = pow(psi, 5.) ;
125 psi5.std_spectral_base() ;
126
127 Scalar spsi7 = pow(psi, -7.) ;
128 spsi7.std_spectral_base() ;
129
130 const Scalar& ener_euler = et[i]->get_ener_euler() ;
131 const Scalar& hacar_auto = et[i]->get_hacar_auto() ;
132 const Scalar& hacar_comp = et[i]->get_hacar_comp() ;
133
134 Scalar source = psi5 % ener_euler
135 + spsi7 % (hacar_auto + hacar_comp)/(4.*qpig) ;
136
137 source.std_spectral_base() ;
138
139 massadm += source.integrale() ;
140 }
141
142 return massadm ;
143}
144
145 //---------------------------------//
146 // Komar mass //
147 //---------------------------------//
148
149double Binary_xcts::mass_kom() const {
150
151 using namespace Unites ;
152
153 if (p_mass_kom == 0x0) { // a new computation is requireed
154
155 p_mass_kom = new double ;
156 *p_mass_kom = 0 ;
157
158 const Scalar& logn = et[0]->get_logn() ;
159 const Metric& flat = et[0]->get_flat() ;
160
161 Map_af map0 (et[0]->get_mp()) ;
162
163 Vector vect = logn.derive_con(flat) ;
164 vect.change_triad(map0.get_bvect_spher()) ;
165
166 Scalar integrant (vect(1)) ;
167
168 *p_mass_kom = map0.integrale_surface_infini (integrant) / qpig ;
169
170 } // End of the case where a new computation was necessary
171
172 return *p_mass_kom ;
173
174}
175
177
178 using namespace Unites ;
179
180 double masskom = 0.;
181
182 for (int i=0; i<=1; i++) { // loop on the stars
183
184 const Scalar& Psi = et[i]->get_Psi() ;
185 const Scalar& psi4 = et[i]->get_psi4() ;
186 const Scalar& chi = et[i]->get_chi() ;
187
188 const Scalar& ener_euler = et[i]->get_ener_euler() ;
189 const Scalar& s_euler = et[i]->get_s_euler() ;
190 const Scalar& hacar_auto = et[i]->get_hacar_auto() ;
191 const Scalar& hacar_comp = et[i]->get_hacar_comp() ;
192
193 Scalar psi4chi = psi4 % chi ;
194 psi4chi.std_spectral_base() ;
195
196 Scalar source = 0.5 * ener_euler * (psi4chi + psi4 % Psi)
197 + psi4chi * s_euler + pow(Psi, -7.) * (7.*chi/Psi + 1.)
198 * (hacar_auto + hacar_comp) / (8.*qpig) ;
199 source.std_spectral_base() ;
200
201 masskom += source.integrale() ;
202
203 }
204
205 return masskom ;
206
207}
208
209 //-------------------------------------//
210 // Total angular momentum (z-axis) //
211 //-------------------------------------//
212
214
215 using namespace Unites ;
216
217 if (p_angu_mom == 0x0) { // a new computation is required
218
219 p_angu_mom = new Tbl(3) ;
220 p_angu_mom->annule_hard() ; // fills the double array with zeros
221
222 // Reference Cartesian vector basis of the Absolute frame
223 Base_vect_cart bvect_ref(0.) ; // 0. = parallel to the Absolute frame
224
225 for (int i=0; i<=1; i++) { // loop on the stars
226
227 const Map& mp = et[i]->get_mp() ;
228
229 // Azimuthal vector d/dphi
230 Vector vphi(mp, CON, bvect_ref) ;
231 Scalar yya (mp) ; yya = mp.ya ;
232 Scalar xxa (mp) ; xxa = mp.xa ;
233 vphi.set(1) = - yya ; // phi^X
234 vphi.set(2) = xxa ;
235 vphi.set(3) = 0 ;
236
237 vphi.std_spectral_base() ;
238 vphi.change_triad(mp.get_bvect_cart()) ;
239
240 // Matter part
241 // -----------
242 const Scalar& ee = et[i]->get_ener_euler() ;
243 const Scalar& pp = et[i]->get_press() ;
244
245 Vector jmom = pow(et[i]->get_Psi(), 10) * (ee + pp)
246 * (et[i]->get_u_euler()) ;
247 jmom.std_spectral_base() ;
248
249 const Metric& flat = et[i]->get_flat() ;
250 Vector vphi_cov = vphi.up_down(flat) ;
251
252 Scalar integrand = contract(jmom, 0, vphi_cov, 0) ;
253
254 p_angu_mom->set(2) += integrand.integrale() ;
255
256 } // End of the loop on the stars
257
258 } // End of the case where a new computation was necessary
259
260 return *p_angu_mom ;
261
262}
263
264 //---------------------------------//
265 // Total linear momentum //
266 //---------------------------------//
267
268const Tbl& Binary_xcts::lin_mom() const {
269
270 using namespace Unites ;
271
272 if (p_lin_mom == 0x0) { // a new computation is required
273
274 p_lin_mom = new Tbl(3) ;
276
277 // Reference Cartesian vector basis of the Absolute frame
278 Base_vect_cart bvect_ref(0.) ; // 0. = parallel to the Absolute frame
279
280 for (int i=0; i<=1; i++) { // loop on the stars
281
282 const Scalar& ee = et[i]->get_ener_euler() ;
283 const Scalar& pp = et[i]->get_press() ;
284 Vector lmom = pow(et[i]->get_Psi(), 10) * (ee + pp)
285 * ( et[i]->get_u_euler() ) ;
286
287 lmom.std_spectral_base() ;
288 lmom.change_triad(bvect_ref) ;
289
290 // loop on the components
291 for (int j=1; j<=2; j++)
292 p_lin_mom->set(j-1) += lmom(j).integrale() ;
293
294
295 } // End of the loop on the stars
296
297 } // End of the case where a new computation was necessary
298
299 return *p_lin_mom ;
300
301}
302
303 //---------------------------------//
304 // Total energy //
305 //---------------------------------//
306
308
309 if (p_total_ener == 0x0) { // a new computation is requireed
310
311 p_total_ener = new double ;
312
314
315 } // End of the case where a new computation was necessary
316
317 return *p_total_ener ;
318
319}
320
321
322 //---------------------------------//
323 // Error on the virial theorem //
324 //---------------------------------//
325
326double Binary_xcts::virial() const {
327
328 if (p_virial == 0x0) { // a new computation is requireed
329
330 p_virial = new double ;
331
332 *p_virial = 1. - mass_kom() / mass_adm() ;
333
334 }
335
336 return *p_virial ;
337
338}
339
340 //---------------------------------//
341 // Error on the virial theorem //
342 // (volume version) //
343 //---------------------------------//
344
346
347 if (p_virial_vol == 0x0) { // a new computation is requireed
348
349 p_virial_vol = new double ;
350
351 *p_virial_vol = 1. - mass_kom_vol() / mass_adm_vol() ;
352
353 }
354
355 return *p_virial_vol ;
356
357}
358}
Cartesian vectorial bases (triads).
Definition base_vect.h:201
const Tbl & lin_mom() const
Total linear momentum.
double total_ener() const
Total energy (excluding the rest mass energy).
double virial_vol() const
Estimates the relative error on the virial theorem (volume version)
double * p_mass_kom
Total Komar mass of the system.
Definition binary_xcts.h:94
Star_bin_xcts star2
Second star of the system.
Definition binary_xcts.h:69
double mass_adm_vol() const
Total ADM mass (computed by a volume integral)
double * p_virial_vol
Virial theorem error (volume version)
Tbl * p_lin_mom
Total linear momentum of the system.
double mass_kom_vol() const
Total Komar mass (computed by a volume integral)
Star_bin_xcts star1
First star of the system.
Definition binary_xcts.h:66
double mass_adm() const
Total ADM mass.
double * p_mass_adm
Total ADM mass of the system.
Definition binary_xcts.h:91
double * p_total_ener
Total energy of the system.
double mass_kom() const
Total Komar mass.
double virial() const
Estimates the relative error on the virial theorem.
Tbl * p_angu_mom
Total angular momentum of the system.
Definition binary_xcts.h:97
Star_bin_xcts * et[2]
Array of the two stars (to perform loops on the stars): et[0] contains the address of star1 and et[1]...
Definition binary_xcts.h:75
double * p_virial
Virial theorem error.
const Tbl & angu_mom() const
Total angular momentum.
Affine radial mapping.
Definition map.h:2027
double integrale_surface_infini(const Cmp &ci) const
Performs the surface integration of ci at infinity.
Base class for coordinate mappings.
Definition map.h:670
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
Coord ya
Absolute y coordinate.
Definition map.h:731
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
Coord xa
Absolute x coordinate.
Definition map.h:730
Metric for tensor calculation.
Definition metric.h:90
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
double integrale() const
Computes the integral over all space of *this .
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
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
const Scalar & get_chi() const
Return the function .
Definition star.h:1392
const Scalar & get_hacar_comp() const
Returns the part of generated by beta_comp.
Definition star.h:1422
const Scalar & get_Psi() const
Return the conformal factor .
Definition star.h:1389
const Scalar & get_hacar_auto() const
Returns the part of generated by beta_auto.
Definition star.h:1417
const Scalar & get_psi4() const
Return the conformal factor .
Definition star.h:1395
const Metric & get_flat() const
Return the flat metric defined on the mapping (Spherical components with respect to the mapping of th...
Definition star.h:1400
virtual double mass_b() const
Baryon mass.
const Map & get_mp() const
Returns the mapping.
Definition star.h:355
const Scalar & get_logn() const
Returns the logarithm of the lapse N.
Definition star.h:396
const Scalar & get_ener_euler() const
Returns the total energy density with respect to the Eulerian observer.
Definition star.h:376
const Scalar & get_s_euler() const
Returns the trace of the stress tensor in the Eulerian frame.
Definition star.h:379
const Scalar & get_press() const
Returns the fluid pressure.
Definition star.h:373
Basic array class.
Definition tbl.h:161
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition tbl.C:372
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition tbl.h:281
Tensor field of valence 1.
Definition vector.h:188
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition vector.C:316
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Scalar & set(int)
Read/write access to a component.
Definition vector.C:296
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.