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 $" ;
78#include "utilitaires.h"
82double fonc_binary_axe(
double ,
const Param& ) ;
83double fonc_binary_orbit(
double ,
const Param& ) ;
87void Binary::orbit(
double fact_omeg_min,
double fact_omeg_max,
double& xgg1,
97 double g00[2], g10[2], g20[2], g11[2], g21[2], g22[2], bx[2], by[2] ;
99 double bz[2], d1sn2[2], unsn2[2] ;
101 double dnulg[2], xgg[2], ori_x[2], dg00[2], dg10[2], dg20[2], dg11[2] ;
103 double dg21[2], dg22[2], dbx[2], dby[2], dbz[2], dbymo[2] ;
105 for (
int i=0; i<2; i++) {
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) ;
131 const Scalar& bbx = shift(1) ;
132 const Scalar& bby = shift(2) ;
133 const Scalar& bbz = shift(3) ;
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 ;
142 cout <<
"bbx" << endl <<
norme(bbx) << endl ;
143 cout <<
"bby" << endl <<
norme(bby) << endl ;
144 cout <<
"bbz" << endl <<
norme(bbz) << endl ;
150 Scalar tmp = logn_auto + logn_comp + loggam ;
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 ;
160 cout <<
"dnulg = " << dnulg[i] << endl ;
198 d1sn2[i] = (1/(nn*nn)).dsdx().val_grid_point(0,0,0,0) ;
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 ;
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 ;
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 ;
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 ;
231 ori_x[i] = (
et[i]->
get_mp()).get_ori_x() ;
244 double ori_x1 = ori_x[0] ;
245 double ori_x2 = ori_x[1] ;
303 int nitmax_axe = 200 ;
305 double precis_axe = 1.e-13 ;
307 x_axe =
zerosec(fonc_binary_axe, paraxe, 0.9*ori_x1, 0.9*ori_x2,
308 precis_axe, nitmax_axe, nit_axe) ;
310 cout <<
"Binary::orbit : Number of iterations in zerosec for x_axe : "
314 cout <<
"Binary::orbit : x_axe [km] : " <<
x_axe / km << endl ;
321 double ori_x0 = (
et[0]->
get_mp()).get_ori_x() ;
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 ;
358 zero_list(fonc_binary_orbit, parf, omega1, omega2, nsub,
363 double omeg_min, omeg_max ;
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 ;
374 "Binary::orbit: WARNING : no zero detected in the interval"
375 << endl <<
" [" << omega1 * f_unit <<
", "
376 << omega2 * f_unit <<
"] rad/s !" << endl ;
381 double dist_min = fabs(omega2 - omega1) ;
382 int i_dist_min = -1 ;
383 for (
int i=0; i<nzer; i++) {
386 double dist = fabs(
omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
387 if (dist < dist_min) {
392 omeg_min = (*azer)(i_dist_min) ;
393 omeg_max = (*bzer)(i_dist_min) ;
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 ;
407 double precis = 1.e-13 ;
409 precis, nitermax, niter) ;
411 cout <<
"Binary::orbit : Number of iterations in zerosec for omega : "
414 cout <<
"Binary::orbit : omega [rad/s] : "
415 <<
omega * f_unit << endl ;
425double fonc_binary_axe(
double x_rot,
const Param& paraxe) {
476 double x1 = ori_x1 - x_rot ;
477 double x2 = ori_x2 - x_rot ;
479 double bymxo_1 = by_1-x1*omega ;
480 double bymxo_2 = by_2-x2*omega ;
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
487 double beta_1 = beta1 + beta2 ;
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 ;
499 double delta_1 = delta1 + delta2 + delta3 + delta4 + delta5 ;
504 om2_star1 = dnulg_1 / (beta_1/(omega*omega)*(dnulg_1*unsn2_1 + d1sn2_1/2.)
505 + unsn2_1*delta_1/(omega*omega)/2.) ;
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
513 double beta_2 = beta3 + beta4 ;
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 ;
526 double delta_2 = delta6 + delta7 + delta8 + delta9 + delta10 ;
531 om2_star2 = dnulg_2 / (beta_2/(omega*omega)*(dnulg_2*unsn2_2 + d1sn2_2/2.)
532 + unsn2_2*delta_2/(omega*omega)/2.) ;
535 return om2_star1 - om2_star2 ;
543double fonc_binary_orbit(
double om,
const Param& parf) {
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) ;
570 double dbymo = dby - om ;
571 double xx = xc - x_axe ;
573 double bymxo = by-xx*om ;
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 ;
582 double alpha = 1 - unsn2*beta ;
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
591 double delta = delta1 + delta2 + delta3 + delta4 + delta5 ;
597 double diff = dnulg + (1/(2.*alpha))*(-d1sn2*beta - unsn2*delta) ;
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]...
double omega
Angular velocity with respect to an asymptotically inertial observer.
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.
Metric for tensor calculation.
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Tensor field of valence 0 (or component of a tensorial field).
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
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...
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.
const Scalar & get_logn_comp() const
Returns the part of the lapse logarithm (gravitational potential at the Newtonian limit) generated pr...
const Scalar & get_nn() const
Returns the lapse function N.
const Map & get_mp() const
Returns the mapping.
const Eos & get_eos() const
Returns the equation of state.
const Metric & get_gamma() const
Returns the 3-metric .
const Scalar & get_ent() const
Returns the enthalpy field.
const Vector & get_beta() const
Returns the shift vector .
int get_taille() const
Gives the total size (ie dim.taille)
Tensor field of valence 1.
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.
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.
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.
virtual void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Standard units of space, time and mass.