LORENE
binaire.C
1 /*
2  * Methods of class Binaire
3  *
4  * (see file binaire.h for documentation)
5  */
6 
7 /*
8  * Copyright (c) 2000-2003 Eric Gourgoulhon
9  *
10  * This file is part of LORENE.
11  *
12  * LORENE is free software; you can redistribute it and/or modify
13  * it under the terms of the GNU General Public License as published by
14  * the Free Software Foundation; either version 2 of the License, or
15  * (at your option) any later version.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 
29 
30 
31 /*
32  * $Id: binaire.C,v 1.9 2016/12/05 16:17:47 j_novak Exp $
33  * $Log: binaire.C,v $
34  * Revision 1.9 2016/12/05 16:17:47 j_novak
35  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
36  *
37  * Revision 1.8 2014/10/13 08:52:44 j_novak
38  * Lorene classes and functions now belong to the namespace Lorene.
39  *
40  * Revision 1.7 2004/03/25 10:28:59 j_novak
41  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
42  *
43  * Revision 1.6 2003/09/15 20:24:24 e_gourgoulhon
44  * Improvements in the function write_global.
45  *
46  * Revision 1.5 2003/09/15 15:10:12 e_gourgoulhon
47  * Added the member function write_global.
48  *
49  * Revision 1.4 2002/12/17 21:18:46 e_gourgoulhon
50  * Suppression of Etoile_bin::set_companion.
51  *
52  * Revision 1.3 2001/12/20 13:03:24 k_taniguchi
53  * Addition of the Komar mass, the virial error by Gourgoulhon and Bonazzola, and the virial error by Friedman, Uryu, and Shibata.
54  *
55  * Revision 1.2 2001/12/04 21:27:52 e_gourgoulhon
56  *
57  * All writing/reading to a binary file are now performed according to
58  * the big endian convention, whatever the system is big endian or
59  * small endian, thanks to the functions fwrite_be and fread_be
60  *
61  * Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
62  * LORENE
63  *
64  * Revision 2.7 2001/06/25 12:55:46 eric
65  * Traitement des etoiles compagnons (appel de Etoile_bin::set_companion dans
66  * les constructeurs).
67  *
68  * Revision 2.6 2000/07/07 14:10:17 eric
69  * AJout de display_poly.
70  *
71  * Revision 2.5 2000/03/13 14:25:52 eric
72  * Ajout des membres p_ham_constr et p_mom_constr.
73  *
74  * Revision 2.4 2000/02/18 14:52:44 eric
75  * Ajout des membres p_virial et p_total_ener.
76  * p_masse_adm --> p_mass_adm.
77  *
78  * Revision 2.3 2000/02/12 17:36:25 eric
79  * Ajout de la fonction separation().
80  *
81  * Revision 2.2 2000/02/04 17:14:47 eric
82  * Ajout du membre ref_triad.
83  *
84  * Revision 2.1 2000/02/02 10:40:36 eric
85  * Modif affichage.
86  *
87  * Revision 2.0 2000/01/31 15:57:49 eric
88  * *** empty log message ***
89  *
90  *
91  * $Header: /cvsroot/Lorene/C++/Source/Binaire/binaire.C,v 1.9 2016/12/05 16:17:47 j_novak Exp $
92  *
93  */
94 
95 // Headers C
96 #include "math.h"
97 
98 // Headers Lorene
99 #include "binaire.h"
100 #include "eos.h"
101 #include "utilitaires.h"
102 #include "unites.h"
103 
104  //--------------//
105  // Constructors //
106  //--------------//
107 
108 // Standard constructor
109 // --------------------
110 
111 namespace Lorene {
112 Binaire::Binaire(Map& mp1, int nzet1, const Eos& eos1, int irrot1,
113  Map& mp2, int nzet2, const Eos& eos2, int irrot2,
114  int relat)
115  : ref_triad(0., "Absolute frame Cartesian basis"),
116  star1(mp1, nzet1, relat, eos1, irrot1, ref_triad),
117  star2(mp2, nzet2, relat, eos2, irrot2, ref_triad)
118 {
119 
120  et[0] = &star1 ;
121  et[1] = &star2 ;
122 
123  omega = 0 ;
124  x_axe = 0 ;
125 
126  // Pointers of derived quantities initialized to zero :
127  set_der_0x0() ;
128 }
129 
130 // Copy constructor
131 // ----------------
132 Binaire::Binaire(const Binaire& bibi)
133  : ref_triad(0., "Absolute frame Cartesian basis"),
134  star1(bibi.star1),
135  star2(bibi.star2),
136  omega(bibi.omega),
137  x_axe(bibi.x_axe)
138 {
139  et[0] = &star1 ;
140  et[1] = &star2 ;
141 
142  // Pointers of derived quantities initialized to zero :
143  set_der_0x0() ;
144 }
145 
146 // Constructor from a file
147 // -----------------------
148 Binaire::Binaire(Map& mp1, const Eos& eos1, Map& mp2, const Eos& eos2,
149  FILE* fich)
150  : ref_triad(0., "Absolute frame Cartesian basis"),
151  star1(mp1, eos1, ref_triad, fich),
152  star2(mp2, eos2, ref_triad, fich)
153 {
154  et[0] = &star1 ;
155  et[1] = &star2 ;
156 
157  // omega and x_axe are read in the file:
158  fread_be(&omega, sizeof(double), 1, fich) ;
159  fread_be(&x_axe, sizeof(double), 1, fich) ;
160 
161  // Pointers of derived quantities initialized to zero :
162  set_der_0x0() ;
163 
164 }
165 
166  //------------//
167  // Destructor //
168  //------------//
169 
170 Binaire::~Binaire(){
171 
172  del_deriv() ;
173 
174 }
175 
176  //----------------------------------//
177  // Management of derived quantities //
178  //----------------------------------//
179 
180 void Binaire::del_deriv() const {
181 
182  if (p_mass_adm != 0x0) delete p_mass_adm ;
183  if (p_mass_kom != 0x0) delete p_mass_kom ;
184  if (p_angu_mom != 0x0) delete p_angu_mom ;
185  if (p_total_ener != 0x0) delete p_total_ener ;
186  if (p_virial != 0x0) delete p_virial ;
187  if (p_virial_gb != 0x0) delete p_virial_gb ;
188  if (p_virial_fus != 0x0) delete p_virial_fus ;
189  if (p_ham_constr != 0x0) delete p_ham_constr ;
190  if (p_mom_constr != 0x0) delete p_mom_constr ;
191 
192  set_der_0x0() ;
193 }
194 
195 
196 
197 
198 void Binaire::set_der_0x0() const {
199 
200  p_mass_adm = 0x0 ;
201  p_mass_kom = 0x0 ;
202  p_angu_mom = 0x0 ;
203  p_total_ener = 0x0 ;
204  p_virial = 0x0 ;
205  p_virial_gb = 0x0 ;
206  p_virial_fus = 0x0 ;
207  p_ham_constr = 0x0 ;
208  p_mom_constr = 0x0 ;
209 
210 }
211 
212 
213  //--------------//
214  // Assignment //
215  //--------------//
216 
217 // Assignment to another Binaire
218 // -----------------------------
219 
220 void Binaire::operator=(const Binaire& bibi) {
221 
222  assert( bibi.ref_triad == ref_triad ) ;
223 
224  star1 = bibi.star1 ;
225  star2 = bibi.star2 ;
226 
227  omega = bibi.omega ;
228  x_axe = bibi.x_axe ;
229 
230  // ref_triad remains unchanged
231 
232  del_deriv() ; // Deletes all derived quantities
233 
234 }
235 
236  //--------------//
237  // Outputs //
238  //--------------//
239 
240 // Save in a file
241 // --------------
242 void Binaire::sauve(FILE* fich) const {
243 
244  star1.sauve(fich) ;
245  star2.sauve(fich) ;
246 
247  fwrite_be(&omega, sizeof(double), 1, fich) ;
248  fwrite_be(&x_axe, sizeof(double), 1, fich) ;
249 
250 }
251 
252 // Printing
253 // --------
254 ostream& operator<<(ostream& ost, const Binaire& bibi) {
255  bibi >> ost ;
256  return ost ;
257 }
258 
259 
260 ostream& Binaire::operator>>(ostream& ost) const {
261 
262  using namespace Unites ;
263 
264  ost << endl ;
265  ost << "Binary system" << endl ;
266  ost << "=============" << endl ;
267  ost << endl <<
268  "Orbital angular velocity : " << omega * f_unit << " rad/s" << endl ;
269  ost << endl <<
270  "Coordinate separation between the two stellar centers : "
271  << separation() / km << " km" << endl ;
272  ost <<
273  "Absolute coordinate X of the rotation axis : " << x_axe / km
274  << " km" << endl ;
275  ost << endl << "Star 1 : " << endl ;
276  ost << "====== " << endl ;
277  ost << star1 << endl ;
278  ost << "Star 2 : " << endl ;
279  ost << "====== " << endl ;
280  ost << star2 << endl ;
281  return ost ;
282 }
283 
284 // Display in polytropic units
285 // ---------------------------
286 
287 void Binaire::display_poly(ostream& ost) const {
288 
289  using namespace Unites ;
290 
291  const Eos* p_eos1 = &( star1.get_eos() ) ;
292  const Eos_poly* p_eos_poly = dynamic_cast<const Eos_poly*>( p_eos1 ) ;
293 
294  if ((p_eos_poly != 0x0) && ( star1.get_eos() == star2.get_eos() )) {
295 
296  double kappa = p_eos_poly->get_kap() ;
297  double gamma = p_eos_poly->get_gam() ; ;
298  double kap_ns2 = pow( kappa, 0.5 /(gamma-1) ) ;
299 
300  // Polytropic unit of length in terms of r_unit :
301  double r_poly = kap_ns2 / sqrt(ggrav) ;
302 
303  // Polytropic unit of time in terms of t_unit :
304  double t_poly = r_poly ;
305 
306  // Polytropic unit of mass in terms of m_unit :
307  double m_poly = r_poly / ggrav ;
308 
309  // Polytropic unit of angular momentum in terms of j_unit :
310  double j_poly = r_poly * r_poly / ggrav ;
311 
312  ost.precision(10) ;
313  ost << endl << "Quantities in polytropic units : " << endl ;
314  ost << "==============================" << endl ;
315  ost << " ( r_poly = " << r_poly / km << " km )" << endl ;
316  ost << " d_e_max : " << separation() / r_poly << endl ;
317  ost << " d_G : "
318  << ( star2.xa_barycenter() - star1.xa_barycenter() ) / r_poly
319  << endl ;
320  ost << " Omega : " << omega * t_poly << endl ;
321  ost << " J : " << angu_mom()(2) / j_poly << endl ;
322  ost << " M_ADM : " << mass_adm() / m_poly << endl ;
323  ost << " M_Komar : " << mass_kom() / m_poly << endl ;
324  ost << " E : " << total_ener() / m_poly << endl ;
325  ost << " M_bar(star 1) : " << star1.mass_b() / m_poly << endl ;
326  ost << " M_bar(star 2) : " << star2.mass_b() / m_poly << endl ;
327  ost << " R_0(star 1) : " <<
328  0.5 * ( star1.ray_eq() + star1.ray_eq_pi() ) / r_poly << endl ;
329  ost << " R_0(star 2) : " <<
330  0.5 * ( star2.ray_eq() + star2.ray_eq_pi() ) / r_poly << endl ;
331 
332  }
333 
334 
335 }
336 
337 void Binaire::write_global(ostream& ost) const {
338 
339  using namespace Unites ;
340 
341  const Map& mp1 = star1.get_mp() ;
342  const Mg3d* mg1 = mp1.get_mg() ;
343  int nz1 = mg1->get_nzone() ;
344 
345  ost.precision(5) ;
346  ost << "# Grid 1 : " << nz1 << "x"
347  << mg1->get_nr(0) << "x" << mg1->get_nt(0) << "x" << mg1->get_np(0)
348  << " R_out(l) [km] : " ;
349  for (int l=0; l<nz1; l++) {
350  ost << " " << mp1.val_r(l, 1., M_PI/2, 0) / km ;
351  }
352  ost << endl ;
353 
354  ost << "# VE(M) "
355  << " VE(GB) "
356  << " VE(FUS) " << endl ;
357 
358  ost.setf(ios::scientific) ;
359  ost.width(14) ;
360  ost << virial() ; ost.width(14) ;
361  ost << virial_gb() ; ost.width(14) ;
362  ost << virial_fus() << endl ;
363 
364  ost << "# d [km] "
365  << " d_G [km] "
366  << " d/(a1 +a1') "
367  << " f [Hz] "
368  << " M_ADM [M_sol] "
369  << " J [G M_sol^2/c] " << endl ;
370 
371  ost.precision(14) ;
372  ost.width(20) ;
373  ost << separation() / km ; ost.width(22) ;
374  ost << ( star2.xa_barycenter() - star1.xa_barycenter() ) / km ; ost.width(22) ;
375  ost << separation() / (star1.ray_eq() + star2.ray_eq()) ; ost.width(22) ;
376  ost << omega / (2*M_PI)* f_unit ; ost.width(22) ;
377  ost << mass_adm() / msol ; ost.width(22) ;
378  ost << angu_mom()(2)/ ( qpig / (4* M_PI) * msol*msol) << endl ;
379 
380  ost << "# H_c(1)[c^2] "
381  << " e_c(1)[rho_nuc] "
382  << " M_B(1) [M_sol] "
383  << " r_eq(1) [km] "
384  << " a2/a1(1) "
385  << " a3/a1(1) " << endl ;
386 
387  ost.width(20) ;
388  ost << star1.get_ent()()(0,0,0,0) ; ost.width(22) ;
389  ost << star1.get_ener()()(0,0,0,0) ; ost.width(22) ;
390  ost << star1.mass_b() / msol ; ost.width(22) ;
391  ost << star1.ray_eq() / km ; ost.width(22) ;
392  ost << star1.ray_eq_pis2() / star1.ray_eq() ; ost.width(22) ;
393  ost << star1.ray_pole() / star1.ray_eq() << endl ;
394 
395  ost << "# H_c(2)[c^2] "
396  << " e_c(2)[rho_nuc] "
397  << " M_B(2) [M_sol] "
398  << " r_eq(2) [km] "
399  << " a2/a1(2) "
400  << " a3/a1(2) " << endl ;
401 
402  ost.width(20) ;
403  ost << star2.get_ent()()(0,0,0,0) ; ost.width(22) ;
404  ost << star2.get_ener()()(0,0,0,0) ; ost.width(22) ;
405  ost << star2.mass_b() / msol ; ost.width(22) ;
406  ost << star2.ray_eq() / km ; ost.width(22) ;
407  ost << star2.ray_eq_pis2() / star1.ray_eq() ; ost.width(22) ;
408  ost << star2.ray_pole() / star1.ray_eq() << endl ;
409 
410  // Quantities in polytropic units if the EOS is a polytropic one
411  // -------------------------------------------------------------
412  const Eos* p_eos1 = &( star1.get_eos() ) ;
413  const Eos_poly* p_eos_poly = dynamic_cast<const Eos_poly*>( p_eos1 ) ;
414 
415  if ((p_eos_poly != 0x0) && ( star1.get_eos() == star2.get_eos() )) {
416 
417  double kappa = p_eos_poly->get_kap() ;
418  double gamma = p_eos_poly->get_gam() ; ;
419  double kap_ns2 = pow( kappa, 0.5 /(gamma-1.) ) ;
420 
421  // Polytropic unit of length in terms of r_unit :
422  double r_poly = kap_ns2 / sqrt(ggrav) ;
423 
424  // Polytropic unit of time in terms of t_unit :
425  double t_poly = r_poly ;
426 
427  // Polytropic unit of mass in terms of m_unit :
428  double m_poly = r_poly / ggrav ;
429 
430  // Polytropic unit of angular momentum in terms of j_unit :
431  double j_poly = r_poly * r_poly / ggrav ;
432 
433  ost << "# d [poly] "
434  << " d_G [poly] "
435  << " Omega [poly] "
436  << " M_ADM [poly] "
437  << " J [poly] "
438  << " M_B(1) [poly] "
439  << " M_B(2) [poly] " << endl ;
440 
441  ost.width(20) ;
442  ost << separation() / r_poly ; ost.width(22) ;
443  ost << ( star2.xa_barycenter() - star1.xa_barycenter() ) / r_poly ; ost.width(22) ;
444  ost << omega * t_poly ; ost.width(22) ;
445  ost << mass_adm() / m_poly ; ost.width(22) ;
446  ost << angu_mom()(2) / j_poly ; ost.width(22) ;
447  ost << star1.mass_b() / m_poly ; ost.width(22) ;
448  ost << star2.mass_b() / m_poly << endl ;
449 
450  }
451 
452 }
453 
454 
455 
456  //-------------------------------//
457  // Miscellaneous //
458  //-------------------------------//
459 
460 double Binaire::separation() const {
461 
462  double dx = star1.get_mp().get_ori_x() - star2.get_mp().get_ori_x() ;
463  double dy = star1.get_mp().get_ori_y() - star2.get_mp().get_ori_y() ;
464  double dz = star1.get_mp().get_ori_z() - star2.get_mp().get_ori_z() ;
465 
466  return sqrt( dx*dx + dy*dy + dz*dz ) ;
467 
468 }
469 }
double total_ener() const
Total energy (excluding the rest mass energy).
virtual void sauve(FILE *) const
Save in a file.
Definition: etoile_bin.C:562
const Map & get_mp() const
Returns the mapping.
Definition: etoile.h:662
Etoile_bin star1
First star of the system.
Definition: binaire.h:118
double mass_adm() const
Total ADM mass.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
double ray_eq_pis2() const
Coordinate radius at , [r_unit].
double get_ori_y() const
Returns the y coordinate of the origin.
Definition: map.h:782
Lorene prototypes.
Definition: app_hor.h:67
Tbl * p_mom_constr
Relative error on the momentum constraint.
Definition: binaire.h:166
void display_poly(ostream &) const
Display in polytropic units.
Definition: binaire.C:287
Standard units of space, time and mass.
Equation of state base class.
Definition: eos.h:209
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:777
double * p_ham_constr
Relative error on the Hamiltonian constraint.
Definition: binaire.h:163
double ray_eq() const
Coordinate radius at , [r_unit].
double * p_virial
Virial theorem error.
Definition: binaire.h:154
Base class for coordinate mappings.
Definition: map.h:682
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:780
const Base_vect_cart ref_triad
Cartesian triad of the absolute reference frame.
Definition: binaire.h:115
double virial_gb() const
Estimates the relative error on the virial theorem calculated by E.Gourgoulhon and S...
ostream & operator>>(ostream &) const
Operator >> (function called by the operator <<).
Definition: binaire.C:260
double x_axe
Absolute X coordinate of the rotation axis.
Definition: binaire.h:136
virtual double xa_barycenter() const
Absolute coordinate X of the barycenter of the baryon density, defined according to the formula wher...
double * p_total_ener
Total energy of the system.
Definition: binaire.h:151
double get_gam() const
Returns the adiabatic index (cf. Eq. (3))
Definition: eos_poly.C:271
Etoile_bin * et[2]
Array of the two stars (to perform loops on the stars): { et[0]} contains the address of { star1} and...
Definition: binaire.h:127
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.
void operator=(const Binaire &)
Assignment to another { Binaire}.
Definition: binaire.C:220
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition: binaire.h:132
double separation() const
Returns the coordinate separation of the two stellar centers [{ r_unit}].
Definition: binaire.C:460
double get_kap() const
Returns the pressure coefficient (cf.
Definition: eos_poly.C:275
double * p_virial_fus
Virial theorem error by J.L.Friedman, K.Uryu, and M.Shibata.
Definition: binaire.h:160
Polytropic equation of state (relativistic case).
Definition: eos.h:812
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
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:73
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
const Tbl & angu_mom() const
Total angular momentum.
void write_global(ostream &) const
Write global quantities in a formatted file.
Definition: binaire.C:337
Tbl * p_angu_mom
Total angular momentum of the system.
Definition: binaire.h:148
Binary systems.
Definition: binaire.h:107
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:72
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
Multi-domain grid.
Definition: grilles.h:279
double ray_pole() const
Coordinate radius at [r_unit].
double ray_eq_pi() const
Coordinate radius at , [r_unit].
double mass_kom() const
Total Komar mass.
double virial_fus() const
Estimates the relative error on the virial theorem calculated by J.L.Friedman, K.Uryu, and M.Shibata (PRD accepted, gr-qc/0108070)
double * p_mass_adm
Total ADM mass of the system.
Definition: binaire.h:142
double * p_mass_kom
Total Komar mass of the system.
Definition: binaire.h:145
double * p_virial_gb
Virial theorem error by E.Gourgoulhon and S.Bonazzola.
Definition: binaire.h:157
void del_deriv() const
Destructor.
Definition: binaire.C:180
const Tenseur & get_ener() const
Returns the proper total energy density.
Definition: etoile.h:682
double virial() const
Estimates the relative error on the virial theorem (for a relativistic one, it returns $|1 - M_{ Koma...
double get_ori_z() const
Returns the z coordinate of the origin.
Definition: map.h:784
void set_der_0x0() const
Sets to { 0x0} all the pointers on derived quantities.
Definition: binaire.C:198
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
Binaire(Map &mp1, int nzet1, const Eos &eos1, int irrot1, Map &mp2, int nzet2, const Eos &eos2, int irrot2, int relat)
Standard constructor.
Definition: binaire.C:112
virtual double mass_b() const
Baryon mass.
const Eos & get_eos() const
Returns the equation of state.
Definition: etoile.h:673
const Tenseur & get_ent() const
Returns the enthalpy field.
Definition: etoile.h:676
Etoile_bin star2
Second star of the system.
Definition: binaire.h:121