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 
26 char 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 
81 namespace Lorene {
82 double 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 
149 double 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 
213 const Tbl& Binary_xcts::angu_mom() const {
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 
268 const 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 
307 double Binary_xcts::total_ener() const {
308 
309  if (p_total_ener == 0x0) { // a new computation is requireed
310 
311  p_total_ener = new double ;
312 
313  *p_total_ener = mass_adm() - star1.mass_b() - star2.mass_b() ;
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 
326 double 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 
345 double Binary_xcts::virial_vol() const {
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)
Definition: binary_xcts.h:109
Tbl * p_lin_mom
Total linear momentum of the system.
Definition: binary_xcts.h:100
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.
Definition: binary_xcts.h:103
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.
Definition: binary_xcts.h:106
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
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
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 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 .
Definition: scalar_integ.C:61
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...
Definition: scalar_deriv.C:402
const Scalar & get_Psi() const
Return the conformal factor .
Definition: star.h:1389
const Scalar & get_psi4() const
Return the conformal factor .
Definition: star.h:1395
const Scalar & get_chi() const
Return the function .
Definition: star.h:1392
const Scalar & get_hacar_auto() const
Returns the part of generated by beta_auto.
Definition: star.h:1417
const Scalar & get_hacar_comp() const
Returns the part of generated by beta_comp.
Definition: star.h:1422
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 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
const Scalar & get_ener_euler() const
Returns the total energy density with respect to the Eulerian observer.
Definition: star.h:376
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
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.