27char binary_orbit_xcts_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Binary_xcts/binary_orbit_xcts.C,v 1.14 2014/10/24 14:10:24 j_novak Exp $" ;
83#include "binary_xcts.h"
86#include "utilitaires.h"
92double fonc_binary_xcts_axe(
double ,
const Param& ) ;
93double fonc_binary_xcts_orbit(
double ,
const Param& ) ;
108 double dnulg[2], asn2[2], dasn2[2], ny[2], dny[2], npn[2], dnpn[2], xgg[2] ;
109 double nyso[2], dnyso[2], npnso2[2], dnpnso2[2], ori_x[2] ;
111 for (
int i=0; i<2; i++) {
146 cout <<
"Binary_xcts::orbit : unknown value of rot_phi !" << endl ;
155 Scalar tgraph = logn -
log( (1. +
et[i]->get_chi_auto()) / (1. +
et[i]->get_Psi_auto()) ) ;
158 save_profile(tgraph, 0., 10., 0.5*M_PI, 0.,
"prof_logn.d") ;
159 save_profile(
et[i]->get_loggam(), 0., 1.8, 0.5*M_PI, 0.,
"prof_loggam.d") ;
165 Scalar Psi6schi2 =
pow(Psi, 6)/(chi % chi) ;
179 ny[i] = factx*shift(2).val_grid_point(0, 0, 0, 0) ;
181 nyso[i] = ny[i] /
omega ;
187 dny[i] = shift(2).dsdx().val_grid_point(0, 0, 0, 0) ;
189 dnyso[i] = dny[i] /
omega ;
209 cout <<
"Binary_xcts::orbit: central d(nu+log(Gam))/dX : "
210 << dnulg[i] << endl ;
211 cout <<
"Binary_xcts::orbit: central A^2/N^2 : "
213 cout <<
"Binary_xcts::orbit: central d(A^2/N^2)/dX : "
214 << dasn2[i] << endl ;
215 cout <<
"Binary_xcts::orbit: central N^Y : "
217 cout <<
"Binary_xcts::orbit: central dN^Y/dX : "
219 cout <<
"Binary_xcts::orbit: central N.N : "
221 cout <<
"Binary_xcts::orbit: central d(N.N)/dX : "
225 ori_x[i] = (
et[i]->
get_mp()).get_ori_x() ;
239 double ori_x1 = ori_x[0] ;
240 double ori_x2 = ori_x[1] ;
243 fabs(
et[0]->get_ent().val_grid_point(0,0,0,0) -
244 et[1]->get_ent().val_grid_point(0,0,0,0) ) < 1.e-14 ) {
269 int nitmax_axe = 200 ;
271 double precis_axe = 1.e-13 ;
273 x_axe =
zerosec(fonc_binary_xcts_axe, paraxe, 0.9*ori_x1, 0.9*ori_x2,
274 precis_axe, nitmax_axe, nit_axe) ;
276 cout <<
"Binary_xcts::orbit : Number of iterations in zerosec for x_axe : "
280 cout <<
"Binary_xcts::orbit: x_axe [km] : " <<
x_axe / km << endl ;
298 double omega1 = fact_omeg_min *
omega ;
299 double omega2 = fact_omeg_max *
omega ;
300 cout <<
"Binary_xcts::orbit: omega1, omega2 [rad/s] : "
301 << omega1 * f_unit <<
" " << omega2 * f_unit << endl ;
308 zero_list(fonc_binary_xcts_orbit, parf, omega1, omega2, nsub,
313 double omeg_min, omeg_max ;
315 cout <<
"Binary_xcts:orbit : " << nzer <<
316 " zero(s) found in the interval [omega1, omega2]." << endl ;
317 cout <<
"omega, omega1, omega2 : " <<
omega <<
" " << omega1
318 <<
" " << omega2 << endl ;
319 cout <<
"azer : " << *azer << endl ;
320 cout <<
"bzer : " << *bzer << endl ;
324 "Binary_xcts::orbit: WARNING : no zero detected in the interval"
325 << endl <<
" [" << omega1 * f_unit <<
", "
326 << omega2 * f_unit <<
"] rad/s !" << endl ;
331 double dist_min = fabs(omega2 - omega1) ;
332 int i_dist_min = -1 ;
333 for (
int i=0; i<nzer; i++) {
336 double dist = fabs(
omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
337 if (dist < dist_min) {
342 omeg_min = (*azer)(i_dist_min) ;
343 omeg_max = (*bzer)(i_dist_min) ;
349 cout <<
"Binary_xcts::orbit : interval selected for the search of the zero : "
350 << endl <<
" [" << omeg_min <<
", " << omeg_max <<
"] = ["
351 << omeg_min * f_unit <<
", " << omeg_max * f_unit <<
"] rad/s " << endl ;
357 double precis = 1.e-13 ;
358 omega =
zerosec_b(fonc_binary_xcts_orbit, parf, omeg_min, omeg_max,
359 precis, nitermax, niter) ;
361 cout <<
"Binary_xcts::orbit : Number of iterations in zerosec for omega : "
364 cout <<
"Binary_xcts::orbit : omega [rad/s] : "
365 <<
omega * f_unit << endl ;
373double fonc_binary_xcts_axe(
double x_rot,
const Param& paraxe) {
397 double x1 = ori_x1 - x_rot ;
398 double x2 = ori_x2 - x_rot ;
402 double andan_1 = 0.5 * dasn2_1 + asn2_1 * dnulg_1 ;
403 double andan_2 = 0.5 * dasn2_2 + asn2_2 * dnulg_2 ;
405 double bpb_1 = x1 * x1 - 2. * nyso_1 * x1 + npnso2_1 ;
406 double bpb_2 = x2 * x2 - 2. * nyso_2 * x2 + npnso2_2 ;
408 double cpc_1 = 0.5 * dnpnso2_1 + x1 * (1. - dnyso_1) - nyso_1 ;
409 double cpc_2 = 0.5 * dnpnso2_2 + x2 * (1. - dnyso_2) - nyso_2 ;
411 om2_star1 = dnulg_1 / (andan_1 * bpb_1 + asn2_1 * cpc_1) ;
412 om2_star2 = dnulg_2 / (andan_2 * bpb_2 + asn2_2 * cpc_2) ;
417 om2_star1 = dnulg_1 / x1 ;
418 om2_star2 = dnulg_2 / x2 ;
422 return om2_star1 - om2_star2 ;
430double fonc_binary_xcts_orbit(
double om,
const Param& parf) {
432 int relat = parf.get_int() ;
434 double xc = parf.get_double(0) ;
435 double dnulg = parf.get_double(1) ;
436 double asn2 = parf.get_double(2) ;
437 double dasn2 = parf.get_double(3) ;
438 double ny = parf.get_double(4) ;
439 double dny = parf.get_double(5) ;
440 double npn = parf.get_double(6) ;
441 double dnpn = parf.get_double(7) ;
442 double x_axe = parf.get_double(8) ;
444 double xx = xc - x_axe ;
450 double bpb = om2 * xx*xx - 2*om * ny * xx + npn ;
452 dphi_cent = ( asn2* ( om* (ny + xx*dny) - om2*xx - 0.5*dnpn )
454 / ( 1 - asn2 * bpb ) ;
457 dphi_cent = - om2*xx ;
460 return dnulg + dphi_cent ;
double x_axe
Absolute X coordinate of the rotation axis.
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.
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]...
double omega
Angular velocity with respect to an asymptotically inertial observer.
Base class for coordinate mappings.
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Metric for tensor calculation.
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.
Tensor field of valence 0 (or component of a tensorial field).
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
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 & dsdx() const
Returns of *this , where .
virtual double xa_barycenter() const
Absolute coordinate X of the barycenter of the baryon density,.
const Scalar & get_chi_auto() const
Returns the scalar field generated principally by the star.
const Scalar & get_chi_comp() const
Returns the scalar field generated principally by the companion star.
const Scalar & get_Psi_comp() const
Returns the scalar field generated principally by the companion star.
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_loggam() const
Returns the logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
const Scalar & get_Psi_auto() const
Returns the scalar field generated principally by the star.
const Map & get_mp() const
Returns the mapping.
const Eos & get_eos() const
Returns the equation of state.
const Vector & get_beta() const
Returns the shift vector .
int get_taille() const
Gives the total size (ie dim.taille)
Tensor field of valence 1.
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Cmp pow(const Cmp &, int)
Power .
Cmp log(const Cmp &)
Neperian logarithm.
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 .
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.
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.
Tensor up_down(const Metric &gam) const
Computes a new tensor by raising or lowering all the indices of *this .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Standard units of space, time and mass.