29char bin_ns_bh_orbit_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Bin_ns_bh/bin_ns_bh_orbit.C,v 1.7 2014/10/24 14:10:24 j_novak Exp $" ;
67#include "utilitaires.h"
71double fonc_bin_ns_bh_orbit(
double ,
const Param& ) ;
83 double dnulg, asn2, dasn2, ny, dny, npn, dnpn ;
92 Cmp confpsi_q =
pow(confpsi, 4.) ;
110 cout <<
"Bin_ns_bh::orbit_omega : unknown value of rot_phi !"
116 Cmp tmp =
log( nnn ) + loggam ;
120 dnulg = factx * tmp.
dsdx()(0, 0, 0, 0) ;
125 double nc = nnn(0, 0, 0, 0) ;
126 double a2c = confpsi_q(0, 0, 0, 0) ;
127 asn2 = a2c / (nc * nc) ;
134 double da2c = factx * confpsi_q.
dsdx()(0, 0, 0, 0) ;
135 double dnc = factx * nnn.
dsdx()(0, 0, 0, 0) ;
137 dasn2 = ( da2c - 2 * a2c / nc * dnc ) / (nc*nc) ;
142 ny = shift(1)(0, 0, 0, 0) ;
147 dny = factx * shift(1).dsdx()(0, 0, 0, 0) ;
154 npn = tmp(0, 0, 0, 0) ;
160 dnpn = factx * tmp.
dsdx()(0, 0, 0, 0) ;
164 cout <<
"Bin_ns_bh::orbit_omega : "
165 <<
"It should be the relativistic calculation !" << endl ;
169 cout <<
"Bin_ns_bh::orbit_omega: central d(nu+log(Gam))/dX : "
171 cout <<
"Bin_ns_bh::orbit_omega: central A^2/N^2 : " << asn2 << endl ;
172 cout <<
"Bin_ns_bh::orbit_omega: central d(A^2/N^2)/dX : "
174 cout <<
"Bin_ns_bh::orbit_omega: central N^Y : " << ny << endl ;
175 cout <<
"Bin_ns_bh::orbit_omega: central dN^Y/dX : " << dny << endl ;
176 cout <<
"Bin_ns_bh::orbit_omega: central N.N : " << npn << endl ;
177 cout <<
"Bin_ns_bh::orbit_omega: central d(N.N)/dX : "
198 double omega1 = fact_omeg_min *
omega ;
199 double omega2 = fact_omeg_max *
omega ;
201 cout <<
"Bin_ns_bh::orbit_omega: omega1, omega2 [rad/s] : "
202 << omega1 * f_unit <<
" " << omega2 * f_unit << endl ;
209 zero_list(fonc_bin_ns_bh_orbit, parf, omega1, omega2, nsub,
214 double omeg_min, omeg_max ;
216 cout <<
"Bin_ns_bh:orbit_omega : " << nzer <<
217 "zero(s) found in the interval [omega1, omega2]." << endl ;
218 cout <<
"omega, omega1, omega2 : " <<
omega <<
" " << omega1
219 <<
" " << omega2 << endl ;
220 cout <<
"azer : " << *azer << endl ;
221 cout <<
"bzer : " << *bzer << endl ;
224 cout <<
"Bin_ns_bh::orbit_omega: WARNING : "
225 <<
"no zero detected in the interval" << endl
226 <<
" [" << omega1 * f_unit <<
", "
227 << omega2 * f_unit <<
"] rad/s !" << endl ;
232 double dist_min = fabs(omega2 - omega1) ;
233 int i_dist_min = -1 ;
234 for (
int i=0; i<nzer; i++) {
237 double dist = fabs(
omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
239 if (dist < dist_min) {
244 omeg_min = (*azer)(i_dist_min) ;
245 omeg_max = (*bzer)(i_dist_min) ;
251 cout <<
"Bin_ns_bh:orbit_omega : "
252 <<
"interval selected for the search of the zero : "
253 << endl <<
" [" << omeg_min <<
", " << omeg_max <<
"] = ["
254 << omeg_min * f_unit <<
", " << omeg_max * f_unit <<
"] rad/s "
262 double precis = 1.e-13 ;
264 precis, nitermax, niter) ;
266 cout <<
"Bin_ns_bh::orbit_omega : "
267 <<
"Number of iterations in zerosec for omega : "
270 cout <<
"Bin_ns_bh::orbit_omega : omega [rad/s] : "
271 <<
omega * f_unit << endl ;
279double fonc_bin_ns_bh_orbit(
double om,
const Param& parf) {
293 double xx = xc - x_axe ;
299 double bpb = om2 * xx*xx - 2*om * ny * xx + npn ;
301 dphi_cent = ( asn2* ( om* (ny + xx*dny) - om2*xx - 0.5*dnpn )
303 / ( 1 - asn2 * bpb ) ;
306 cout <<
"Bin_ns_bh::orbit_omega : "
307 <<
"It should be the relativistic calculation !" << endl ;
311 return dnulg + dphi_cent ;
double x_axe
Absolute X coordinate of the rotation axis.
Et_bin_nsbh star
The neutron star.
double omega
Angular velocity with respect to an asymptotically inertial observer.
void orbit_omega(double fact_omeg_min, double fact_omeg_max)
Computes the orbital angular velocity {\tt omega}.
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
const Cmp & dsdx() const
Returns of *this , where .
const Tenseur & get_confpsi() const
Returns the part of the conformal factor $\Psi$.
const Tenseur & get_loggam() const
Returns the logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
const Tenseur & get_nnn() const
Returns the total lapse function N.
const Tenseur & get_shift() const
Returns the total shift vector .
const Map & get_mp() const
Returns the mapping.
bool is_relativistic() const
Returns true for a relativistic star, false for a Newtonian one.
Base class for coordinate mappings.
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
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 ***.
Cmp pow(const Cmp &, int)
Power .
Cmp log(const Cmp &)
Neperian logarithm.
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...
Standard units of space, time and mass.