61#include "bin_bhns_extr.h"
64#include "utilitaires.h"
68double fonc_bhns_orbit_ks(
double,
const Param& ) ;
69double fonc_bhns_orbit_cf(
double,
const Param& ) ;
74 double fact_omeg_max) {
82 double dnulg, asn2, dasn2, nx, dnx, ny, dny, npn, dnpn ;
87 const Cmp& logn_auto_regu =
star.get_logn_auto_regu()() ;
88 const Cmp& loggam =
star.get_loggam()() ;
89 const Cmp& nnn =
star.get_nnn()() ;
90 const Cmp& a_car =
star.get_a_car()() ;
92 const Tenseur& d_logn_auto_div =
star.get_d_logn_auto_div() ;
94 Tenseur dln_auto_div = d_logn_auto_div ;
106 const Tenseur& d_logn_comp =
star.get_d_logn_comp() ;
108 Tenseur dln_comp = d_logn_comp ;
133 if (fabs(mp.get_rot_phi()) < 1.e-14) {
137 if (fabs(mp.get_rot_phi() - M_PI) < 1.e-14) {
141 cout <<
"Bin_bhns_extr::orbit_omega : unknown value of rot_phi !"
147 Cmp tmp = logn_auto_regu + loggam ;
150 dnulg = dln_auto_div(0)(0, 0, 0, 0) + dln_comp(0)(0, 0, 0, 0)
151 + factx * tmp.
dsdx()(0, 0, 0, 0) ;
156 double nc = nnn(0, 0, 0, 0) ;
157 double a2c = a_car(0, 0, 0, 0) ;
158 asn2 = a2c / (nc * nc) ;
160 if (
star.is_relativistic() ) {
166 double da2c = factx * a_car.
dsdx()(0, 0, 0, 0) ;
167 double dnc = factx * nnn.
dsdx()(0, 0, 0, 0) ;
169 dasn2 = ( da2c - 2. * a2c / nc * dnc ) / (nc*nc) ;
174 nx = shift(0)(0, 0, 0, 0) ;
179 dnx = factx * shift(0).dsdx()(0, 0, 0, 0) ;
184 ny = shift(1)(0, 0, 0, 0) ;
189 dny = factx * shift(1).dsdx()(0, 0, 0, 0) ;
196 npn = tmp(0, 0, 0, 0) ;
202 dnpn = factx * tmp.
dsdx()(0, 0, 0, 0) ;
206 cout <<
"Bin_bhns_extr::orbit_omega : "
207 <<
"It should be the relativistic calculation !" << endl ;
211 cout <<
"Bin_bhns_extr::orbit_omega: coord. ori. d(nu+log(Gam))/dX : "
213 cout <<
"Bin_bhns_extr::orbit_omega: coord. ori. A^2/N^2 : "
215 cout <<
"Bin_bhns_extr::orbit_omega: coord. ori. d(A^2/N^2)/dX : "
217 cout <<
"Bin_bhns_extr::orbit_omega: coord. ori. N^X : " << nx << endl ;
218 cout <<
"Bin_bhns_extr::orbit_omega: coord. ori. dN^X/dX : "
220 cout <<
"Bin_bhns_extr::orbit_omega: coord. ori. N^Y : " << ny << endl ;
221 cout <<
"Bin_bhns_extr::orbit_omega: coord. ori. dN^Y/dX : "
223 cout <<
"Bin_bhns_extr::orbit_omega: coord. ori. N.N : " << npn << endl ;
224 cout <<
"Bin_bhns_extr::orbit_omega: coord. ori. d(N.N)/dX : "
230 int relat = (
star.is_relativistic() ) ? 1 : 0 ;
232 double ori_x = (
star.get_mp()).get_ori_x() ;
247 double omega1 = fact_omeg_min *
omega ;
248 double omega2 = fact_omeg_max *
omega ;
250 cout <<
"Bin_bhns_extr::orbit_omega: omega1, omega2 [rad/s] : "
251 << omega1 * f_unit <<
" " << omega2 * f_unit << endl ;
258 zero_list(fonc_bhns_orbit_ks, parf, omega1, omega2, nsub,
263 double omeg_min, omeg_max ;
265 cout <<
"Bin_bhns_extr::orbit_omega : " << nzer <<
266 "zero(s) found in the interval [omega1, omega2]." << endl ;
267 cout <<
"omega, omega1, omega2 : " <<
omega <<
" " << omega1
268 <<
" " << omega2 << endl ;
269 cout <<
"azer : " << *azer << endl ;
270 cout <<
"bzer : " << *bzer << endl ;
273 cout <<
"Bin_bhns_extr::orbit_omega: WARNING : "
274 <<
"no zero detected in the interval" << endl
275 <<
" [" << omega1 * f_unit <<
", "
276 << omega2 * f_unit <<
"] rad/s !" << endl ;
281 double dist_min = fabs(omega2 - omega1) ;
282 int i_dist_min = -1 ;
283 for (
int i=0; i<nzer; i++) {
286 double dist = fabs(
omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
288 if (dist < dist_min) {
293 omeg_min = (*azer)(i_dist_min) ;
294 omeg_max = (*bzer)(i_dist_min) ;
300 cout <<
"Bin_bhns_extr:orbit_omega : "
301 <<
"interval selected for the search of the zero : "
302 << endl <<
" [" << omeg_min <<
", " << omeg_max <<
"] = ["
303 << omeg_min * f_unit <<
", " << omeg_max * f_unit <<
"] rad/s "
311 double precis = 1.e-13 ;
313 precis, nitermax, niter) ;
315 cout <<
"Bin_bhns_extr::orbit_omega : "
316 <<
"Number of iterations in zerosec for omega : "
319 cout <<
"Bin_bhns_extr::orbit_omega : omega [rad/s] : "
320 <<
omega * f_unit << endl ;
326 double fact_omeg_max) {
334 double dnulg, asn2, dasn2, ny, dny, npn, dnpn ;
336 const Map& mp =
star.get_mp() ;
338 const Cmp& logn_auto_regu =
star.get_logn_auto_regu()() ;
339 const Cmp& logn_comp =
star.get_logn_comp()() ;
340 const Cmp& loggam =
star.get_loggam()() ;
341 const Cmp& nnn =
star.get_nnn()() ;
342 const Cmp& a_car =
star.get_a_car()() ;
344 const Tenseur& d_logn_auto_div =
star.get_d_logn_auto_div() ;
346 Tenseur dln_auto_div = d_logn_auto_div ;
365 if (fabs(mp.get_rot_phi()) < 1.e-14) {
369 if (fabs(mp.get_rot_phi() - M_PI) < 1.e-14) {
374 "Bin_bhns_extr::orbit_omega_cfc : unknown value of rot_phi !"
380 Cmp tmp = logn_auto_regu + logn_comp + loggam ;
383 dnulg = dln_auto_div(0)(0, 0, 0, 0)
384 + factx * tmp.
dsdx()(0, 0, 0, 0) ;
389 double nc = nnn(0, 0, 0, 0) ;
390 double a2c = a_car(0, 0, 0, 0) ;
391 asn2 = a2c / (nc * nc) ;
393 if (
star.is_relativistic() ) {
399 double da2c = factx * a_car.
dsdx()(0, 0, 0, 0) ;
400 double dnc = factx * nnn.
dsdx()(0, 0, 0, 0) ;
402 dasn2 = ( da2c - 2. * a2c / nc * dnc ) / (nc*nc) ;
407 ny = shift(1)(0, 0, 0, 0) ;
412 dny = factx * shift(1).dsdx()(0, 0, 0, 0) ;
419 npn = tmp(0, 0, 0, 0) ;
425 dnpn = factx * tmp.
dsdx()(0, 0, 0, 0) ;
429 cout <<
"Bin_bhns_extr::orbit_omega_cfc : "
430 <<
"It should be the relativistic calculation !" << endl ;
434 cout <<
"Bin_bhns_extr::orbit_omega_cfc: coord. ori. d(nu+log(Gam))/dX : "
436 cout <<
"Bin_bhns_extr::orbit_omega_cfc: coord. ori. A^2/N^2 : "
438 cout <<
"Bin_bhns_extr::orbit_omega_cfc: coord. ori. d(A^2/N^2)/dX : "
440 cout <<
"Bin_bhns_extr::orbit_omega_cfc: coord. ori. N^Y : "
442 cout <<
"Bin_bhns_extr::orbit_omega_cfc: coord. ori. dN^Y/dX : "
444 cout <<
"Bin_bhns_extr::orbit_omega_cfc: coord. ori. N.N : "
446 cout <<
"Bin_bhns_extr::orbit_omega_cfc: coord. ori. d(N.N)/dX : "
452 int relat = (
star.is_relativistic() ) ? 1 : 0 ;
453 double ori_x = (
star.get_mp()).get_ori_x() ;
466 double omega1 = fact_omeg_min *
omega ;
467 double omega2 = fact_omeg_max *
omega ;
469 cout <<
"Bin_bhns_extr::orbit_omega_cfc: omega1, omega2 [rad/s] : "
470 << omega1 * f_unit <<
" " << omega2 * f_unit << endl ;
477 zero_list(fonc_bhns_orbit_cf, parf, omega1, omega2, nsub,
482 double omeg_min, omeg_max ;
484 cout <<
"Bin_bhns_extr::orbit_omega_cfc : " << nzer <<
485 "zero(s) found in the interval [omega1, omega2]." << endl ;
486 cout <<
"omega, omega1, omega2 : " <<
omega <<
" " << omega1
487 <<
" " << omega2 << endl ;
488 cout <<
"azer : " << *azer << endl ;
489 cout <<
"bzer : " << *bzer << endl ;
492 cout <<
"Bin_bhns_extr::orbit_omega_cfc: WARNING : "
493 <<
"no zero detected in the interval" << endl
494 <<
" [" << omega1 * f_unit <<
", "
495 << omega2 * f_unit <<
"] rad/s !" << endl ;
500 double dist_min = fabs(omega2 - omega1) ;
501 int i_dist_min = -1 ;
502 for (
int i=0; i<nzer; i++) {
505 double dist = fabs(
omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
507 if (dist < dist_min) {
512 omeg_min = (*azer)(i_dist_min) ;
513 omeg_max = (*bzer)(i_dist_min) ;
519 cout <<
"Bin_bhns_extr:orbit_omega_cfc : "
520 <<
"interval selected for the search of the zero : "
521 << endl <<
" [" << omeg_min <<
", " << omeg_max <<
"] = ["
522 << omeg_min * f_unit <<
", " << omeg_max * f_unit <<
"] rad/s "
530 double precis = 1.e-13 ;
532 precis, nitermax, niter) ;
534 cout <<
"Bin_bhns_extr::orbit_omega_cfc : "
535 <<
"Number of iterations in zerosec for omega : "
538 cout <<
"Bin_bhns_extr::orbit_omega_cfc : omega [rad/s] : "
539 <<
omega * f_unit << endl ;
548double fonc_bhns_orbit_ks(
double om,
const Param& parf) {
569 double bpb = om2 * xc*xc - 2.*om * ny * xc + npn + 2.*msr * nx*nx;
571 dphi_cent = ( asn2* ( om* (ny + xc*dny) - om2*xc - 0.5*dnpn
572 - 2.*msr * nx * dnx + msr * nx*nx / xc )
573 - 0.5*bpb* dasn2 ) / ( 1 - asn2 * bpb ) ;
576 cout <<
"Bin_bhns_extr::orbit_omega_ks : "
577 <<
"It should be the relativistic calculation !" << endl ;
581 return dnulg + dphi_cent ;
585double fonc_bhns_orbit_cf(
double om,
const Param& parf) {
587 int relat = parf.get_int() ;
589 double xc = parf.get_double(0) ;
590 double dnulg = parf.get_double(1) ;
591 double asn2 = parf.get_double(2) ;
592 double dasn2 = parf.get_double(3) ;
593 double ny = parf.get_double(4) ;
594 double dny = parf.get_double(5) ;
595 double npn = parf.get_double(6) ;
596 double dnpn = parf.get_double(7) ;
603 double bpb = om2 * xc*xc - 2.*om * ny * xc + npn ;
605 dphi_cent = ( asn2* ( om* (ny + xc*dny) - om2*xc - 0.5*dnpn )
606 - 0.5*bpb* dasn2 ) / ( 1 - asn2 * bpb ) ;
609 cout <<
"Bin_bhns_extr::orbit_omega_cf : "
610 <<
"It should be the relativistic calculation !" << endl ;
614 return dnulg + dphi_cent ;
void orbit_omega_ks(double fact_omeg_min, double fact_omeg_max)
Computes the orbital angular velocity {\tt omega} in the Kerr-Schild background metric.
const Base_vect_cart ref_triad
Cartesian triad of the absolute reference frame.
Et_bin_bhns_extr star
Neutron star.
double separ
Absolute orbital separation between two centers of BH and NS.
void orbit_omega_cf(double fact_omeg_min, double fact_omeg_max)
Computes the orbital angular velocity {\tt omega} in the conformally flat background metric.
double mass_bh
Gravitational mass of BH.
double omega
Angular velocity with respect to an asymptotically inertial observer.
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
const Cmp & dsdx() const
Returns of *this , where .
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
int get_taille() const
Gives the total size (ie dim.taille).
Tensor handling *** DEPRECATED : use class Tensor instead ***.
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
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.
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
Map(const Mg3d &)
Constructor from a multi-domain 3D grid.
Standard units of space, time and mass.