LORENE
binary_xcts.C
1/*
2 * Methods of class Binary_xcts
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_xcts_C[] = "$Header: /cvsroot/Lorene/C++/Source/Binary_xcts/binary_xcts.C,v 1.6 2014/10/13 08:52:45 j_novak Exp $" ;
27
28/*
29 * $Id: binary_xcts.C,v 1.6 2014/10/13 08:52:45 j_novak Exp $
30 * $Log: binary_xcts.C,v $
31 * Revision 1.6 2014/10/13 08:52:45 j_novak
32 * Lorene classes and functions now belong to the namespace Lorene.
33 *
34 * Revision 1.5 2014/10/06 15:12:59 j_novak
35 * Modified #include directives to use c++ syntax.
36 *
37 * Revision 1.4 2010/12/20 09:56:02 m_bejger
38 * Pointer to the linear momentum added
39 *
40 * Revision 1.3 2010/12/09 10:38:46 m_bejger
41 * virial_vol() added, fait_decouple() removed
42 *
43 * Revision 1.2 2010/10/28 12:00:07 m_bejger
44 * Mass-shedding indicators added to the output in Binary_xcts::write_global
45 *
46 * Revision 1.1 2010/05/04 07:35:54 m_bejger
47 * Initial version
48 *
49 * $Header: /cvsroot/Lorene/C++/Source/Binary_xcts/binary_xcts.C,v 1.6 2014/10/13 08:52:45 j_novak Exp $
50 *
51 */
52
53// Headers C
54#include <cmath>
55
56// Headers Lorene
57#include "binary_xcts.h"
58#include "eos.h"
59#include "utilitaires.h"
60#include "graphique.h"
61#include "param.h"
62#include "unites.h"
63
64 //--------------//
65 // Constructors //
66 //--------------//
67
68// Standard constructor
69// --------------------
70
71namespace Lorene {
73 int nzet1,
74 const Eos& eos1,
75 int irrot1,
76 Map& mp2,
77 int nzet2,
78 const Eos& eos2,
79 int irrot2)
80 : star1(mp1, nzet1, eos1, irrot1),
81 star2(mp2, nzet2, eos2, irrot2)
82{
83
84 et[0] = &star1 ;
85 et[1] = &star2 ;
86
87 omega = 0 ;
88 x_axe = 0 ;
89
90 // Pointers of derived quantities initialized to zero :
91 set_der_0x0() ;
92}
93
94// Copy constructor
95// ----------------
97 : star1(bibi.star1),
98 star2(bibi.star2),
99 omega(bibi.omega),
100 x_axe(bibi.x_axe)
101{
102 et[0] = &star1 ;
103 et[1] = &star2 ;
104
105 // Pointers of derived quantities initialized to zero :
106 set_der_0x0() ;
107}
108
109// Constructor from a file
110// -----------------------
112 const Eos& eos1,
113 Map& mp2,
114 const Eos& eos2,
115 FILE* fich)
116 : star1(mp1, eos1, fich),
117 star2(mp2, eos2, fich)
118{
119
120 et[0] = &star1 ;
121 et[1] = &star2 ;
122
123 // omega and x_axe are read in the file:
124 fread_be(&omega, sizeof(double), 1, fich) ;
125 fread_be(&x_axe, sizeof(double), 1, fich) ;
126
127 // Pointers of derived quantities initialized to zero :
128 set_der_0x0() ;
129
130}
131
132 //------------//
133 // Destructor //
134 //------- -----//
135
137
138 del_deriv() ;
139
140}
141
142 //----------------------------------//
143 // Management of derived quantities //
144 //----------------------------------//
145
147
148 if (p_mass_adm != 0x0) delete p_mass_adm ;
149 if (p_mass_kom != 0x0) delete p_mass_kom ;
150 if (p_angu_mom != 0x0) delete p_angu_mom ;
151 if (p_lin_mom != 0x0) delete p_lin_mom ;
152 if (p_total_ener != 0x0) delete p_total_ener ;
153 if (p_virial != 0x0) delete p_virial ;
154 if (p_virial_vol != 0x0) delete p_virial_vol ;
155
156 set_der_0x0() ;
157}
158
159
160
161
163
164 p_mass_adm = 0x0 ;
165 p_mass_kom = 0x0 ;
166 p_angu_mom = 0x0 ;
167 p_lin_mom = 0x0 ;
168 p_total_ener = 0x0 ;
169 p_virial = 0x0 ;
170 p_virial_vol = 0x0 ;
171
172}
173
174
175 //--------------//
176 // Assignment //
177 //--------------//
178
179// Assignment to another Binary_xcts
180// --------------------------------
181
183
184 star1 = bibi.star1 ;
185 star2 = bibi.star2 ;
186
187 omega = bibi.omega ;
188 x_axe = bibi.x_axe ;
189
190 del_deriv() ; // Deletes all derived quantities
191
192}
193
194 //--------------//
195 // Outputs //
196 //--------------//
197
198// Save in a file
199// --------------
200void Binary_xcts::sauve(FILE* fich) const {
201
202 star1.sauve(fich) ;
203 star2.sauve(fich) ;
204
205 fwrite_be(&omega, sizeof(double), 1, fich) ;
206 fwrite_be(&x_axe, sizeof(double), 1, fich) ;
207
208}
209
210// Printing
211// --------
212ostream& operator<<(ostream& ost, const Binary_xcts& bibi) {
213
214 bibi >> ost ;
215 return ost ;
216}
217
218
219ostream& Binary_xcts::operator>>(ostream& ost) const {
220
221 using namespace Unites ;
222
223 ost << endl ;
224 ost << "Binary neutron stars" << endl ;
225 ost << "=============" << endl ;
226 ost << endl <<
227 "Orbital angular velocity : " << omega * f_unit << " rad/s" << endl ;
228 ost << endl <<
229 "Coordinate separation between the two stellar centers : "
230 << separation() / km << " km" << endl ;
231 ost <<
232 "Absolute coordinate X of the rotation axis : " << x_axe / km
233 << " km" << endl ;
234 ost << endl << "Star 1 : " << endl ;
235 ost << "====== " << endl ;
236 ost << star1 << endl ;
237 ost << "Star 2 : " << endl ;
238 ost << "====== " << endl ;
239 ost << star2 << endl ;
240 return ost ;
241}
242
243// Display in polytropic units
244// ---------------------------
245
246void Binary_xcts::display_poly(ostream& ost) const {
247
248 using namespace Unites ;
249
250 const Eos* p_eos1 = &( star1.get_eos() ) ;
251 const Eos_poly* p_eos_poly = dynamic_cast<const Eos_poly*>( p_eos1 ) ;
252
253 if (p_eos_poly != 0x0) {
254
255 assert( star1.get_eos() == star2.get_eos() ) ;
256
257 double kappa = p_eos_poly->get_kap() ;
258 double gamma = p_eos_poly->get_gam() ; ;
259 double kap_ns2 = pow( kappa, 0.5 /(gamma-1) ) ;
260
261 // Polytropic unit of length in terms of r_unit :
262 double r_poly = kap_ns2 / sqrt(ggrav) ;
263
264 // Polytropic unit of time in terms of t_unit :
265 double t_poly = r_poly ;
266
267 // Polytropic unit of mass in terms of m_unit :
268 double m_poly = r_poly / ggrav ;
269
270 // Polytropic unit of angular momentum in terms of j_unit :
271 double j_poly = r_poly * r_poly / ggrav ;
272
273 ost.precision(10) ;
274 ost << endl << "Quantities in polytropic units : " << endl ;
275 ost << "==============================" << endl ;
276 ost << " ( r_poly = " << r_poly / km << " km )" << endl ;
277 ost << " d_e_max : " << separation() / r_poly << endl ;
278 ost << " d_G : "
279 << ( star2.xa_barycenter() - star1.xa_barycenter() ) / r_poly
280 << endl ;
281 ost << " Omega : " << omega * t_poly << endl ;
282 ost << " J : " << angu_mom()(2) / j_poly << endl ;
283 ost << " M_ADM : " << mass_adm() / m_poly << endl ;
284 ost << " M_Komar : " << mass_kom() / m_poly << endl ;
285 ost << " M_bar(star 1) : " << star1.mass_b() / m_poly << endl ;
286 ost << " M_bar(star 2) : " << star2.mass_b() / m_poly << endl ;
287 ost << " R_0(star 1) : " <<
288 0.5 * ( star1.ray_eq() + star1.ray_eq_pi() ) / r_poly << endl ;
289 ost << " R_0(star 2) : " <<
290 0.5 * ( star2.ray_eq() + star2.ray_eq_pi() ) / r_poly << endl ;
291
292 }
293
294}
295
296void Binary_xcts::write_global(ostream& ost) const {
297
298 using namespace Unites ;
299
300 const Map& mp1 = star1.get_mp() ;
301 const Mg3d* mg1 = mp1.get_mg() ;
302 int nz1 = mg1->get_nzone() ;
303
304 // Mass-shedding indicators
305 double dent1_eq = (star1.ent).dsdr().val_point(star1.ray_eq(),M_PI/2.,0.) ;
306 double dent1_pole = (star1.ent).dsdr().val_point(star1.ray_pole(),0.,0.) ;
307 double chi1 = fabs( dent1_eq / dent1_pole ) ;
308
309 double dent2_eq = (star2.ent).dsdr().val_point(star2.ray_eq(),M_PI/2.,0.) ;
310 double dent2_pole = (star2.ent).dsdr().val_point(star2.ray_pole(),0.,0.) ;
311 double chi2 = fabs( dent2_eq / dent2_pole ) ;
312
313 ost.precision(5) ;
314 ost << "# Grid 1 : " << nz1 << "x"
315 << mg1->get_nr(0) << "x" << mg1->get_nt(0) << "x" << mg1->get_np(0)
316 << " R_out(l) [km] : " ;
317 for (int l=0; l<nz1; l++) {
318 ost << " " << mp1.val_r(l, 1., M_PI/2, 0) / km ;
319 }
320 ost << endl ;
321
322 ost << "# VE(M) VE(M) [vol]" << endl ;
323
324
325 ost.setf(ios::scientific) ;
326 ost.width(14) ;
327 ost << virial() << " " << virial_vol() << endl ;
328
329 ost << "# d [km] "
330 << " d_G [km] "
331 << " d/(a1 +a1') "
332 << " f [Hz] "
333 << " M_ADM [M_sol] "
334 << " M_ADM_vol [M_sol] "
335 << " M_Komar [M_sol] "
336 << " M_Komar_vol [M_sol] "
337 << " J [G M_sol^2/c] " << endl ;
338
339 ost.precision(14) ;
340 ost.width(20) ;
341 ost << separation() / km ; ost.width(22) ;
342 ost << ( star2.xa_barycenter() - star1.xa_barycenter() ) / km ; ost.width(22) ;
343 ost << separation() / (star1.ray_eq() + star2.ray_eq()) ; ost.width(22) ;
344 ost << omega / (2*M_PI)* f_unit ; ost.width(22) ;
345 ost << mass_adm() / msol ; ost.width(22) ;
346 ost << mass_adm_vol() / msol ; ost.width(22) ;
347 ost << mass_kom() / msol ; ost.width(22) ;
348 ost << mass_kom_vol() / msol ; ost.width(22) ;
349 ost << angu_mom()(2)/ ( qpig / (4* M_PI) * msol*msol) << endl ;
350
351 ost << "# H_c(1)[c^2] "
352 << " e_c(1)[rho_nuc] "
353 << " M_B(1) [M_sol] "
354 << " r_eq(1) [km] "
355 << " a2/a1(1) "
356 << " a3/a1(1) "
357 << " chi1 " << endl ;
358
359 ost.width(20) ;
360 ost << star1.get_ent().val_grid_point(0,0,0,0) ; ost.width(22) ;
361 ost << star1.get_ener().val_grid_point(0,0,0,0) ; ost.width(22) ;
362 ost << star1.mass_b() / msol ; ost.width(22) ;
363 ost << star1.ray_eq() / km ; ost.width(22) ;
364 ost << star1.ray_eq_pis2() / star1.ray_eq() ; ost.width(22) ;
365 ost << star1.ray_pole() / star1.ray_eq() ; ost.width(22) ;
366 ost << chi1 << endl ;
367
368 ost << "# H_c(2)[c^2] "
369 << " e_c(2)[rho_nuc] "
370 << " M_B(2) [M_sol] "
371 << " r_eq(2) [km] "
372 << " a2/a1(2) "
373 << " a3/a1(2) "
374 << " chi2 " << endl ;
375
376
377 ost.width(20) ;
378 ost << star2.get_ent().val_grid_point(0,0,0,0) ; ost.width(22) ;
379 ost << star2.get_ener().val_grid_point(0,0,0,0) ; ost.width(22) ;
380 ost << star2.mass_b() / msol ; ost.width(22) ;
381 ost << star2.ray_eq() / km ; ost.width(22) ;
382 ost << star2.ray_eq_pis2() / star1.ray_eq() ; ost.width(22) ;
383 ost << star2.ray_pole() / star1.ray_eq() ; ost.width(22) ;
384 ost << chi2 << endl ;
385
386 // Quantities in polytropic units if the EOS is a polytropic one
387 // -------------------------------------------------------------
388 const Eos* p_eos1 = &( star1.get_eos() ) ;
389 const Eos_poly* p_eos_poly = dynamic_cast<const Eos_poly*>( p_eos1 ) ;
390
391 if ((p_eos_poly != 0x0) && ( star1.get_eos() == star2.get_eos() )) {
392
393 double kappa = p_eos_poly->get_kap() ;
394 double gamma = p_eos_poly->get_gam() ; ;
395 double kap_ns2 = pow( kappa, 0.5 /(gamma-1.) ) ;
396
397 // Polytropic unit of length in terms of r_unit :
398 double r_poly = kap_ns2 / sqrt(ggrav) ;
399
400 // Polytropic unit of time in terms of t_unit :
401 double t_poly = r_poly ;
402
403 // Polytropic unit of mass in terms of m_unit :
404 double m_poly = r_poly / ggrav ;
405
406 // Polytropic unit of angular momentum in terms of j_unit :
407 double j_poly = r_poly * r_poly / ggrav ;
408
409 ost << "# d [poly] "
410 << " d_G [poly] "
411 << " Omega [poly] "
412 << " M_ADM [poly] "
413 << " J [poly] "
414 << " M_B(1) [poly] "
415 << " M_B(2) [poly] " << endl ;
416
417 ost.width(20) ;
418 ost << separation() / r_poly ; ost.width(22) ;
419 ost << ( star2.xa_barycenter() - star1.xa_barycenter() ) / r_poly ; ost.width(22) ;
420 ost << omega * t_poly ; ost.width(22) ;
421 ost << mass_adm() / m_poly ; ost.width(22) ;
422 ost << angu_mom()(2) / j_poly ; ost.width(22) ;
423 ost << star1.mass_b() / m_poly ; ost.width(22) ;
424 ost << star2.mass_b() / m_poly << endl ;
425
426 }
427
428}
429
430 //-------------------------------//
431 // Miscellaneous //
432 //-------------------------------//
433
435
436 double dx = star1.mp.get_ori_x() - star2.mp.get_ori_x() ;
437 double dy = star1.mp.get_ori_y() - star2.mp.get_ori_y() ;
438 double dz = star1.mp.get_ori_z() - star2.mp.get_ori_z() ;
439
440 return sqrt( dx*dx + dy*dy + dz*dz ) ;
441
442}
443}
Binary systems in eXtended Conformal Thin Sandwich formulation.
Definition binary_xcts.h:59
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
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
double x_axe
Absolute X coordinate of the rotation axis.
Definition binary_xcts.h:84
void write_global(ostream &) const
Write global quantities in a formatted file.
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.
void operator=(const Binary_xcts &)
Assignment to another Binary_xcts.
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
void sauve(FILE *) const
Save in a file.
double mass_adm() const
Total ADM mass.
double separation() const
Returns the coordinate separation of the two stellar centers [r_unit].
void del_deriv() const
Deletes all the derived quantities.
Binary_xcts(Map &mp1, int nzet1, const Eos &eos1, int irrot1, Map &mp2, int nzet2, const Eos &eos2, int irrot2)
Standard constructor.
Definition binary_xcts.C:72
double * p_mass_adm
Total ADM mass of the system.
Definition binary_xcts.h:91
double * p_total_ener
Total energy of the system.
void display_poly(ostream &) const
Display in polytropic units.
ostream & operator>>(ostream &) const
Operator >> (function called by the operator <<).
double mass_kom() const
Total Komar mass.
~Binary_xcts()
Destructor.
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.
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition binary_xcts.h:80
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
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
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition scalar.h:637
virtual double xa_barycenter() const
Absolute coordinate X of the barycenter of the baryon density,.
virtual void sauve(FILE *) const
Save in a file.
virtual double mass_b() const
Baryon mass.
const Map & get_mp() const
Returns the mapping.
Definition star.h:355
const Eos & get_eos() const
Returns the equation of state.
Definition star.h:361
const Scalar & get_ener() const
Returns the proper total energy density.
Definition star.h:370
double ray_eq() const
Coordinate radius at , [r_unit].
double ray_eq_pis2() const
Coordinate radius at , [r_unit].
const Scalar & get_ent() const
Returns the enthalpy field.
Definition star.h:364
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Scalar ent
Log-enthalpy.
Definition star.h:190
Map & mp
Mapping associated with the star.
Definition star.h:180
double ray_pole() const
Coordinate radius at [r_unit].
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.