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 
26 char 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 
71 namespace 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 // --------------
200 void 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 // --------
212 ostream& operator<<(ostream& ost, const Binary_xcts& bibi) {
213 
214  bibi >> ost ;
215  return ost ;
216 }
217 
218 
219 ostream& 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 
246 void 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 
296 void 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 
434 double Binary_xcts::separation() const {
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 }
void display_poly(ostream &) const
Display in polytropic units.
Definition: binary_xcts.C:246
ostream & operator>>(ostream &) const
Operator >> (function called by the operator <<).
Definition: binary_xcts.C:219
Map & mp
Mapping associated with the star.
Definition: star.h:180
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
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 get_ori_y() const
Returns the y coordinate of the origin.
Definition: map.h:770
Binary systems in eXtended Conformal Thin Sandwich formulation.
Definition: binary_xcts.h:59
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.
Equation of state base class.
Definition: eos.h:190
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
void write_global(ostream &) const
Write global quantities in a formatted file.
Definition: binary_xcts.C:296
const Scalar & get_ener() const
Returns the proper total energy density.
Definition: star.h:370
Base class for coordinate mappings.
Definition: map.h:670
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:768
double separation() const
Returns the coordinate separation of the two stellar centers [r_unit].
Definition: binary_xcts.C:434
const Map & get_mp() const
Returns the mapping.
Definition: star.h:355
Scalar ent
Log-enthalpy.
Definition: star.h:190
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 get_gam() const
Returns the adiabatic index (cf. Eq. (3))
Definition: eos_poly.C:256
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
Tbl * p_angu_mom
Total angular momentum of the system.
Definition: binary_xcts.h:97
virtual void sauve(FILE *) const
Save in a file.
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.
double x_axe
Absolute X coordinate of the rotation axis.
Definition: binary_xcts.h:84
double get_kap() const
Returns the pressure coefficient (cf.
Definition: eos_poly.C:260
double ray_eq_pis2() const
Coordinate radius at , [r_unit].
Definition: star_global.C:138
Star_bin_xcts star2
Second star of the system.
Definition: binary_xcts.h:69
const Tbl & angu_mom() const
Total angular momentum.
Polytropic equation of state (relativistic case).
Definition: eos.h:757
double mass_adm() const
Total ADM mass.
double * p_mass_kom
Total Komar mass of the system.
Definition: binary_xcts.h:94
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Definition: star_global.C:186
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
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition: binary_xcts.h:80
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
double ray_pole() const
Coordinate radius at [r_unit].
Definition: star_global.C:278
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 get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Multi-domain grid.
Definition: grilles.h:273
double mass_kom_vol() const
Total Komar mass (computed by a volume integral)
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition: binary_xcts.C:162
friend ostream & operator<<(ostream &, const Binary_xcts &)
Display.
Definition: binary_xcts.C:212
~Binary_xcts()
Destructor.
Definition: binary_xcts.C:136
double * p_total_ener
Total energy of the system.
Definition: binary_xcts.h:103
Star_bin_xcts star1
First star of the system.
Definition: binary_xcts.h:66
const Scalar & get_ent() const
Returns the enthalpy field.
Definition: star.h:364
double virial_vol() const
Estimates the relative error on the virial theorem (volume version)
void del_deriv() const
Deletes all the derived quantities.
Definition: binary_xcts.C:146
double ray_eq() const
Coordinate radius at , [r_unit].
Definition: star_global.C:108
virtual double xa_barycenter() const
Absolute coordinate X of the barycenter of the baryon density,.
double * p_virial
Virial theorem error.
Definition: binary_xcts.h:106
virtual double mass_b() const
Baryon mass.
double get_ori_z() const
Returns the z coordinate of the origin.
Definition: map.h:772
double mass_kom() const
Total Komar mass.
const Eos & get_eos() const
Returns the equation of state.
Definition: star.h:361
double mass_adm_vol() const
Total ADM mass (computed by a volume integral)
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
double * p_mass_adm
Total ADM mass of the system.
Definition: binary_xcts.h:91
double * p_virial_vol
Virial theorem error (volume version)
Definition: binary_xcts.h:109
void operator=(const Binary_xcts &)
Assignment to another Binary_xcts.
Definition: binary_xcts.C:182
void sauve(FILE *) const
Save in a file.
Definition: binary_xcts.C:200
Tbl * p_lin_mom
Total linear momentum of the system.
Definition: binary_xcts.h:100
double virial() const
Estimates the relative error on the virial theorem.