67 #include "bin_ns_bh.h" 70 #include "utilitaires.h" 74 double fonc_bin_ns_bh_orbit(
double ,
const Param& ) ;
86 double dnulg, asn2, dasn2, ny, dny, npn, dnpn ;
95 Cmp confpsi_q =
pow(confpsi, 4.) ;
113 cout <<
"Bin_ns_bh::orbit_omega : unknown value of rot_phi !" 119 Cmp tmp =
log( nnn ) + loggam ;
123 dnulg = factx * tmp.
dsdx()(0, 0, 0, 0) ;
128 double nc = nnn(0, 0, 0, 0) ;
129 double a2c = confpsi_q(0, 0, 0, 0) ;
130 asn2 = a2c / (nc * nc) ;
137 double da2c = factx * confpsi_q.
dsdx()(0, 0, 0, 0) ;
138 double dnc = factx * nnn.
dsdx()(0, 0, 0, 0) ;
140 dasn2 = ( da2c - 2 * a2c / nc * dnc ) / (nc*nc) ;
145 ny = shift(1)(0, 0, 0, 0) ;
150 dny = factx * shift(1).dsdx()(0, 0, 0, 0) ;
157 npn = tmp(0, 0, 0, 0) ;
163 dnpn = factx * tmp.
dsdx()(0, 0, 0, 0) ;
167 cout <<
"Bin_ns_bh::orbit_omega : " 168 <<
"It should be the relativistic calculation !" << endl ;
172 cout <<
"Bin_ns_bh::orbit_omega: central d(nu+log(Gam))/dX : " 174 cout <<
"Bin_ns_bh::orbit_omega: central A^2/N^2 : " << asn2 << endl ;
175 cout <<
"Bin_ns_bh::orbit_omega: central d(A^2/N^2)/dX : " 177 cout <<
"Bin_ns_bh::orbit_omega: central N^Y : " << ny << endl ;
178 cout <<
"Bin_ns_bh::orbit_omega: central dN^Y/dX : " << dny << endl ;
179 cout <<
"Bin_ns_bh::orbit_omega: central N.N : " << npn << endl ;
180 cout <<
"Bin_ns_bh::orbit_omega: central d(N.N)/dX : " 201 double omega1 = fact_omeg_min *
omega ;
202 double omega2 = fact_omeg_max *
omega ;
204 cout <<
"Bin_ns_bh::orbit_omega: omega1, omega2 [rad/s] : " 205 << omega1 * f_unit <<
" " << omega2 * f_unit << endl ;
212 zero_list(fonc_bin_ns_bh_orbit, parf, omega1, omega2, nsub,
217 double omeg_min, omeg_max ;
219 cout <<
"Bin_ns_bh:orbit_omega : " << nzer <<
220 "zero(s) found in the interval [omega1, omega2]." << endl ;
221 cout <<
"omega, omega1, omega2 : " <<
omega <<
" " << omega1
222 <<
" " << omega2 << endl ;
223 cout <<
"azer : " << *azer << endl ;
224 cout <<
"bzer : " << *bzer << endl ;
227 cout <<
"Bin_ns_bh::orbit_omega: WARNING : " 228 <<
"no zero detected in the interval" << endl
229 <<
" [" << omega1 * f_unit <<
", " 230 << omega2 * f_unit <<
"] rad/s !" << endl ;
235 double dist_min = fabs(omega2 - omega1) ;
236 int i_dist_min = -1 ;
237 for (
int i=0; i<nzer; i++) {
240 double dist = fabs(
omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
242 if (dist < dist_min) {
247 omeg_min = (*azer)(i_dist_min) ;
248 omeg_max = (*bzer)(i_dist_min) ;
254 cout <<
"Bin_ns_bh:orbit_omega : " 255 <<
"interval selected for the search of the zero : " 256 << endl <<
" [" << omeg_min <<
", " << omeg_max <<
"] = [" 257 << omeg_min * f_unit <<
", " << omeg_max * f_unit <<
"] rad/s " 265 double precis = 1.e-13 ;
267 precis, nitermax, niter) ;
269 cout <<
"Bin_ns_bh::orbit_omega : " 270 <<
"Number of iterations in zerosec for omega : " 273 cout <<
"Bin_ns_bh::orbit_omega : omega [rad/s] : " 274 <<
omega * f_unit << endl ;
282 double fonc_bin_ns_bh_orbit(
double om,
const Param& parf) {
296 double xx = xc -
x_axe ;
302 double bpb = om2 * xx*xx - 2*om * ny * xx + npn ;
304 dphi_cent = ( asn2* ( om* (ny + xx*dny) - om2*xx - 0.5*dnpn )
306 / ( 1 - asn2 * bpb ) ;
309 cout <<
"Bin_ns_bh::orbit_omega : " 310 <<
"It should be the relativistic calculation !" << endl ;
314 return dnulg + dphi_cent ;
Cmp log(const Cmp &)
Neperian logarithm.
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
const Cmp & dsdx() const
Returns of *this , where .
const Map & get_mp() const
Returns the mapping.
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Standard units of space, time and mass.
const Tenseur & get_shift() const
Returns the total shift vector .
Base class for coordinate mappings.
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...
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 get_rot_phi() const
Returns the angle between the x –axis and X –axis.
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
const Tenseur & get_loggam() const
Returns the logarithm of the Lorentz factor between the fluid and the co-orbiting observer...
double x_axe
Absolute X coordinate of the rotation axis.
Et_bin_nsbh star
The neutron star.
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.
const Tenseur & get_nnn() const
Returns the total lapse function N.
double omega
Angular velocity with respect to an asymptotically inertial observer.
const Tenseur & get_confpsi() const
Returns the part of the conformal factor $$.
Cmp pow(const Cmp &, int)
Power .
void orbit_omega(double fact_omeg_min, double fact_omeg_max)
Computes the orbital angular velocity { omega}.
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
int get_taille() const
Gives the total size (ie dim.taille)
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.
bool is_relativistic() const
Returns true for a relativistic star, false for a Newtonian one.
Tensor handling *** DEPRECATED : use class Tensor instead ***.