LORENE
binary_orbite.C
1 /*
2  * Method of class Binary to compute the orbital angular velocity {\tt omega}
3  * and the position of the rotation axis {\tt x_axe}.
4  *
5  * (See file binary.h for documentation)
6  *
7  */
8 
9 /*
10  * Copyright (c) 2004 Francois Limousin
11  *
12  * This file is part of LORENE.
13  *
14  * LORENE is free software; you can redistribute it and/or modify
15  * it under the terms of the GNU General Public License as published by
16  * the Free Software Foundation; either version 2 of the License, or
17  * (at your option) any later version.
18  *
19  * LORENE is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU General Public License for more details.
23  *
24  * You should have received a copy of the GNU General Public License
25  * along with LORENE; if not, write to the Free Software
26  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27  *
28  */
29 
30 
31 
32 
33 /*
34  * $Id: binary_orbite.C,v 1.10 2016/12/05 16:17:47 j_novak Exp $
35  * $Log: binary_orbite.C,v $
36  * Revision 1.10 2016/12/05 16:17:47 j_novak
37  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
38  *
39  * Revision 1.9 2015/08/10 15:32:26 j_novak
40  * Better calls to Param::add_int(), to avoid weird problems (e.g. with g++ 4.8).
41  *
42  * Revision 1.8 2014/10/13 08:52:45 j_novak
43  * Lorene classes and functions now belong to the namespace Lorene.
44  *
45  * Revision 1.7 2014/10/06 15:12:59 j_novak
46  * Modified #include directives to use c++ syntax.
47  *
48  * Revision 1.6 2005/09/13 19:38:31 f_limousin
49  * Reintroduction of the resolution of the equations in cartesian coordinates.
50  *
51  * Revision 1.5 2005/02/17 17:35:35 f_limousin
52  * Change the name of some quantities to be consistent with other classes
53  * (for instance nnn is changed to nn, shift to beta, beta to lnq...)
54  *
55  * Revision 1.4 2004/03/25 10:29:02 j_novak
56  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
57  *
58  * Revision 1.3 2004/02/24 12:39:30 f_limousin
59  * Change fonc_bin_ncp_orbit to fonc_binary_orbit and fonc_bin_ncp_axe
60  * to fonc_binary_axe.
61  *
62  * Revision 1.2 2004/02/21 17:05:13 e_gourgoulhon
63  * Method Scalar::point renamed Scalar::val_grid_point.
64  * Method Scalar::set_point renamed Scalar::set_grid_point.
65  *
66  * Revision 1.1 2004/01/20 15:22:19 f_limousin
67  * First version
68  *
69  *
70  * $Header: /cvsroot/Lorene/C++/Source/Binary/binary_orbite.C,v 1.10 2016/12/05 16:17:47 j_novak Exp $
71  *
72  */
73 
74 // Headers C
75 #include <cmath>
76 
77 // Headers Lorene
78 #include "binary.h"
79 #include "eos.h"
80 #include "param.h"
81 #include "utilitaires.h"
82 #include "unites.h"
83 
84 namespace Lorene {
85 double fonc_binary_axe(double , const Param& ) ;
86 double fonc_binary_orbit(double , const Param& ) ;
87 
88 //******************************************************************************
89 
90 void Binary::orbit(double fact_omeg_min, double fact_omeg_max, double& xgg1,
91  double& xgg2) {
92 
93 using namespace Unites ;
94 
95  //-------------------------------------------------------------
96  // Evaluation of various quantities at the center of each star
97  //-------------------------------------------------------------
98 
99 
100  double g00[2], g10[2], g20[2], g11[2], g21[2], g22[2], bx[2], by[2] ;
101 
102  double bz[2], d1sn2[2], unsn2[2] ;
103 
104  double dnulg[2], xgg[2], ori_x[2], dg00[2], dg10[2], dg20[2], dg11[2] ;
105 
106  double dg21[2], dg22[2], dbx[2], dby[2], dbz[2], dbymo[2] ;
107 
108  for (int i=0; i<2; i++) {
109 
110  const Scalar& logn_auto = et[i]->get_logn_auto() ;
111  const Scalar& logn_comp = et[i]->get_logn_comp() ;
112  const Scalar& loggam = et[i]->get_loggam() ;
113  const Scalar& nn = et[i]->get_nn() ;
114  Vector shift = et[i]->get_beta() ;
115  const Metric& gamma = et[i]->get_gamma() ;
116 
117  Tensor gamma_cov = gamma.cov() ;
118 
119  // With the new convention for shift (beta^i = - N^i)
120  shift = - shift ;
121 
122  // All tensors must be in the cartesian triad
123 
124  shift.change_triad(et[i]->mp.get_bvect_cart()) ;
125  gamma_cov.change_triad(et[i]->mp.get_bvect_cart()) ;
126 
127  const Scalar& gg00 = gamma_cov(1,1) ;
128  const Scalar& gg10 = gamma_cov(2,1) ;
129  const Scalar& gg20 = gamma_cov(3,1) ;
130  const Scalar& gg11 = gamma_cov(2,2) ;
131  const Scalar& gg21 = gamma_cov(3,2) ;
132  const Scalar& gg22 = gamma_cov(3,3) ;
133 
134  const Scalar& bbx = shift(1) ;
135  const Scalar& bby = shift(2) ;
136  const Scalar& bbz = shift(3) ;
137 
138  cout << "gg00" << endl << norme(gg00) << endl ;
139  cout << "gg10" << endl << norme(gg10) << endl ;
140  cout << "gg20" << endl << norme(gg20) << endl ;
141  cout << "gg11" << endl << norme(gg11) << endl ;
142  cout << "gg21" << endl << norme(gg21) << endl ;
143  cout << "gg22" << endl << norme(gg22) << endl ;
144 
145  cout << "bbx" << endl << norme(bbx) << endl ;
146  cout << "bby" << endl << norme(bby) << endl ;
147  cout << "bbz" << endl << norme(bbz) << endl ;
148 
149  //----------------------------------
150  // Calcul de d/dX( nu + ln(Gamma) ) au centre de l'etoile ---> dnulg[i]
151  //----------------------------------
152 
153  Scalar tmp = logn_auto + logn_comp + loggam ;
154 
155  cout << "logn" << endl << norme(logn_auto + logn_comp) << endl ;
156  cout << "loggam" << endl << norme(loggam) << endl ;
157  cout << "dnulg" << endl << norme(tmp.dsdx()) << endl ;
158 
159  // ... gradient suivant X :
160  dnulg[i] = tmp.dsdx().val_grid_point(0, 0, 0, 0) ;
161 
162  cout.precision(8) ;
163  cout << "dnulg = " << dnulg[i] << endl ;
164 
165 
166  //----------------------------------
167  // Calcul de gij, lapse et shift au centre de l'etoile
168  //----------------------------------
169 
170  g00[i] = gg00.val_grid_point(0,0,0,0) ;
171  g10[i] = gg10.val_grid_point(0,0,0,0) ;
172  g20[i] = gg20.val_grid_point(0,0,0,0) ;
173  g11[i] = gg11.val_grid_point(0,0,0,0) ;
174  g21[i] = gg21.val_grid_point(0,0,0,0) ;
175  g22[i] = gg22.val_grid_point(0,0,0,0) ;
176 
177  bx[i] = bbx.val_grid_point(0,0,0,0) ;
178  by[i] = bby.val_grid_point(0,0,0,0) ;
179  bz[i] = bbz.val_grid_point(0,0,0,0) ;
180 
181  unsn2[i] = 1/(nn.val_grid_point(0,0,0,0)*nn.val_grid_point(0,0,0,0)) ;
182 
183  //----------------------------------
184  // Calcul de d/dX(gij), d/dX(shift) au centre de l'etoile
185  //----------------------------------
186 
187  dg00[i] = gg00.dsdx().val_grid_point(0,0,0,0) ;
188  dg10[i] = gg10.dsdx().val_grid_point(0,0,0,0) ;
189  dg20[i] = gg20.dsdx().val_grid_point(0,0,0,0) ;
190  dg11[i] = gg11.dsdx().val_grid_point(0,0,0,0) ;
191  dg21[i] = gg21.dsdx().val_grid_point(0,0,0,0) ;
192  dg22[i] = gg22.dsdx().val_grid_point(0,0,0,0) ;
193 
194  dbx[i] = bbx.dsdx().val_grid_point(0,0,0,0) ;
195  dby[i] = bby.dsdx().val_grid_point(0,0,0,0) ;
196  dbz[i] = bbz.dsdx().val_grid_point(0,0,0,0) ;
197 
198  dbymo[i] = bby.dsdx().val_grid_point(0,0,0,0) - omega ;
199 
200 
201  d1sn2[i] = (1/(nn*nn)).dsdx().val_grid_point(0,0,0,0) ;
202 
203 
204  cout << "Binary::orbit: central d(nu+log(Gam))/dX : "
205  << dnulg[i] << endl ;
206  cout << "Binary::orbit: central g00 :" << g00[i] << endl ;
207  cout << "Binary::orbit: central g10 :" << g10[i] << endl ;
208  cout << "Binary::orbit: central g20 :" << g20[i] << endl ;
209  cout << "Binary::orbit: central g11 :" << g11[i] << endl ;
210  cout << "Binary::orbit: central g21 :" << g21[i] << endl ;
211  cout << "Binary::orbit: central g22 :" << g22[i] << endl ;
212 
213  cout << "Binary::orbit: central shift_x :" << bx[i] << endl ;
214  cout << "Binary::orbit: central shift_y :" << by[i] << endl ;
215  cout << "Binary::orbit: central shift_z :" << bz[i] << endl ;
216 
217  cout << "Binary::orbit: central d/dX(g00) :" << dg00[i] << endl ;
218  cout << "Binary::orbit: central d/dX(g10) :" << dg10[i] << endl ;
219  cout << "Binary::orbit: central d/dX(g20) :" << dg20[i] << endl ;
220  cout << "Binary::orbit: central d/dX(g11) :" << dg11[i] << endl ;
221  cout << "Binary::orbit: central d/dX(g21) :" << dg21[i] << endl ;
222  cout << "Binary::orbit: central d/dX(g22) :" << dg22[i] << endl ;
223 
224  cout << "Binary::orbit: central d/dX(shift_x) :" << dbx[i] << endl ;
225  cout << "Binary::orbit: central d/dX(shift_y) :" << dby[i] << endl ;
226  cout << "Binary::orbit: central d/dX(shift_z) :" << dbz[i] << endl ;
227 
228  //----------------------
229  // Pour information seulement : 1/ calcul des positions des "centres de
230  // de masse"
231  // 2/ calcul de dH/dX en r=0
232  //-----------------------
233 
234  ori_x[i] = (et[i]->get_mp()).get_ori_x() ;
235 
236  xgg[i] = (et[i]->xa_barycenter() - ori_x[i]) ;
237 
238  } // fin de la boucle sur les etoiles
239 
240  xgg1 = xgg[0] ;
241  xgg2 = xgg[1] ;
242 
243  //---------------------------------
244  // Position de l'axe de rotation
245  //---------------------------------
246 
247  double ori_x1 = ori_x[0] ;
248  double ori_x2 = ori_x[1] ;
249 
250  if ( et[0]->get_eos() == et[1]->get_eos() &&
251  et[0]->get_ent().val_grid_point(0,0,0,0) ==
252  et[1]->get_ent().val_grid_point(0,0,0,0) ) {
253 
254  x_axe = 0. ;
255 
256  }
257  else {
258 
259  Param paraxe ;
260  paraxe.add_double( ori_x1, 0) ;
261  paraxe.add_double( ori_x2, 1) ;
262  paraxe.add_double( dnulg[0], 2) ;
263  paraxe.add_double( dnulg[1], 3) ;
264  paraxe.add_double( g00[0], 4) ;
265  paraxe.add_double( g00[1], 5) ;
266  paraxe.add_double( g10[0], 6) ;
267  paraxe.add_double( g10[1], 7) ;
268  paraxe.add_double( g20[0], 8) ;
269  paraxe.add_double( g20[1], 9) ;
270  paraxe.add_double( g11[0], 10) ;
271  paraxe.add_double( g11[1], 11) ;
272  paraxe.add_double( g21[0], 12) ;
273  paraxe.add_double( g21[1], 13) ;
274  paraxe.add_double( g22[0], 14) ;
275  paraxe.add_double( g22[1], 15) ;
276  paraxe.add_double( bx[0], 16) ;
277  paraxe.add_double( bx[1], 17) ;
278  paraxe.add_double( by[0], 18) ;
279  paraxe.add_double( by[1], 19) ;
280  paraxe.add_double( bz[0], 20) ;
281  paraxe.add_double( bz[1], 21) ;
282  paraxe.add_double( dg00[0], 22) ;
283  paraxe.add_double( dg00[1], 23) ;
284  paraxe.add_double( dg10[0], 24) ;
285  paraxe.add_double( dg10[1], 25) ;
286  paraxe.add_double( dg20[0], 26) ;
287  paraxe.add_double( dg20[1], 27) ;
288  paraxe.add_double( dg11[0], 28) ;
289  paraxe.add_double( dg11[1], 29) ;
290  paraxe.add_double( dg21[0], 30) ;
291  paraxe.add_double( dg21[1], 31) ;
292  paraxe.add_double( dg22[0], 32) ;
293  paraxe.add_double( dg22[1], 33) ;
294  paraxe.add_double( dbx[0], 34) ;
295  paraxe.add_double( dbx[1], 35) ;
296  paraxe.add_double( dbz[0], 36) ;
297  paraxe.add_double( dbz[1], 37) ;
298  paraxe.add_double( dbymo[0], 38) ;
299  paraxe.add_double( dbymo[1], 39) ;
300  paraxe.add_double( d1sn2[0], 40) ;
301  paraxe.add_double( d1sn2[1], 41) ;
302  paraxe.add_double( unsn2[0], 42) ;
303  paraxe.add_double( unsn2[1], 43) ;
304  paraxe.add_double( omega, 44) ;
305 
306  int nitmax_axe = 200 ;
307  int nit_axe ;
308  double precis_axe = 1.e-13 ;
309 
310  x_axe = zerosec(fonc_binary_axe, paraxe, 0.9*ori_x1, 0.9*ori_x2,
311  precis_axe, nitmax_axe, nit_axe) ;
312 
313  cout << "Binary::orbit : Number of iterations in zerosec for x_axe : "
314  << nit_axe << endl ;
315  }
316 
317  cout << "Binary::orbit : x_axe [km] : " << x_axe / km << endl ;
318 
319 //-------------------------------------
320 // Calcul de la vitesse orbitale
321 //-------------------------------------
322 
323  Param parf ;
324  double ori_x0 = (et[0]->get_mp()).get_ori_x() ;
325  parf.add_double( ori_x0, 0) ;
326  parf.add_double( dnulg[0], 1) ;
327  parf.add_double( g00[0], 2) ;
328  parf.add_double( g10[0], 3) ;
329  parf.add_double( g20[0], 4) ;
330  parf.add_double( g11[0], 5) ;
331  parf.add_double( g21[0], 6) ;
332  parf.add_double( g22[0], 7) ;
333  parf.add_double( bx[0], 8) ;
334  parf.add_double( by[0], 9) ;
335  parf.add_double( bz[0], 10) ;
336  parf.add_double( dg00[0], 11) ;
337  parf.add_double( dg10[0], 12) ;
338  parf.add_double( dg20[0], 13) ;
339  parf.add_double( dg11[0], 14) ;
340  parf.add_double( dg21[0], 15) ;
341  parf.add_double( dg22[0], 16) ;
342  parf.add_double( dbx[0], 17) ;
343  parf.add_double( dbz[0], 18) ;
344  parf.add_double( dby[0], 19) ;
345  parf.add_double( d1sn2[0], 20) ;
346  parf.add_double( unsn2[0], 21) ;
347  parf.add_double( x_axe, 22) ;
348 
349 
350  double omega1 = fact_omeg_min * omega ;
351  double omega2 = fact_omeg_max * omega ;
352  cout << "Binary::orbit: omega1, omega2 [rad/s] : "
353  << omega1 * f_unit << " " << omega2 * f_unit << endl ;
354 
355 
356  // Search for the various zeros in the interval [omega1,omega2]
357  // ------------------------------------------------------------
358  int nsub = 50 ; // total number of subdivisions of the interval
359  Tbl* azer = 0x0 ;
360  Tbl* bzer = 0x0 ;
361  zero_list(fonc_binary_orbit, parf, omega1, omega2, nsub,
362  azer, bzer) ;
363 
364  // Search for the zero closest to the previous value of omega
365  // ----------------------------------------------------------
366  double omeg_min, omeg_max ;
367  int nzer = azer->get_taille() ; // number of zeros found by zero_list
368  cout << "Binary:orbit : " << nzer <<
369  " zero(s) found in the interval [omega1, omega2]." << endl ;
370  cout << "omega, omega1, omega2 : " << omega << " " << omega1
371  << " " << omega2 << endl ;
372  cout << "azer : " << *azer << endl ;
373  cout << "bzer : " << *bzer << endl ;
374 
375  if (nzer == 0) {
376  cout <<
377  "Binary::orbit: WARNING : no zero detected in the interval"
378  << endl << " [" << omega1 * f_unit << ", "
379  << omega2 * f_unit << "] rad/s !" << endl ;
380  omeg_min = omega1 ;
381  omeg_max = omega2 ;
382  }
383  else {
384  double dist_min = fabs(omega2 - omega1) ;
385  int i_dist_min = -1 ;
386  for (int i=0; i<nzer; i++) {
387  // Distance of previous value of omega from the center of the
388  // interval [azer(i), bzer(i)]
389  double dist = fabs( omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
390  if (dist < dist_min) {
391  dist_min = dist ;
392  i_dist_min = i ;
393  }
394  }
395  omeg_min = (*azer)(i_dist_min) ;
396  omeg_max = (*bzer)(i_dist_min) ;
397  }
398 
399  delete azer ; // Tbl allocated by zero_list
400  delete bzer ; //
401  cout << "Binary::orbit : interval selected for the search of the zero : "
402  << endl << " [" << omeg_min << ", " << omeg_max << "] = ["
403  << omeg_min * f_unit << ", " << omeg_max * f_unit << "] rad/s " << endl ;
404 
405  // Computation of the zero in the selected interval by the secant method
406  // ---------------------------------------------------------------------
407 
408  int nitermax = 200 ;
409  int niter ;
410  double precis = 1.e-13 ;
411  omega = zerosec_b(fonc_binary_orbit, parf, omeg_min, omeg_max,
412  precis, nitermax, niter) ;
413 
414  cout << "Binary::orbit : Number of iterations in zerosec for omega : "
415  << niter << endl ;
416 
417  cout << "Binary::orbit : omega [rad/s] : "
418  << omega * f_unit << endl ;
419 
420 
421 }
422 
423 
424 //-------------------------------------------------
425 // Function used for search of the rotation axis
426 //-------------------------------------------------
427 
428 double fonc_binary_axe(double x_rot, const Param& paraxe) {
429 
430  double ori_x1 = paraxe.get_double(0) ;
431  double ori_x2 = paraxe.get_double(1) ;
432  double dnulg_1 = paraxe.get_double(2) ;
433  double dnulg_2 = paraxe.get_double(3) ;
434  double g00_1 = paraxe.get_double(4) ;
435  double g00_2 = paraxe.get_double(5) ;
436  double g10_1 = paraxe.get_double(6) ;
437  double g10_2 = paraxe.get_double(7) ;
438  double g20_1 = paraxe.get_double(8) ;
439  double g20_2 = paraxe.get_double(9) ;
440  double g11_1 = paraxe.get_double(10) ;
441  double g11_2 = paraxe.get_double(11) ;
442  double g21_1 = paraxe.get_double(12) ;
443  double g21_2 = paraxe.get_double(13) ;
444  double g22_1 = paraxe.get_double(14) ;
445  double g22_2 = paraxe.get_double(15) ;
446  double bx_1 = paraxe.get_double(16) ;
447  double bx_2 = paraxe.get_double(17) ;
448  double by_1 = paraxe.get_double(18) ;
449  double by_2 = paraxe.get_double(19) ;
450  double bz_1 = paraxe.get_double(20) ;
451  double bz_2 = paraxe.get_double(21) ;
452  double dg00_1 = paraxe.get_double(22) ;
453  double dg00_2 = paraxe.get_double(23) ;
454  double dg10_1 = paraxe.get_double(24) ;
455  double dg10_2 = paraxe.get_double(25) ;
456  double dg20_1 = paraxe.get_double(26) ;
457  double dg20_2 = paraxe.get_double(27) ;
458  double dg11_1 = paraxe.get_double(28) ;
459  double dg11_2 = paraxe.get_double(29) ;
460  double dg21_1 = paraxe.get_double(30) ;
461  double dg21_2 = paraxe.get_double(31) ;
462  double dg22_1 = paraxe.get_double(32) ;
463  double dg22_2 = paraxe.get_double(33) ;
464  double dbx_1 = paraxe.get_double(34) ;
465  double dbx_2 = paraxe.get_double(35) ;
466  double dbz_1 = paraxe.get_double(36) ;
467  double dbz_2 = paraxe.get_double(37) ;
468  double dbymo_1 = paraxe.get_double(38) ;
469  double dbymo_2 = paraxe.get_double(39) ;
470  double d1sn2_1 = paraxe.get_double(40) ;
471  double d1sn2_2 = paraxe.get_double(41) ;
472  double unsn2_1 = paraxe.get_double(42) ;
473  double unsn2_2 = paraxe.get_double(43) ;
474  double omega = paraxe.get_double(44) ;
475 
476  double om2_star1 ;
477  double om2_star2 ;
478 
479  double x1 = ori_x1 - x_rot ;
480  double x2 = ori_x2 - x_rot ;
481 
482  double bymxo_1 = by_1-x1*omega ;
483  double bymxo_2 = by_2-x2*omega ;
484 
485 
486  double beta1 = g00_1*bx_1*bx_1 + 2*g10_1*bx_1*bymxo_1 + 2*g20_1*bx_1*bz_1 ;
487  double beta2 = g11_1*bymxo_1*bymxo_1 + 2*g21_1*bz_1*bymxo_1
488  + g22_1*bz_1*bz_1 ;
489 
490  double beta_1 = beta1 + beta2 ;
491 
492 
493  double delta1 = dg00_1*bx_1*bx_1 + 2*g00_1*dbx_1*bx_1
494  + 2*dg10_1*bx_1*bymxo_1 ;
495  double delta2 = 2*g10_1*bymxo_1*dbx_1 + 2*g10_1*bx_1*dbymo_1
496  + 2*dg20_1*bx_1*bz_1 ;
497  double delta3 = 2*g20_1*bx_1*dbz_1 +2*g20_1*bz_1*dbx_1 + dg11_1*bymxo_1*bymxo_1 ;
498  double delta4 = 2*g11_1*bymxo_1*dbymo_1 + 2*dg21_1*bz_1*bymxo_1;
499  double delta5 = 2*g21_1*bymxo_1*dbz_1 +2*g21_1*bz_1*dbymo_1
500  + dg22_1*bz_1*bz_1 + 2*g22_1*bz_1*dbz_1 ;
501 
502  double delta_1 = delta1 + delta2 + delta3 + delta4 + delta5 ;
503 
504  // Computation of omega for star 1
505  //---------------------------------
506 
507  om2_star1 = dnulg_1 / (beta_1/(omega*omega)*(dnulg_1*unsn2_1 + d1sn2_1/2.)
508  + unsn2_1*delta_1/(omega*omega)/2.) ;
509 
510 
511 
512  double beta3 = g00_2*bx_2*bx_2 + 2*g10_2*bx_2*bymxo_2 + 2*g20_2*bx_2*bz_2 ;
513  double beta4 = g11_2*bymxo_2*bymxo_2 + 2*g21_2*bz_2*bymxo_2
514  + g22_2*bz_2*bz_2 ;
515 
516  double beta_2 = beta3 + beta4 ;
517 
518 
519  double delta6 = dg00_2*bx_2*bx_2 + 2*g00_2*dbx_2*bx_2
520  + 2*dg10_2*bx_2*bymxo_2 ;
521  double delta7 = 2*g10_2*bymxo_2*dbx_2 + 2*g10_2*bx_2*dbymo_2
522  + 2*dg20_2*bx_2*bz_2 ;
523  double delta8 = 2*g20_2*bx_2*dbz_2 +2*g20_2*bz_2*dbx_2
524  + dg11_2*bymxo_2*bymxo_2 ;
525  double delta9 = 2*g11_2*bymxo_2*dbymo_2 + 2*dg21_2*bz_2*bymxo_2;
526  double delta10 = 2*g21_2*bymxo_2*dbz_2 +2*g21_2*bz_2*dbymo_2
527  + dg22_2*bz_2*bz_2 + 2*g22_2*bz_2*dbz_2 ;
528 
529  double delta_2 = delta6 + delta7 + delta8 + delta9 + delta10 ;
530 
531  // Computation of omega for star 2
532  //---------------------------------
533 
534  om2_star2 = dnulg_2 / (beta_2/(omega*omega)*(dnulg_2*unsn2_2 + d1sn2_2/2.)
535  + unsn2_2*delta_2/(omega*omega)/2.) ;
536  ;
537 
538  return om2_star1 - om2_star2 ;
539 
540 }
541 
542 //----------------------------------------------------------------------------
543 // Fonction utilisee pour la recherche de omega par la methode de la secante
544 //----------------------------------------------------------------------------
545 
546 double fonc_binary_orbit(double om, const Param& parf) {
547 
548  double xc = parf.get_double(0) ;
549  double dnulg = parf.get_double(1) ;
550  double g00 = parf.get_double(2) ;
551  double g10 = parf.get_double(3) ;
552  double g20 = parf.get_double(4) ;
553  double g11 = parf.get_double(5) ;
554  double g21 = parf.get_double(6) ;
555  double g22 = parf.get_double(7) ;
556  double bx = parf.get_double(8) ;
557  double by = parf.get_double(9) ;
558  double bz = parf.get_double(10) ;
559  double dg00 = parf.get_double(11) ;
560  double dg10 = parf.get_double(12) ;
561  double dg20 = parf.get_double(13) ;
562  double dg11 = parf.get_double(14) ;
563  double dg21 = parf.get_double(15) ;
564  double dg22 = parf.get_double(16) ;
565  double dbx = parf.get_double(17) ;
566  double dbz = parf.get_double(18) ;
567  double dby = parf.get_double(19) ;
568  double d1sn2 = parf.get_double(20) ;
569  double unsn2 = parf.get_double(21) ;
570  double x_axe = parf.get_double(22) ;
571 
572 
573  double dbymo = dby - om ;
574  double xx = xc - x_axe ;
575 
576  double bymxo = by-xx*om ;
577 
578  // bymxo = - bymxo ;
579  //dbymo = - dbymo ;
580 
581  double beta1 = g00*bx*bx + 2*g10*bx*bymxo + 2*g20*bx*bz ;
582  double beta2 = g11*bymxo*bymxo + 2*g21*bz*bymxo+ g22*bz*bz ;
583  double beta = beta1 + beta2 ;
584 
585  double alpha = 1 - unsn2*beta ;
586 
587  double delta1 = dg00*bx*bx + 2*g00*dbx*bx + 2*dg10*bx*bymxo ;
588  double delta2 = 2*g10*bymxo*dbx + 2*g10*bx*dbymo + 2*dg20*bx*bz ;
589  double delta3 = 2*g20*bx*dbz +2*g20*bz*dbx + dg11*bymxo*bymxo ;
590  double delta4 = 2*g11*bymxo*dbymo + 2*dg21*bz*bymxo;
591  double delta5 = 2*g21*bymxo*dbz +2*g21*bz*dbymo + dg22*bz*bz
592  + 2*g22*bz*dbz ;
593 
594  double delta = delta1 + delta2 + delta3 + delta4 + delta5 ;
595 
596  // Difference entre les 2 termes de l'eq.(95) de Gourgoulhon et.al (2001)
597  //centre de l'etoile
598  //-----------------------------------------------------------------------
599 
600  double diff = dnulg + (1/(2.*alpha))*(-d1sn2*beta - unsn2*delta) ;
601 
602  /*
603  double bpb = om*om*xx*xx - 2*om*by*xx + by*by ;
604  double dphi_cent = (g00*unsn2*( om*(by+xx*dby) - om*om*xx - by*dby )
605  - 0.5*bpb*(dg00 + g00*d1sn2/unsn2)*unsn2 )
606  / (1 - g00*unsn2*bpb) ;
607  double diff = dnulg + dphi_cent ;
608 
609 
610  cout.precision(8) ;
611  cout << "bpb = " << bpb << endl ;
612  cout << "om = " << om << endl ;
613  cout << "by = " << by << endl ;
614  cout << "xx = " << xx << endl ;
615  cout << "dby = " << dby << endl ;
616  cout << "part11 = " << g00*unsn2 << endl ;
617  cout << "part12 = " << om*(by+xx*dby) << endl ;
618  cout << "part13 = " <<- om*om*xx << endl ;
619  cout << "part14 = " << - by*dby << endl ;
620  cout << "part1 = " << g00*unsn2*( om*(by+xx*dby) - om*om*xx - by*dby ) << endl ;
621  cout << "part2 = " << - 0.5*bpb*(dg00 + g00*d1sn2/unsn2)*unsn2 << endl ;
622  cout << "part3 = " << (1 - g00*unsn2*bpb) << endl ;
623  cout << "dnulg = " << dnulg << endl ;
624  cout << "dphi_cent = " << dphi_cent << endl ;
625  cout << "diff = " << diff << endl ;
626  cout << endl ;
627  */
628 
629  return diff ;
630 
631 
632 
633 }
634 
635 
636 
637 }
Metric for tensor calculation.
Definition: metric.h:90
const Vector & get_beta() const
Returns the shift vector .
Definition: star.h:402
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition: binary.h:94
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
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 & get_logn_comp() const
Returns the part of the lapse logarithm (gravitational potential at the Newtonian limit) generated pr...
Definition: star.h:811
const Scalar & dsdx() const
Returns of *this , where .
Definition: scalar_deriv.C:266
const Metric & get_gamma() const
Returns the 3-metric .
Definition: star.h:409
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.
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 Scalar & get_loggam() const
Returns the logarithm of the Lorentz factor between the fluid and the co-orbiting observer...
Definition: star.h:792
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition: cmp_math.C:484
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
Parameter storage.
Definition: param.h:125
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...
Definition: binary_orbite.C:90
double x_axe
Absolute X coordinate of the rotation axis.
Definition: binary.h:98
Tensor handling.
Definition: tensor.h:294
Star_bin * 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.h:89
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:283
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:417
const Scalar & get_ent() const
Returns the enthalpy field.
Definition: star.h:364
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 void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
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_nn() const
Returns the lapse function N.
Definition: star.h:399
const Scalar & get_logn_auto() const
Returns the part of the lapse logarithm (gravitational potential at the Newtonian limit) generated pr...
Definition: star.h:806
virtual double xa_barycenter() const
Absolute coordinate X of the barycenter of the baryon density,.