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
31char 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
81namespace Lorene {
82double fonc_binary_axe(double , const Param& ) ;
83double fonc_binary_orbit(double , const Param& ) ;
84
85//******************************************************************************
86
87void Binary::orbit(double fact_omeg_min, double fact_omeg_max, double& xgg1,
88 double& xgg2) {
89
90using 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
425double 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
543double 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}
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
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition binary.h:94
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.
double x_axe
Absolute X coordinate of the rotation axis.
Definition binary.h:98
Metric for tensor calculation.
Definition metric.h:90
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:280
Parameter storage.
Definition param.h:125
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition param.C:315
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition param.C:361
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
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 & dsdx() const
Returns of *this , where .
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,.
const Scalar & get_loggam() const
Returns the logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition star.h:792
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 & get_nn() const
Returns the lapse function N.
Definition star.h:399
const Map & get_mp() const
Returns the mapping.
Definition star.h:355
const Eos & get_eos() const
Returns the equation of state.
Definition star.h:361
const Metric & get_gamma() const
Returns the 3-metric .
Definition star.h:409
const Scalar & get_ent() const
Returns the enthalpy field.
Definition star.h:364
const Vector & get_beta() const
Returns the shift vector .
Definition star.h:402
Basic array class.
Definition tbl.h:161
int get_taille() const
Gives the total size (ie dim.taille)
Definition tbl.h:397
Tensor handling.
Definition tensor.h:288
Tensor field of valence 1.
Definition vector.h:188
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
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
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.
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
virtual void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.