LORENE
binary_orbit_xcts.C
1 /*
2  * Method of class Binary_xcts to compute the orbital angular velocity
3  * {\tt omega} and the position of the rotation axis {\tt x_axe}.
4  * (See file binary_xcts.h for documentation)
5  */
6 
7 /*
8  * Copyright (c) 2010 Michal Bejger
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 version 2
14  * as published by the Free Software Foundation.
15  *
16  * LORENE is distributed in the hope that it will be useful,
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19  * GNU General Public License for more details.
20  *
21  * You should have received a copy of the GNU General Public License
22  * along with LORENE; if not, write to the Free Software
23  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24  *
25  */
26 
27 
28 
29 /*
30  * $Id: binary_orbit_xcts.C,v 1.15 2016/12/05 16:17:47 j_novak Exp $
31  * $Log: binary_orbit_xcts.C,v $
32  * Revision 1.15 2016/12/05 16:17:47 j_novak
33  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
34  *
35  * Revision 1.14 2014/10/24 14:10:24 j_novak
36  * Minor change to prevent weird error from g++-4.8...
37  *
38  * Revision 1.13 2014/10/13 08:52:45 j_novak
39  * Lorene classes and functions now belong to the namespace Lorene.
40  *
41  * Revision 1.12 2014/10/06 15:12:59 j_novak
42  * Modified #include directives to use c++ syntax.
43  *
44  * Revision 1.11 2011/03/30 13:14:27 m_bejger
45  * Psi and chi rewritten using auto and comp parts to improve the convergence (in all the remaining fields, not only logn)
46  *
47  * Revision 1.10 2011/03/28 17:13:37 m_bejger
48  * logn in dnulg stated using Psi1,2 and chi1,2
49  *
50  * Revision 1.9 2011/03/27 16:41:19 e_gourgoulhon
51  * -- Corrected sign of ny and dny for star no. 2
52  * -- Added output via new function save_profile for graphics
53  *
54  * Revision 1.8 2011/03/27 14:58:48 m_bejger
55  * dnulg by means of dsdx(); rearrangements to use primary variables
56  *
57  * Revision 1.7 2011/03/25 16:28:36 e_gourgoulhon
58  * Still in progress
59  *
60  * Revision 1.6 2010/12/09 10:41:20 m_bejger
61  * For testing; not sure if working properly
62  *
63  * Revision 1.5 2010/10/26 19:45:45 m_bejger
64  * Cleanup
65  *
66  * Revision 1.4 2010/07/16 16:27:19 m_bejger
67  * This version is basically a copy of the one used by Binaire (binaire_orbite.C)
68  *
69  * Revision 1.3 2010/06/17 14:15:41 m_bejger
70  * Using method get_Psi()
71  *
72  * Revision 1.2 2010/06/15 07:57:30 m_bejger
73  * Minor corrections
74  *
75  * Revision 1.1 2010/05/04 07:35:54 m_bejger
76  * Initial version
77  *
78  * $Header: /cvsroot/Lorene/C++/Source/Binary_xcts/binary_orbit_xcts.C,v 1.15 2016/12/05 16:17:47 j_novak Exp $
79  *
80  */
81 
82 // Headers C
83 #include <cmath>
84 
85 // Headers Lorene
86 #include "binary_xcts.h"
87 #include "eos.h"
88 #include "param.h"
89 #include "utilitaires.h"
90 #include "unites.h"
91 
92 #include "graphique.h"
93 
94 namespace Lorene {
95 double fonc_binary_xcts_axe(double , const Param& ) ;
96 double fonc_binary_xcts_orbit(double , const Param& ) ;
97 
98 //******************************************************************************
99 
100 void Binary_xcts::orbit(double fact_omeg_min,
101  double fact_omeg_max,
102  double& xgg1,
103  double& xgg2) {
104 
105  using namespace Unites ;
106 
107  //-------------------------------------------------------------
108  // Evaluation of various quantities at the center of each star
109  //-------------------------------------------------------------
110 
111  double dnulg[2], asn2[2], dasn2[2], ny[2], dny[2], npn[2], dnpn[2], xgg[2] ;
112  double nyso[2], dnyso[2], npnso2[2], dnpnso2[2], ori_x[2] ;
113 
114  for (int i=0; i<2; i++) {
115 
116  const Map& mp = et[i]->get_mp() ;
117  const Metric& flat = et[i]->get_flat() ;
118 
119  //------------------------------------------------------------------
120  // Recasting Phi and chi to manifestly equal auto and comp part
121  // - more fortunate from the point of view of Omega computation
122  //------------------------------------------------------------------
123  Scalar chi = et[i]->get_chi_auto() + et[i]->get_chi_comp() + 1 ;
124  Scalar Psi = et[i]->get_Psi_auto() + et[i]->get_Psi_comp() + 1 ;
125 
126  Scalar logn = log(chi) - log(Psi) ;
127  logn.std_spectral_base() ;
128 
129  // Sign convention for shift (beta^i = - N^i)
130  Vector shift = - ( et[i]->get_beta() ) ;
131  shift.change_triad(et[i]->mp.get_bvect_cart()) ;
132 
133  //------------------------------------------------------------------
134  // d/dX( log(N) + log(Gamma) )
135  // in the center of the star ---> dnulg[i]
136  //------------------------------------------------------------------
137 
138  // Facteur de passage x --> X :
139  double factx ;
140 
141  if (fabs(mp.get_rot_phi()) < 1.e-14) factx = 1. ;
142  else {
143 
144  if (fabs(mp.get_rot_phi() - M_PI) < 1.e-14) {
145  factx = - 1. ;
146 
147  } else {
148 
149  cout << "Binary_xcts::orbit : unknown value of rot_phi !" << endl ;
150  abort() ;
151  }
152  }
153 
154  Scalar tmp = logn + et[i]->get_loggam() ;
155  dnulg[i] = factx*tmp.dsdx().val_grid_point(0, 0, 0, 0) ;
156 
157  // For graphical outputs:
158  Scalar tgraph = logn - log( (1. + et[i]->get_chi_auto()) / (1. + et[i]->get_Psi_auto()) ) ;
159  // tmp = log( (1. + et[i]->get_chi_comp()) / (1. + et[i]->get_Psi_comp()) ) ;
160  tgraph.std_spectral_base() ;
161  save_profile(tgraph, 0., 10., 0.5*M_PI, 0., "prof_logn.d") ;
162  save_profile(et[i]->get_loggam(), 0., 1.8, 0.5*M_PI, 0., "prof_loggam.d") ;
163 
164  //------------------------------------------------------------------
165  // Psi^4/N^2 = in the center of the star ---> asn2[i]
166  //------------------------------------------------------------------
167 
168  Scalar Psi6schi2 = pow(Psi, 6)/(chi % chi) ;
169  Psi6schi2.std_spectral_base() ;
170  asn2[i] = Psi6schi2.val_grid_point(0, 0, 0, 0) ;
171 
172  //------------------------------------------------------------------
173  // d/dX(A^2/N^2) in the center of the star ---> dasn2[i]
174  //------------------------------------------------------------------
175 
176  dasn2[i] = Psi6schi2.dsdx().val_grid_point(0, 0, 0, 0) ;
177 
178  //------------------------------------------------------------------
179  // N^Y in the center of the star ---> ny[i]
180  //------------------------------------------------------------------
181 
182  ny[i] = factx*shift(2).val_grid_point(0, 0, 0, 0) ;
183 
184  nyso[i] = ny[i] / omega ;
185 
186  //------------------------------------------------------------------
187  // dN^Y/dX in the center of the star ---> dny[i]
188  //------------------------------------------------------------------
189 
190  dny[i] = shift(2).dsdx().val_grid_point(0, 0, 0, 0) ;
191 
192  dnyso[i] = dny[i] / omega ;
193 
194  //------------------------------------------------------------------
195  // (N^X)^2 + (N^Y)^2 + (N^Z)^2
196  // in the center of the star ---> npn[i]
197  //------------------------------------------------------------------
198 
199  tmp = contract(shift, 0, shift.up_down(flat), 0) ;
200 
201  npn[i] = tmp.val_grid_point(0, 0, 0, 0) ;
202  npnso2[i] = npn[i] / omega / omega ;
203 
204  //------------------------------------------------------------------
205  // d/dX( (N^X)^2 + (N^Y)^2 + (N^Z)^2 )
206  // in the center of the star ---> dnpn[i]
207  //------------------------------------------------------------------
208 
209  dnpn[i] = factx * tmp.dsdx().val_grid_point(0, 0, 0, 0) ;
210  dnpnso2[i] = dnpn[i] / omega / omega ;
211 
212  cout << "Binary_xcts::orbit: central d(nu+log(Gam))/dX : "
213  << dnulg[i] << endl ;
214  cout << "Binary_xcts::orbit: central A^2/N^2 : "
215  << asn2[i] << endl ;
216  cout << "Binary_xcts::orbit: central d(A^2/N^2)/dX : "
217  << dasn2[i] << endl ;
218  cout << "Binary_xcts::orbit: central N^Y : "
219  << ny[i] << endl ;
220  cout << "Binary_xcts::orbit: central dN^Y/dX : "
221  << dny[i] << endl ;
222  cout << "Binary_xcts::orbit: central N.N : "
223  << npn[i] << endl ;
224  cout << "Binary_xcts::orbit: central d(N.N)/dX : "
225  << dnpn[i] << endl ;
226 
227 
228  ori_x[i] = (et[i]->get_mp()).get_ori_x() ;
229  xgg[i] = factx * (et[i]->xa_barycenter() - ori_x[i]) ;
230 
231  }
232 
233  xgg1 = xgg[0] ;
234  xgg2 = xgg[1] ;
235 
236 //---------------------------------
237 // axis of rotation
238 //---------------------------------
239 
240  int relat = 1 ;
241 
242  double ori_x1 = ori_x[0] ;
243  double ori_x2 = ori_x[1] ;
244 
245  if ( et[0]->get_eos() == et[1]->get_eos() &&
246  fabs( et[0]->get_ent().val_grid_point(0,0,0,0) -
247  et[1]->get_ent().val_grid_point(0,0,0,0) ) < 1.e-14 ) {
248 
249  x_axe = 0. ;
250 
251  } else {
252 
253  Param paraxe ;
254  paraxe.add_int(relat) ;
255  paraxe.add_double( ori_x1, 0) ;
256  paraxe.add_double( ori_x2, 1) ;
257  paraxe.add_double( dnulg[0], 2) ;
258  paraxe.add_double( dnulg[1], 3) ;
259  paraxe.add_double( asn2[0], 4) ;
260  paraxe.add_double( asn2[1], 5) ;
261  paraxe.add_double( dasn2[0], 6) ;
262  paraxe.add_double( dasn2[1], 7) ;
263  paraxe.add_double( nyso[0], 8) ;
264  paraxe.add_double( nyso[1], 9) ;
265  paraxe.add_double( dnyso[0], 10) ;
266  paraxe.add_double( dnyso[1], 11) ;
267  paraxe.add_double( npnso2[0], 12) ;
268  paraxe.add_double( npnso2[1], 13) ;
269  paraxe.add_double( dnpnso2[0], 14) ;
270  paraxe.add_double( dnpnso2[1], 15) ;
271 
272  int nitmax_axe = 200 ;
273  int nit_axe ;
274  double precis_axe = 1.e-13 ;
275 
276  x_axe = zerosec(fonc_binary_xcts_axe, paraxe, 0.9*ori_x1, 0.9*ori_x2,
277  precis_axe, nitmax_axe, nit_axe) ;
278 
279  cout << "Binary_xcts::orbit : Number of iterations in zerosec for x_axe : "
280  << nit_axe << endl ;
281  }
282 
283  cout << "Binary_xcts::orbit: x_axe [km] : " << x_axe / km << endl ;
284 
285 //-------------------------------------
286 // Calcul de la vitesse orbitale
287 //-------------------------------------
288 
289  Param parf ;
290  parf.add_int(relat) ;
291  parf.add_double( ori_x1, 0) ;
292  parf.add_double( dnulg[0], 1) ;
293  parf.add_double( asn2[0], 2) ;
294  parf.add_double( dasn2[0], 3) ;
295  parf.add_double( ny[0], 4) ;
296  parf.add_double( dny[0], 5) ;
297  parf.add_double( npn[0], 6) ;
298  parf.add_double( dnpn[0], 7) ;
299  parf.add_double( x_axe, 8) ;
300 
301  double omega1 = fact_omeg_min * omega ;
302  double omega2 = fact_omeg_max * omega ;
303  cout << "Binary_xcts::orbit: omega1, omega2 [rad/s] : "
304  << omega1 * f_unit << " " << omega2 * f_unit << endl ;
305 
306  // Search for the various zeros in the interval [omega1, omega2]
307  // ------------------------------------------------------------
308  int nsub = 50 ; // total number of subdivisions of the interval
309  Tbl* azer = 0x0 ;
310  Tbl* bzer = 0x0 ;
311  zero_list(fonc_binary_xcts_orbit, parf, omega1, omega2, nsub,
312  azer, bzer) ;
313 
314  // Search for the zero closest to the previous value of omega
315  // ----------------------------------------------------------
316  double omeg_min, omeg_max ;
317  int nzer = azer->get_taille() ; // number of zeros found by zero_list
318  cout << "Binary_xcts:orbit : " << nzer <<
319  " zero(s) found in the interval [omega1, omega2]." << endl ;
320  cout << "omega, omega1, omega2 : " << omega << " " << omega1
321  << " " << omega2 << endl ;
322  cout << "azer : " << *azer << endl ;
323  cout << "bzer : " << *bzer << endl ;
324 
325  if (nzer == 0) {
326  cout <<
327  "Binary_xcts::orbit: WARNING : no zero detected in the interval"
328  << endl << " [" << omega1 * f_unit << ", "
329  << omega2 * f_unit << "] rad/s !" << endl ;
330  omeg_min = omega1 ;
331  omeg_max = omega2 ;
332  }
333  else {
334  double dist_min = fabs(omega2 - omega1) ;
335  int i_dist_min = -1 ;
336  for (int i=0; i<nzer; i++) {
337  // Distance of previous value of omega from the center of the
338  // interval [azer(i), bzer(i)]
339  double dist = fabs( omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
340  if (dist < dist_min) {
341  dist_min = dist ;
342  i_dist_min = i ;
343  }
344  }
345  omeg_min = (*azer)(i_dist_min) ;
346  omeg_max = (*bzer)(i_dist_min) ;
347  }
348 
349  delete azer ; // Tbl allocated by zero_list
350  delete bzer ; //
351 
352  cout << "Binary_xcts::orbit : interval selected for the search of the zero : "
353  << endl << " [" << omeg_min << ", " << omeg_max << "] = ["
354  << omeg_min * f_unit << ", " << omeg_max * f_unit << "] rad/s " << endl ;
355 
356  // Computation of the zero in the selected interval by the secant method
357  // ---------------------------------------------------------------------
358  int nitermax = 200 ;
359  int niter ;
360  double precis = 1.e-13 ;
361  omega = zerosec_b(fonc_binary_xcts_orbit, parf, omeg_min, omeg_max,
362  precis, nitermax, niter) ;
363 
364  cout << "Binary_xcts::orbit : Number of iterations in zerosec for omega : "
365  << niter << endl ;
366 
367  cout << "Binary_xcts::orbit : omega [rad/s] : "
368  << omega * f_unit << endl ;
369 
370 }
371 
372 //*************************************************
373 // Function used for search of the rotation axis
374 //*************************************************
375 
376 double fonc_binary_xcts_axe(double x_rot, const Param& paraxe) {
377 
378  int relat = paraxe.get_int() ;
379 
380  double ori_x1 = paraxe.get_double(0) ;
381  double ori_x2 = paraxe.get_double(1) ;
382  double dnulg_1 = paraxe.get_double(2) ;
383  double dnulg_2 = paraxe.get_double(3) ;
384  double asn2_1 = paraxe.get_double(4) ;
385  double asn2_2 = paraxe.get_double(5) ;
386  double dasn2_1 = paraxe.get_double(6) ;
387  double dasn2_2 = paraxe.get_double(7) ;
388  double nyso_1 = paraxe.get_double(8) ;
389  double nyso_2 = paraxe.get_double(9) ;
390  double dnyso_1 = paraxe.get_double(10) ;
391  double dnyso_2 = paraxe.get_double(11) ;
392  double npnso2_1 = paraxe.get_double(12) ;
393  double npnso2_2 = paraxe.get_double(13) ;
394  double dnpnso2_1 = paraxe.get_double(14) ;
395  double dnpnso2_2 = paraxe.get_double(15) ;
396 
397  double om2_star1 ;
398  double om2_star2 ;
399 
400  double x1 = ori_x1 - x_rot ;
401  double x2 = ori_x2 - x_rot ;
402 
403  if (relat == 1) {
404 
405  double andan_1 = 0.5 * dasn2_1 + asn2_1 * dnulg_1 ;
406  double andan_2 = 0.5 * dasn2_2 + asn2_2 * dnulg_2 ;
407 
408  double bpb_1 = x1 * x1 - 2. * nyso_1 * x1 + npnso2_1 ;
409  double bpb_2 = x2 * x2 - 2. * nyso_2 * x2 + npnso2_2 ;
410 
411  double cpc_1 = 0.5 * dnpnso2_1 + x1 * (1. - dnyso_1) - nyso_1 ;
412  double cpc_2 = 0.5 * dnpnso2_2 + x2 * (1. - dnyso_2) - nyso_2 ;
413 
414  om2_star1 = dnulg_1 / (andan_1 * bpb_1 + asn2_1 * cpc_1) ;
415  om2_star2 = dnulg_2 / (andan_2 * bpb_2 + asn2_2 * cpc_2) ;
416 
417  }
418  else {
419 
420  om2_star1 = dnulg_1 / x1 ;
421  om2_star2 = dnulg_2 / x2 ;
422 
423  }
424 
425  return om2_star1 - om2_star2 ;
426 
427 }
428 
429 //*****************************************************************************
430 // Fonction utilisee pour la recherche de omega par la methode de la secante
431 //*****************************************************************************
432 
433 double fonc_binary_xcts_orbit(double om, const Param& parf) {
434 
435  int relat = parf.get_int() ;
436 
437  double xc = parf.get_double(0) ;
438  double dnulg = parf.get_double(1) ;
439  double asn2 = parf.get_double(2) ;
440  double dasn2 = parf.get_double(3) ;
441  double ny = parf.get_double(4) ;
442  double dny = parf.get_double(5) ;
443  double npn = parf.get_double(6) ;
444  double dnpn = parf.get_double(7) ;
445  double x_axe = parf.get_double(8) ;
446 
447  double xx = xc - x_axe ;
448  double om2 = om*om ;
449 
450  double dphi_cent ;
451 
452  if (relat == 1) {
453  double bpb = om2 * xx*xx - 2*om * ny * xx + npn ;
454 
455  dphi_cent = ( asn2* ( om* (ny + xx*dny) - om2*xx - 0.5*dnpn )
456  - 0.5*bpb* dasn2 )
457  / ( 1 - asn2 * bpb ) ;
458  }
459  else {
460  dphi_cent = - om2*xx ;
461  }
462 
463  return dnulg + dphi_cent ;
464 
465 }
466 
467 
468 }
Metric for tensor calculation.
Definition: metric.h:90
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:249
const Vector & get_beta() const
Returns the shift vector .
Definition: star.h:402
const Scalar & get_loggam() const
Returns the logarithm of the Lorentz factor between the fluid and the co-orbiting observer...
Definition: star.h:1336
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
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
void orbit(double fact_omeg_min, double fact_omeg_max, double &xgg1, double &xgg2)
Computes the orbital angular velocity omega and the position of the rotation axis x_axe...
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Base class for coordinate mappings.
Definition: map.h:688
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
double zerosec_b(double(*f)(double, const Param &), const Param &par, double a, double b, double precis, int nitermax, int &niter)
Finding the zero a function on a bounded domain.
Definition: zerosec_borne.C:71
const Scalar & dsdx() const
Returns of *this , where .
Definition: scalar_deriv.C:266
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
const Map & get_mp() const
Returns the mapping.
Definition: star.h:355
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition: map.h:793
Tensor field of valence 1.
Definition: vector.h:188
double zerosec(double(*f)(double, const Param &), const Param &par, double a, double b, double precis, int nitermax, int &niter, bool abort=true)
Finding the zero a function.
Definition: zerosec.C:92
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:643
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:295
double x_axe
Absolute X coordinate of the rotation axis.
Definition: binary_xcts.h:84
const Scalar & get_Psi_auto() const
Returns the scalar field generated principally by the star.
Definition: star.h:1362
void zero_list(double(*f)(double, const Param &), const Param &par, double xmin, double xmax, int nsub, Tbl *&az, Tbl *&bz)
Locates approximatively all the zeros of a function in a given interval.
Definition: zero_list.C:60
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
Parameter storage.
Definition: param.h:125
const Scalar & get_chi_comp() const
Returns the scalar field generated principally by the companion star.
Definition: star.h:1377
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:351
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
const Scalar & get_Psi_comp() const
Returns the scalar field generated principally by the companion star.
Definition: star.h:1367
void save_profile(const Scalar &uu, double r_min, double r_max, double theta, double phi, const char *filename)
Saves in a file the profile of a Scalar along some radial axis determined by a fixed value of ...
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:809
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:417
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:318
virtual double xa_barycenter() const
Absolute coordinate X of the barycenter of the baryon density,.
const Eos & get_eos() const
Returns the equation of state.
Definition: star.h:361
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition: param.C:364
Basic array class.
Definition: tbl.h:164
const Scalar & get_chi_auto() const
Returns the scalar field generated principally by the star.
Definition: star.h:1372