81 #include "utilitaires.h" 85 double fonc_binary_axe(
double ,
const Param& ) ;
86 double fonc_binary_orbit(
double ,
const Param& ) ;
90 void Binary::orbit(
double fact_omeg_min,
double fact_omeg_max,
double& xgg1,
100 double g00[2], g10[2], g20[2], g11[2], g21[2], g22[2], bx[2], by[2] ;
102 double bz[2], d1sn2[2], unsn2[2] ;
104 double dnulg[2], xgg[2], ori_x[2], dg00[2], dg10[2], dg20[2], dg11[2] ;
106 double dg21[2], dg22[2], dbx[2], dby[2], dbz[2], dbymo[2] ;
108 for (
int i=0; i<2; i++) {
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) ;
134 const Scalar& bbx = shift(1) ;
135 const Scalar& bby = shift(2) ;
136 const Scalar& bbz = shift(3) ;
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 ;
145 cout <<
"bbx" << endl <<
norme(bbx) << endl ;
146 cout <<
"bby" << endl <<
norme(bby) << endl ;
147 cout <<
"bbz" << endl <<
norme(bbz) << endl ;
153 Scalar tmp = logn_auto + logn_comp + loggam ;
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 ;
163 cout <<
"dnulg = " << dnulg[i] << endl ;
201 d1sn2[i] = (1/(nn*nn)).dsdx().val_grid_point(0,0,0,0) ;
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 ;
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 ;
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 ;
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 ;
234 ori_x[i] = (
et[i]->
get_mp()).get_ori_x() ;
247 double ori_x1 = ori_x[0] ;
248 double ori_x2 = ori_x[1] ;
306 int nitmax_axe = 200 ;
308 double precis_axe = 1.e-13 ;
310 x_axe =
zerosec(fonc_binary_axe, paraxe, 0.9*ori_x1, 0.9*ori_x2,
311 precis_axe, nitmax_axe, nit_axe) ;
313 cout <<
"Binary::orbit : Number of iterations in zerosec for x_axe : " 317 cout <<
"Binary::orbit : x_axe [km] : " <<
x_axe / km << endl ;
324 double ori_x0 = (
et[0]->
get_mp()).get_ori_x() ;
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 ;
361 zero_list(fonc_binary_orbit, parf, omega1, omega2, nsub,
366 double omeg_min, omeg_max ;
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 ;
377 "Binary::orbit: WARNING : no zero detected in the interval" 378 << endl <<
" [" << omega1 * f_unit <<
", " 379 << omega2 * f_unit <<
"] rad/s !" << endl ;
384 double dist_min = fabs(omega2 - omega1) ;
385 int i_dist_min = -1 ;
386 for (
int i=0; i<nzer; i++) {
389 double dist = fabs(
omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
390 if (dist < dist_min) {
395 omeg_min = (*azer)(i_dist_min) ;
396 omeg_max = (*bzer)(i_dist_min) ;
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 ;
410 double precis = 1.e-13 ;
412 precis, nitermax, niter) ;
414 cout <<
"Binary::orbit : Number of iterations in zerosec for omega : " 417 cout <<
"Binary::orbit : omega [rad/s] : " 418 <<
omega * f_unit << endl ;
428 double fonc_binary_axe(
double x_rot,
const Param& paraxe) {
479 double x1 = ori_x1 - x_rot ;
480 double x2 = ori_x2 - x_rot ;
482 double bymxo_1 = by_1-x1*
omega ;
483 double bymxo_2 = by_2-x2*
omega ;
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
490 double beta_1 = beta1 + beta2 ;
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 ;
502 double delta_1 = delta1 + delta2 + delta3 + delta4 + delta5 ;
507 om2_star1 = dnulg_1 / (beta_1/(
omega*
omega)*(dnulg_1*unsn2_1 + d1sn2_1/2.)
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
516 double beta_2 = beta3 + beta4 ;
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 ;
529 double delta_2 = delta6 + delta7 + delta8 + delta9 + delta10 ;
534 om2_star2 = dnulg_2 / (beta_2/(
omega*
omega)*(dnulg_2*unsn2_2 + d1sn2_2/2.)
538 return om2_star1 - om2_star2 ;
546 double fonc_binary_orbit(
double om,
const Param& parf) {
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) ;
573 double dbymo = dby - om ;
574 double xx = xc -
x_axe ;
576 double bymxo = by-xx*om ;
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 ;
585 double alpha = 1 - unsn2*beta ;
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
594 double delta = delta1 + delta2 + delta3 + delta4 + delta5 ;
600 double diff = dnulg + (1/(2.*alpha))*(-d1sn2*beta - unsn2*delta) ;
Metric for tensor calculation.
const Vector & get_beta() const
Returns the shift vector .
Standard units of space, time and mass.
double omega
Angular velocity with respect to an asymptotically inertial observer.
Tensor field of valence 0 (or component of a tensorial field).
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.
const Scalar & get_logn_comp() const
Returns the part of the lapse logarithm (gravitational potential at the Newtonian limit) generated pr...
const Scalar & dsdx() const
Returns of *this , where .
const Metric & get_gamma() const
Returns the 3-metric .
const Map & get_mp() const
Returns the mapping.
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.
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.
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 & get_loggam() const
Returns the logarithm of the Lorentz factor between the fluid and the co-orbiting observer...
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.
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.
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]...
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
int get_taille() const
Gives the total size (ie dim.taille)
const Scalar & get_ent() const
Returns the enthalpy field.
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
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.
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
const Scalar & get_nn() const
Returns the lapse function N.
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,.