32char binaire_orbite_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Binaire/binaire_orbite.C,v 1.8 2014/10/24 11:27:49 j_novak Exp $" ;
110#include "utilitaires.h"
114#include "graphique.h"
117double fonc_binaire_axe(
double ,
const Param& ) ;
118double fonc_binaire_orbit(
double ,
const Param& ) ;
131 double dnulg[2], asn2[2], dasn2[2], ny[2], dny[2], npn[2], dnpn[2], xgg[2] ;
132 double nyso[2], dnyso[2], npnso2[2], dnpnso2[2], ori_x[2] ;
134 for (
int i=0; i<2; i++) {
146 Tenseur dln_auto_div = d_logn_auto_div ;
172 cout <<
"Binaire::orbit : unknown value of rot_phi !" << endl ;
177 Cmp tmp = logn_auto_regu + logn_comp + loggam ;
180 dnulg[i] = dln_auto_div(0)(0, 0, 0, 0)
181 + factx * tmp.
dsdx()(0, 0, 0, 0) ;
183 tmp = logn_auto_regu + logn_comp ;
184 cout <<
"dlnndx_div_c : " << dln_auto_div(0)(0, 0, 0, 0) << endl ;
185 cout <<
"dlnndx_c : " << dln_auto_div(0)(0, 0, 0, 0) + factx*tmp.
dsdx()(0, 0, 0, 0) << endl ;
187 cout <<
"dloggamdx_c : " << factx*loggam.
dsdx()(0, 0, 0, 0) << endl ;
189 save_profile(stmp, 0., 10., 0.5*M_PI, 0.,
"prof_logn.d") ;
191 save_profile(stmp, 0., 1.8, 0.5*M_PI, 0.,
"prof_loggam.d") ;
197 double nc = nnn(0, 0, 0, 0) ;
198 double a2c = a_car(0, 0, 0, 0) ;
199 asn2[i] = a2c / (nc * nc) ;
201 if (
et[i]->is_relativistic() ) {
207 double da2c = factx * a_car.
dsdx()(0, 0, 0, 0) ;
208 double dnc = factx * nnn.
dsdx()(0, 0, 0, 0) ;
210 dasn2[i] = ( da2c - 2 * a2c / nc * dnc ) / (nc*nc) ;
216 ny[i] = shift(1)(0, 0, 0, 0) ;
217 nyso[i] = ny[i] /
omega ;
223 dny[i] = factx * shift(1).dsdx()(0, 0, 0, 0) ;
224 dnyso[i] = dny[i] /
omega ;
233 npn[i] = tmp(0, 0, 0, 0) ;
241 dnpn[i] = factx * tmp.
dsdx()(0, 0, 0, 0) ;
257 cout <<
"Binaire::orbit: central d(nu+log(Gam))/dX : "
258 << dnulg[i] << endl ;
259 cout <<
"Binaire::orbit: central A^2/N^2 : " << asn2[i] << endl ;
260 cout <<
"Binaire::orbit: central d(A^2/N^2)/dX : " << dasn2[i] << endl ;
261 cout <<
"Binaire::orbit: central N^Y : " << ny[i] << endl ;
262 cout <<
"Binaire::orbit: central dN^Y/dX : " << dny[i] << endl ;
263 cout <<
"Binaire::orbit: central N.N : " << npn[i] << endl ;
264 cout <<
"Binaire::orbit: central d(N.N)/dX : " << dnpn[i] << endl ;
272 ori_x[i] = (
et[i]->
get_mp()).get_ori_x() ;
286 double ori_x1 = ori_x[0] ;
287 double ori_x2 = ori_x[1] ;
289 if (
et[0]->get_eos() ==
et[1]->get_eos() &&
290 et[0]->get_ent()()(0,0,0,0) ==
et[1]->get_ent()()(0,0,0,0) ) {
316 int nitmax_axe = 200 ;
318 double precis_axe = 1.e-13 ;
320 x_axe =
zerosec(fonc_binaire_axe, paraxe, 0.9*ori_x1, 0.9*ori_x2,
321 precis_axe, nitmax_axe, nit_axe) ;
323 cout <<
"Binaire::orbit : Number of iterations in zerosec for x_axe : "
327 cout <<
"Binaire::orbit : x_axe [km] : " <<
x_axe / km << endl ;
345 double omega1 = fact_omeg_min *
omega ;
346 double omega2 = fact_omeg_max *
omega ;
347 cout <<
"Binaire::orbit: omega1, omega2 [rad/s] : "
348 << omega1 * f_unit <<
" " << omega2 * f_unit << endl ;
355 zero_list(fonc_binaire_orbit, parf, omega1, omega2, nsub,
360 double omeg_min, omeg_max ;
362 cout <<
"Binaire:orbit : " << nzer <<
363 " zero(s) found in the interval [omega1, omega2]." << endl ;
364 cout <<
"omega, omega1, omega2 : " <<
omega <<
" " << omega1
365 <<
" " << omega2 << endl ;
366 cout <<
"azer : " << *azer << endl ;
367 cout <<
"bzer : " << *bzer << endl ;
371 "Binaire::orbit: WARNING : no zero detected in the interval"
372 << endl <<
" [" << omega1 * f_unit <<
", "
373 << omega2 * f_unit <<
"] rad/s !" << endl ;
378 double dist_min = fabs(omega2 - omega1) ;
379 int i_dist_min = -1 ;
380 for (
int i=0; i<nzer; i++) {
383 double dist = fabs(
omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
384 if (dist < dist_min) {
389 omeg_min = (*azer)(i_dist_min) ;
390 omeg_max = (*bzer)(i_dist_min) ;
396 cout <<
"Binaire:orbit : interval selected for the search of the zero : "
397 << endl <<
" [" << omeg_min <<
", " << omeg_max <<
"] = ["
398 << omeg_min * f_unit <<
", " << omeg_max * f_unit <<
"] rad/s " << endl ;
404 double precis = 1.e-13 ;
406 precis, nitermax, niter) ;
408 cout <<
"Binaire::orbit : Number of iterations in zerosec for omega : "
411 cout <<
"Binaire::orbit : omega [rad/s] : "
412 <<
omega * f_unit << endl ;
465 double mass1,
double mass2,
466 double& xgg1,
double& xgg2) {
474 double dnulg[2], asn2[2], dasn2[2], ny[2], dny[2], npn[2], dnpn[2], xgg[2] ;
475 double nyso[2], dnyso[2], npnso2[2], dnpnso2[2], ori_x[2] ;
477 for (
int i=0; i<2; i++) {
489 Tenseur dln_auto_div = d_logn_auto_div ;
515 cout <<
"Binaire::orbit : unknown value of rot_phi !" << endl ;
520 Cmp tmp = logn_auto_regu + logn_comp + loggam ;
523 dnulg[i] = dln_auto_div(0)(0, 0, 0, 0)
524 + factx * tmp.
dsdx()(0, 0, 0, 0) ;
530 double nc = nnn(0, 0, 0, 0) ;
531 double a2c = a_car(0, 0, 0, 0) ;
532 asn2[i] = a2c / (nc * nc) ;
534 if (
et[i]->is_relativistic() ) {
540 double da2c = factx * a_car.
dsdx()(0, 0, 0, 0) ;
541 double dnc = factx * nnn.
dsdx()(0, 0, 0, 0) ;
543 dasn2[i] = ( da2c - 2 * a2c / nc * dnc ) / (nc*nc) ;
549 ny[i] = shift(1)(0, 0, 0, 0) ;
550 nyso[i] = ny[i] /
omega ;
556 dny[i] = factx * shift(1).dsdx()(0, 0, 0, 0) ;
557 dnyso[i] = dny[i] /
omega ;
566 npn[i] = tmp(0, 0, 0, 0) ;
574 dnpn[i] = factx * tmp.
dsdx()(0, 0, 0, 0) ;
590 cout <<
"Binaire::orbit: central d(nu+log(Gam))/dX : "
591 << dnulg[i] << endl ;
592 cout <<
"Binaire::orbit: central A^2/N^2 : " << asn2[i] << endl ;
593 cout <<
"Binaire::orbit: central d(A^2/N^2)/dX : " << dasn2[i] << endl ;
594 cout <<
"Binaire::orbit: central N^Y : " << ny[i] << endl ;
595 cout <<
"Binaire::orbit: central dN^Y/dX : " << dny[i] << endl ;
596 cout <<
"Binaire::orbit: central N.N : " << npn[i] << endl ;
597 cout <<
"Binaire::orbit: central d(N.N)/dX : " << dnpn[i] << endl ;
605 ori_x[i] = (
et[i]->
get_mp()).get_ori_x() ;
619 double ori_x1 = ori_x[0] ;
620 double ori_x2 = ori_x[1] ;
622 if (
et[0]->get_eos() ==
et[1]->get_eos() && mass1 == mass2 ) {
648 int nitmax_axe = 200 ;
650 double precis_axe = 1.e-13 ;
652 x_axe =
zerosec(fonc_binaire_axe, paraxe, 0.9*ori_x1, 0.9*ori_x2,
653 precis_axe, nitmax_axe, nit_axe) ;
655 cout <<
"Binaire::orbit : Number of iterations in zerosec for x_axe : "
659 cout <<
"Binaire::orbit : x_axe [km] : " <<
x_axe / km << endl ;
677 double omega1 = fact_omeg_min *
omega ;
678 double omega2 = fact_omeg_max *
omega ;
679 cout <<
"Binaire::orbit: omega1, omega2 [rad/s] : "
680 << omega1 * f_unit <<
" " << omega2 * f_unit << endl ;
687 zero_list(fonc_binaire_orbit, parf, omega1, omega2, nsub,
692 double omeg_min, omeg_max ;
694 cout <<
"Binaire:orbit : " << nzer <<
695 " zero(s) found in the interval [omega1, omega2]." << endl ;
696 cout <<
"omega, omega1, omega2 : " <<
omega <<
" " << omega1
697 <<
" " << omega2 << endl ;
698 cout <<
"azer : " << *azer << endl ;
699 cout <<
"bzer : " << *bzer << endl ;
703 "Binaire::orbit: WARNING : no zero detected in the interval"
704 << endl <<
" [" << omega1 * f_unit <<
", "
705 << omega2 * f_unit <<
"] rad/s !" << endl ;
710 double dist_min = fabs(omega2 - omega1) ;
711 int i_dist_min = -1 ;
712 for (
int i=0; i<nzer; i++) {
715 double dist = fabs(
omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
716 if (dist < dist_min) {
721 omeg_min = (*azer)(i_dist_min) ;
722 omeg_max = (*bzer)(i_dist_min) ;
728 cout <<
"Binaire:orbit : interval selected for the search of the zero : "
729 << endl <<
" [" << omeg_min <<
", " << omeg_max <<
"] = ["
730 << omeg_min * f_unit <<
", " << omeg_max * f_unit <<
"] rad/s " << endl ;
736 double precis = 1.e-13 ;
738 precis, nitermax, niter) ;
740 cout <<
"Binaire::orbit : Number of iterations in zerosec for omega : "
743 cout <<
"Binaire::orbit : omega [rad/s] : "
744 <<
omega * f_unit << endl ;
800double fonc_binaire_axe(
double x_rot,
const Param& paraxe) {
824 double x1 = ori_x1 - x_rot ;
825 double x2 = ori_x2 - x_rot ;
829 double andan_1 = 0.5 * dasn2_1 + asn2_1 * dnulg_1 ;
830 double andan_2 = 0.5 * dasn2_2 + asn2_2 * dnulg_2 ;
832 double bpb_1 = x1 * x1 - 2. * nyso_1 * x1 + npnso2_1 ;
833 double bpb_2 = x2 * x2 - 2. * nyso_2 * x2 + npnso2_2 ;
835 double cpc_1 = 0.5 * dnpnso2_1 + x1 * (1. - dnyso_1) - nyso_1 ;
836 double cpc_2 = 0.5 * dnpnso2_2 + x2 * (1. - dnyso_2) - nyso_2 ;
838 om2_star1 = dnulg_1 / (andan_1 * bpb_1 + asn2_1 * cpc_1) ;
839 om2_star2 = dnulg_2 / (andan_2 * bpb_2 + asn2_2 * cpc_2) ;
844 om2_star1 = dnulg_1 / x1 ;
845 om2_star2 = dnulg_2 / x2 ;
849 return om2_star1 - om2_star2 ;
857double fonc_binaire_orbit(
double om,
const Param& parf) {
859 int relat = parf.get_int() ;
861 double xc = parf.get_double(0) ;
862 double dnulg = parf.get_double(1) ;
863 double asn2 = parf.get_double(2) ;
864 double dasn2 = parf.get_double(3) ;
865 double ny = parf.get_double(4) ;
866 double dny = parf.get_double(5) ;
867 double npn = parf.get_double(6) ;
868 double dnpn = parf.get_double(7) ;
869 double x_axe = parf.get_double(8) ;
871 double xx = xc - x_axe ;
877 double bpb = om2 * xx*xx - 2*om * ny * xx + npn ;
879 dphi_cent = ( asn2* ( om* (ny + xx*dny) - om2*xx - 0.5*dnpn )
881 / ( 1 - asn2 * bpb ) ;
884 dphi_cent = - om2*xx ;
887 return dnulg + dphi_cent ;
Etoile_bin * et[2]
Array of the two stars (to perform loops on the stars): {\tt et[0]} contains the address of {\tt star...
void orbit_eqmass(double fact_omeg_min, double fact_omeg_max, double mass1, double mass2, double &xgg1, double &xgg2)
Computes the orbital angular velocity {\tt omega} and the position of the rotation axis {\tt x_axe}.
const Base_vect_cart ref_triad
Cartesian triad of the absolute reference frame.
double omega
Angular velocity with respect to an asymptotically inertial observer.
void orbit(double fact_omeg_min, double fact_omeg_max, double &xgg1, double &xgg2)
Computes the orbital angular velocity {\tt omega} and the position of the rotation axis {\tt x_axe}.
double x_axe
Absolute X coordinate of the rotation axis.
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
const Cmp & dsdx() const
Returns of *this , where .
const Tenseur & get_logn_comp() const
Returns the part of the lapse logarithm (gravitational potential at the Newtonian limit) generated pr...
const Tenseur & get_loggam() const
Returns the logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
virtual double xa_barycenter() const
Absolute coordinate X of the barycenter of the baryon density, defined according to the formula.
const Tenseur & get_nnn() const
Returns the total lapse function N.
const Tenseur & get_shift() const
Returns the total shift vector .
const Tenseur & get_logn_auto_regu() const
Returns the regular part of the logarithm of the part of the lapse N generated principaly by the star...
const Map & get_mp() const
Returns the mapping.
bool is_relativistic() const
Returns true for a relativistic star, false for a Newtonian one.
const Tenseur & get_d_logn_auto_div() const
Returns the gradient of logn_auto_div.
const Tenseur & get_a_car() const
Returns the total conformal factor .
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.
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).
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 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.
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.