86 #include "binary_xcts.h" 89 #include "utilitaires.h" 92 #include "graphique.h" 95 double fonc_binary_xcts_axe(
double ,
const Param& ) ;
96 double fonc_binary_xcts_orbit(
double ,
const Param& ) ;
101 double fact_omeg_max,
111 double dnulg[2], asn2[2], dasn2[2], ny[2], dny[2], npn[2], dnpn[2], xgg[2] ;
112 double nyso[2], dnyso[2], npnso2[2], dnpnso2[2], ori_x[2] ;
114 for (
int i=0; i<2; i++) {
149 cout <<
"Binary_xcts::orbit : unknown value of rot_phi !" << endl ;
158 Scalar tgraph = logn -
log( (1. +
et[i]->get_chi_auto()) / (1. +
et[i]->get_Psi_auto()) ) ;
161 save_profile(tgraph, 0., 10., 0.5*M_PI, 0.,
"prof_logn.d") ;
162 save_profile(
et[i]->get_loggam(), 0., 1.8, 0.5*M_PI, 0.,
"prof_loggam.d") ;
168 Scalar Psi6schi2 =
pow(Psi, 6)/(chi % chi) ;
182 ny[i] = factx*shift(2).val_grid_point(0, 0, 0, 0) ;
184 nyso[i] = ny[i] /
omega ;
190 dny[i] = shift(2).dsdx().val_grid_point(0, 0, 0, 0) ;
192 dnyso[i] = dny[i] /
omega ;
212 cout <<
"Binary_xcts::orbit: central d(nu+log(Gam))/dX : " 213 << dnulg[i] << endl ;
214 cout <<
"Binary_xcts::orbit: central A^2/N^2 : " 216 cout <<
"Binary_xcts::orbit: central d(A^2/N^2)/dX : " 217 << dasn2[i] << endl ;
218 cout <<
"Binary_xcts::orbit: central N^Y : " 220 cout <<
"Binary_xcts::orbit: central dN^Y/dX : " 222 cout <<
"Binary_xcts::orbit: central N.N : " 224 cout <<
"Binary_xcts::orbit: central d(N.N)/dX : " 228 ori_x[i] = (
et[i]->
get_mp()).get_ori_x() ;
242 double ori_x1 = ori_x[0] ;
243 double ori_x2 = ori_x[1] ;
246 fabs(
et[0]->get_ent().val_grid_point(0,0,0,0) -
247 et[1]->get_ent().val_grid_point(0,0,0,0) ) < 1.e-14 ) {
272 int nitmax_axe = 200 ;
274 double precis_axe = 1.e-13 ;
276 x_axe =
zerosec(fonc_binary_xcts_axe, paraxe, 0.9*ori_x1, 0.9*ori_x2,
277 precis_axe, nitmax_axe, nit_axe) ;
279 cout <<
"Binary_xcts::orbit : Number of iterations in zerosec for x_axe : " 283 cout <<
"Binary_xcts::orbit: x_axe [km] : " <<
x_axe / km << endl ;
301 double omega1 = fact_omeg_min *
omega ;
302 double omega2 = fact_omeg_max *
omega ;
303 cout <<
"Binary_xcts::orbit: omega1, omega2 [rad/s] : " 304 << omega1 * f_unit <<
" " << omega2 * f_unit << endl ;
311 zero_list(fonc_binary_xcts_orbit, parf, omega1, omega2, nsub,
316 double omeg_min, omeg_max ;
318 cout <<
"Binary_xcts:orbit : " << nzer <<
319 " zero(s) found in the interval [omega1, omega2]." << endl ;
320 cout <<
"omega, omega1, omega2 : " <<
omega <<
" " << omega1
321 <<
" " << omega2 << endl ;
322 cout <<
"azer : " << *azer << endl ;
323 cout <<
"bzer : " << *bzer << endl ;
327 "Binary_xcts::orbit: WARNING : no zero detected in the interval" 328 << endl <<
" [" << omega1 * f_unit <<
", " 329 << omega2 * f_unit <<
"] rad/s !" << endl ;
334 double dist_min = fabs(omega2 - omega1) ;
335 int i_dist_min = -1 ;
336 for (
int i=0; i<nzer; i++) {
339 double dist = fabs(
omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
340 if (dist < dist_min) {
345 omeg_min = (*azer)(i_dist_min) ;
346 omeg_max = (*bzer)(i_dist_min) ;
352 cout <<
"Binary_xcts::orbit : interval selected for the search of the zero : " 353 << endl <<
" [" << omeg_min <<
", " << omeg_max <<
"] = [" 354 << omeg_min * f_unit <<
", " << omeg_max * f_unit <<
"] rad/s " << endl ;
360 double precis = 1.e-13 ;
361 omega =
zerosec_b(fonc_binary_xcts_orbit, parf, omeg_min, omeg_max,
362 precis, nitermax, niter) ;
364 cout <<
"Binary_xcts::orbit : Number of iterations in zerosec for omega : " 367 cout <<
"Binary_xcts::orbit : omega [rad/s] : " 368 <<
omega * f_unit << endl ;
376 double fonc_binary_xcts_axe(
double x_rot,
const Param& paraxe) {
400 double x1 = ori_x1 - x_rot ;
401 double x2 = ori_x2 - x_rot ;
405 double andan_1 = 0.5 * dasn2_1 + asn2_1 * dnulg_1 ;
406 double andan_2 = 0.5 * dasn2_2 + asn2_2 * dnulg_2 ;
408 double bpb_1 = x1 * x1 - 2. * nyso_1 * x1 + npnso2_1 ;
409 double bpb_2 = x2 * x2 - 2. * nyso_2 * x2 + npnso2_2 ;
411 double cpc_1 = 0.5 * dnpnso2_1 + x1 * (1. - dnyso_1) - nyso_1 ;
412 double cpc_2 = 0.5 * dnpnso2_2 + x2 * (1. - dnyso_2) - nyso_2 ;
414 om2_star1 = dnulg_1 / (andan_1 * bpb_1 + asn2_1 * cpc_1) ;
415 om2_star2 = dnulg_2 / (andan_2 * bpb_2 + asn2_2 * cpc_2) ;
420 om2_star1 = dnulg_1 / x1 ;
421 om2_star2 = dnulg_2 / x2 ;
425 return om2_star1 - om2_star2 ;
433 double fonc_binary_xcts_orbit(
double om,
const Param& parf) {
435 int relat = parf.get_int() ;
437 double xc = parf.get_double(0) ;
438 double dnulg = parf.get_double(1) ;
439 double asn2 = parf.get_double(2) ;
440 double dasn2 = parf.get_double(3) ;
441 double ny = parf.get_double(4) ;
442 double dny = parf.get_double(5) ;
443 double npn = parf.get_double(6) ;
444 double dnpn = parf.get_double(7) ;
445 double x_axe = parf.get_double(8) ;
447 double xx = xc -
x_axe ;
453 double bpb = om2 * xx*xx - 2*om * ny * xx + npn ;
455 dphi_cent = ( asn2* ( om* (ny + xx*dny) - om2*xx - 0.5*dnpn )
457 / ( 1 - asn2 * bpb ) ;
460 dphi_cent = - om2*xx ;
463 return dnulg + dphi_cent ;
Metric for tensor calculation.
Cmp log(const Cmp &)
Neperian logarithm.
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
const Vector & get_beta() const
Returns the shift vector .
const Scalar & get_loggam() const
Returns the logarithm of the Lorentz factor between the fluid and the co-orbiting observer...
Star_bin_xcts * et[2]
Array of the two stars (to perform loops on the stars): et[0] contains the address of star1 and et[1]...
Standard units of space, time and mass.
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...
Tensor field of valence 0 (or component of a tensorial field).
Base class for coordinate mappings.
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
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 & dsdx() const
Returns of *this , where .
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
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.
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
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 int & get_int(int position=0) const
Returns the reference of a int stored in the list.
double x_axe
Absolute X coordinate of the rotation axis.
const Scalar & get_Psi_auto() const
Returns the scalar field generated principally by the 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 Metric & get_flat() const
Return the flat metric defined on the mapping (Spherical components with respect to the mapping of th...
const Scalar & get_chi_comp() const
Returns the scalar field generated principally by the companion star.
double omega
Angular velocity with respect to an asymptotically inertial observer.
Cmp pow(const Cmp &, int)
Power .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
const Scalar & get_Psi_comp() const
Returns the scalar field generated principally by the companion star.
void save_profile(const Scalar &uu, double r_min, double r_max, double theta, double phi, const char *filename)
Saves in a file the profile of a Scalar along some radial axis determined by a fixed value of ...
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
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.
virtual double xa_barycenter() const
Absolute coordinate X of the barycenter of the baryon density,.
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_chi_auto() const
Returns the scalar field generated principally by the star.