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