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