binary_orbite.C
/*
* Method of class Binary to compute the orbital angular velocity {\tt omega}
* and the position of the rotation axis {\tt x_axe}.
*
* (See file binary.h for documentation)
*
*/
/*
* Copyright (c) 2004 Francois Limousin
*
* This file is part of LORENE.
*
* LORENE is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* LORENE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with LORENE; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
char binary_orbite_C[] = "$Header: /cvsroot/Lorene/C++/Source/Binary/binary_orbite.C,v 1.9 2015/08/10 15:32:26 j_novak Exp $" ;
/*
* $Id: binary_orbite.C,v 1.9 2015/08/10 15:32:26 j_novak Exp $
* $Log: binary_orbite.C,v $
* Revision 1.9 2015/08/10 15:32:26 j_novak
* Better calls to Param::add_int(), to avoid weird problems (e.g. with g++ 4.8).
*
* Revision 1.8 2014/10/13 08:52:45 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.7 2014/10/06 15:12:59 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.6 2005/09/13 19:38:31 f_limousin
* Reintroduction of the resolution of the equations in cartesian coordinates.
*
* Revision 1.5 2005/02/17 17:35:35 f_limousin
* Change the name of some quantities to be consistent with other classes
* (for instance nnn is changed to nn, shift to beta, beta to lnq...)
*
* Revision 1.4 2004/03/25 10:29:02 j_novak
* All LORENE's units are now defined in the namespace Unites (in file unites.h).
*
* Revision 1.3 2004/02/24 12:39:30 f_limousin
* Change fonc_bin_ncp_orbit to fonc_binary_orbit and fonc_bin_ncp_axe
* to fonc_binary_axe.
*
* Revision 1.2 2004/02/21 17:05:13 e_gourgoulhon
* Method Scalar::point renamed Scalar::val_grid_point.
* Method Scalar::set_point renamed Scalar::set_grid_point.
*
* Revision 1.1 2004/01/20 15:22:19 f_limousin
* First version
*
*
* $Header: /cvsroot/Lorene/C++/Source/Binary/binary_orbite.C,v 1.9 2015/08/10 15:32:26 j_novak Exp $
*
*/
// Headers C
#include <cmath>
// Headers Lorene
#include "binary.h"
#include "eos.h"
#include "param.h"
#include "utilitaires.h"
#include "unites.h"
namespace Lorene {
double fonc_binary_axe(double , const Param& ) ;
double fonc_binary_orbit(double , const Param& ) ;
//******************************************************************************
void Binary::orbit(double fact_omeg_min, double fact_omeg_max, double& xgg1,
double& xgg2) {
using namespace Unites ;
//-------------------------------------------------------------
// Evaluation of various quantities at the center of each star
//-------------------------------------------------------------
double g00[2], g10[2], g20[2], g11[2], g21[2], g22[2], bx[2], by[2] ;
double bz[2], d1sn2[2], unsn2[2] ;
double dnulg[2], xgg[2], ori_x[2], dg00[2], dg10[2], dg20[2], dg11[2] ;
double dg21[2], dg22[2], dbx[2], dby[2], dbz[2], dbymo[2] ;
for (int i=0; i<2; i++) {
const Scalar& logn_auto = et[i]->get_logn_auto() ;
const Scalar& logn_comp = et[i]->get_logn_comp() ;
const Scalar& loggam = et[i]->get_loggam() ;
const Scalar& nn = et[i]->get_nn() ;
Vector shift = et[i]->get_beta() ;
const Metric& gamma = et[i]->get_gamma() ;
Tensor gamma_cov = gamma.cov() ;
// With the new convention for shift (beta^i = - N^i)
shift = - shift ;
// All tensors must be in the cartesian triad
shift.change_triad(et[i]->mp.get_bvect_cart()) ;
gamma_cov.change_triad(et[i]->mp.get_bvect_cart()) ;
const Scalar& gg00 = gamma_cov(1,1) ;
const Scalar& gg10 = gamma_cov(2,1) ;
const Scalar& gg20 = gamma_cov(3,1) ;
const Scalar& gg11 = gamma_cov(2,2) ;
const Scalar& gg21 = gamma_cov(3,2) ;
const Scalar& gg22 = gamma_cov(3,3) ;
const Scalar& bbx = shift(1) ;
const Scalar& bby = shift(2) ;
const Scalar& bbz = shift(3) ;
cout << "gg00" << endl << norme(gg00) << endl ;
cout << "gg10" << endl << norme(gg10) << endl ;
cout << "gg20" << endl << norme(gg20) << endl ;
cout << "gg11" << endl << norme(gg11) << endl ;
cout << "gg21" << endl << norme(gg21) << endl ;
cout << "gg22" << endl << norme(gg22) << endl ;
cout << "bbx" << endl << norme(bbx) << endl ;
cout << "bby" << endl << norme(bby) << endl ;
cout << "bbz" << endl << norme(bbz) << endl ;
//----------------------------------
// Calcul de d/dX( nu + ln(Gamma) ) au centre de l'etoile ---> dnulg[i]
//----------------------------------
Scalar tmp = logn_auto + logn_comp + loggam ;
cout << "logn" << endl << norme(logn_auto + logn_comp) << endl ;
cout << "loggam" << endl << norme(loggam) << endl ;
cout << "dnulg" << endl << norme(tmp.dsdx()) << endl ;
// ... gradient suivant X :
dnulg[i] = tmp.dsdx().val_grid_point(0, 0, 0, 0) ;
cout.precision(8) ;
cout << "dnulg = " << dnulg[i] << endl ;
//----------------------------------
// Calcul de gij, lapse et shift au centre de l'etoile
//----------------------------------
g00[i] = gg00.val_grid_point(0,0,0,0) ;
g10[i] = gg10.val_grid_point(0,0,0,0) ;
g20[i] = gg20.val_grid_point(0,0,0,0) ;
g11[i] = gg11.val_grid_point(0,0,0,0) ;
g21[i] = gg21.val_grid_point(0,0,0,0) ;
g22[i] = gg22.val_grid_point(0,0,0,0) ;
bx[i] = bbx.val_grid_point(0,0,0,0) ;
by[i] = bby.val_grid_point(0,0,0,0) ;
bz[i] = bbz.val_grid_point(0,0,0,0) ;
unsn2[i] = 1/(nn.val_grid_point(0,0,0,0)*nn.val_grid_point(0,0,0,0)) ;
//----------------------------------
// Calcul de d/dX(gij), d/dX(shift) au centre de l'etoile
//----------------------------------
dg00[i] = gg00.dsdx().val_grid_point(0,0,0,0) ;
dg10[i] = gg10.dsdx().val_grid_point(0,0,0,0) ;
dg20[i] = gg20.dsdx().val_grid_point(0,0,0,0) ;
dg11[i] = gg11.dsdx().val_grid_point(0,0,0,0) ;
dg21[i] = gg21.dsdx().val_grid_point(0,0,0,0) ;
dg22[i] = gg22.dsdx().val_grid_point(0,0,0,0) ;
dbx[i] = bbx.dsdx().val_grid_point(0,0,0,0) ;
dby[i] = bby.dsdx().val_grid_point(0,0,0,0) ;
dbz[i] = bbz.dsdx().val_grid_point(0,0,0,0) ;
dbymo[i] = bby.dsdx().val_grid_point(0,0,0,0) - omega ;
d1sn2[i] = (1/(nn*nn)).dsdx().val_grid_point(0,0,0,0) ;
cout << "Binary::orbit: central d(nu+log(Gam))/dX : "
<< dnulg[i] << endl ;
cout << "Binary::orbit: central g00 :" << g00[i] << endl ;
cout << "Binary::orbit: central g10 :" << g10[i] << endl ;
cout << "Binary::orbit: central g20 :" << g20[i] << endl ;
cout << "Binary::orbit: central g11 :" << g11[i] << endl ;
cout << "Binary::orbit: central g21 :" << g21[i] << endl ;
cout << "Binary::orbit: central g22 :" << g22[i] << endl ;
cout << "Binary::orbit: central shift_x :" << bx[i] << endl ;
cout << "Binary::orbit: central shift_y :" << by[i] << endl ;
cout << "Binary::orbit: central shift_z :" << bz[i] << endl ;
cout << "Binary::orbit: central d/dX(g00) :" << dg00[i] << endl ;
cout << "Binary::orbit: central d/dX(g10) :" << dg10[i] << endl ;
cout << "Binary::orbit: central d/dX(g20) :" << dg20[i] << endl ;
cout << "Binary::orbit: central d/dX(g11) :" << dg11[i] << endl ;
cout << "Binary::orbit: central d/dX(g21) :" << dg21[i] << endl ;
cout << "Binary::orbit: central d/dX(g22) :" << dg22[i] << endl ;
cout << "Binary::orbit: central d/dX(shift_x) :" << dbx[i] << endl ;
cout << "Binary::orbit: central d/dX(shift_y) :" << dby[i] << endl ;
cout << "Binary::orbit: central d/dX(shift_z) :" << dbz[i] << endl ;
//----------------------
// Pour information seulement : 1/ calcul des positions des "centres de
// de masse"
// 2/ calcul de dH/dX en r=0
//-----------------------
ori_x[i] = (et[i]->get_mp()).get_ori_x() ;
xgg[i] = (et[i]->xa_barycenter() - ori_x[i]) ;
} // fin de la boucle sur les etoiles
xgg1 = xgg[0] ;
xgg2 = xgg[1] ;
//---------------------------------
// Position de l'axe de rotation
//---------------------------------
double ori_x1 = ori_x[0] ;
double ori_x2 = ori_x[1] ;
if ( et[0]->get_eos() == et[1]->get_eos() &&
et[0]->get_ent().val_grid_point(0,0,0,0) ==
et[1]->get_ent().val_grid_point(0,0,0,0) ) {
x_axe = 0. ;
}
else {
Param paraxe ;
paraxe.add_double( ori_x1, 0) ;
paraxe.add_double( ori_x2, 1) ;
paraxe.add_double( dnulg[0], 2) ;
paraxe.add_double( dnulg[1], 3) ;
paraxe.add_double( g00[0], 4) ;
paraxe.add_double( g00[1], 5) ;
paraxe.add_double( g10[0], 6) ;
paraxe.add_double( g10[1], 7) ;
paraxe.add_double( g20[0], 8) ;
paraxe.add_double( g20[1], 9) ;
paraxe.add_double( g11[0], 10) ;
paraxe.add_double( g11[1], 11) ;
paraxe.add_double( g21[0], 12) ;
paraxe.add_double( g21[1], 13) ;
paraxe.add_double( g22[0], 14) ;
paraxe.add_double( g22[1], 15) ;
paraxe.add_double( bx[0], 16) ;
paraxe.add_double( bx[1], 17) ;
paraxe.add_double( by[0], 18) ;
paraxe.add_double( by[1], 19) ;
paraxe.add_double( bz[0], 20) ;
paraxe.add_double( bz[1], 21) ;
paraxe.add_double( dg00[0], 22) ;
paraxe.add_double( dg00[1], 23) ;
paraxe.add_double( dg10[0], 24) ;
paraxe.add_double( dg10[1], 25) ;
paraxe.add_double( dg20[0], 26) ;
paraxe.add_double( dg20[1], 27) ;
paraxe.add_double( dg11[0], 28) ;
paraxe.add_double( dg11[1], 29) ;
paraxe.add_double( dg21[0], 30) ;
paraxe.add_double( dg21[1], 31) ;
paraxe.add_double( dg22[0], 32) ;
paraxe.add_double( dg22[1], 33) ;
paraxe.add_double( dbx[0], 34) ;
paraxe.add_double( dbx[1], 35) ;
paraxe.add_double( dbz[0], 36) ;
paraxe.add_double( dbz[1], 37) ;
paraxe.add_double( dbymo[0], 38) ;
paraxe.add_double( dbymo[1], 39) ;
paraxe.add_double( d1sn2[0], 40) ;
paraxe.add_double( d1sn2[1], 41) ;
paraxe.add_double( unsn2[0], 42) ;
paraxe.add_double( unsn2[1], 43) ;
paraxe.add_double( omega, 44) ;
int nitmax_axe = 200 ;
int nit_axe ;
double precis_axe = 1.e-13 ;
x_axe = zerosec(fonc_binary_axe, paraxe, 0.9*ori_x1, 0.9*ori_x2,
precis_axe, nitmax_axe, nit_axe) ;
cout << "Binary::orbit : Number of iterations in zerosec for x_axe : "
<< nit_axe << endl ;
}
cout << "Binary::orbit : x_axe [km] : " << x_axe / km << endl ;
//-------------------------------------
// Calcul de la vitesse orbitale
//-------------------------------------
Param parf ;
double ori_x0 = (et[0]->get_mp()).get_ori_x() ;
parf.add_double( ori_x0, 0) ;
parf.add_double( dnulg[0], 1) ;
parf.add_double( g00[0], 2) ;
parf.add_double( g10[0], 3) ;
parf.add_double( g20[0], 4) ;
parf.add_double( g11[0], 5) ;
parf.add_double( g21[0], 6) ;
parf.add_double( g22[0], 7) ;
parf.add_double( bx[0], 8) ;
parf.add_double( by[0], 9) ;
parf.add_double( bz[0], 10) ;
parf.add_double( dg00[0], 11) ;
parf.add_double( dg10[0], 12) ;
parf.add_double( dg20[0], 13) ;
parf.add_double( dg11[0], 14) ;
parf.add_double( dg21[0], 15) ;
parf.add_double( dg22[0], 16) ;
parf.add_double( dbx[0], 17) ;
parf.add_double( dbz[0], 18) ;
parf.add_double( dby[0], 19) ;
parf.add_double( d1sn2[0], 20) ;
parf.add_double( unsn2[0], 21) ;
parf.add_double( x_axe, 22) ;
double omega1 = fact_omeg_min * omega ;
double omega2 = fact_omeg_max * omega ;
cout << "Binary::orbit: omega1, omega2 [rad/s] : "
<< omega1 * f_unit << " " << omega2 * f_unit << endl ;
// Search for the various zeros in the interval [omega1,omega2]
// ------------------------------------------------------------
int nsub = 50 ; // total number of subdivisions of the interval
Tbl* azer = 0x0 ;
Tbl* bzer = 0x0 ;
zero_list(fonc_binary_orbit, parf, omega1, omega2, nsub,
azer, bzer) ;
// Search for the zero closest to the previous value of omega
// ----------------------------------------------------------
double omeg_min, omeg_max ;
int nzer = azer->get_taille() ; // number of zeros found by zero_list
cout << "Binary:orbit : " << nzer <<
" zero(s) found in the interval [omega1, omega2]." << endl ;
cout << "omega, omega1, omega2 : " << omega << " " << omega1
<< " " << omega2 << endl ;
cout << "azer : " << *azer << endl ;
cout << "bzer : " << *bzer << endl ;
if (nzer == 0) {
cout <<
"Binary::orbit: WARNING : no zero detected in the interval"
<< endl << " [" << omega1 * f_unit << ", "
<< omega2 * f_unit << "] rad/s !" << endl ;
omeg_min = omega1 ;
omeg_max = omega2 ;
}
else {
double dist_min = fabs(omega2 - omega1) ;
int i_dist_min = -1 ;
for (int i=0; i<nzer; i++) {
// Distance of previous value of omega from the center of the
// interval [azer(i), bzer(i)]
double dist = fabs( omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
if (dist < dist_min) {
dist_min = dist ;
i_dist_min = i ;
}
}
omeg_min = (*azer)(i_dist_min) ;
omeg_max = (*bzer)(i_dist_min) ;
}
delete azer ; // Tbl allocated by zero_list
delete bzer ; //
cout << "Binary::orbit : interval selected for the search of the zero : "
<< endl << " [" << omeg_min << ", " << omeg_max << "] = ["
<< omeg_min * f_unit << ", " << omeg_max * f_unit << "] rad/s " << endl ;
// Computation of the zero in the selected interval by the secant method
// ---------------------------------------------------------------------
int nitermax = 200 ;
int niter ;
double precis = 1.e-13 ;
omega = zerosec_b(fonc_binary_orbit, parf, omeg_min, omeg_max,
precis, nitermax, niter) ;
cout << "Binary::orbit : Number of iterations in zerosec for omega : "
<< niter << endl ;
cout << "Binary::orbit : omega [rad/s] : "
<< omega * f_unit << endl ;
}
//-------------------------------------------------
// Function used for search of the rotation axis
//-------------------------------------------------
double fonc_binary_axe(double x_rot, const Param& paraxe) {
double ori_x1 = paraxe.get_double(0) ;
double ori_x2 = paraxe.get_double(1) ;
double dnulg_1 = paraxe.get_double(2) ;
double dnulg_2 = paraxe.get_double(3) ;
double g00_1 = paraxe.get_double(4) ;
double g00_2 = paraxe.get_double(5) ;
double g10_1 = paraxe.get_double(6) ;
double g10_2 = paraxe.get_double(7) ;
double g20_1 = paraxe.get_double(8) ;
double g20_2 = paraxe.get_double(9) ;
double g11_1 = paraxe.get_double(10) ;
double g11_2 = paraxe.get_double(11) ;
double g21_1 = paraxe.get_double(12) ;
double g21_2 = paraxe.get_double(13) ;
double g22_1 = paraxe.get_double(14) ;
double g22_2 = paraxe.get_double(15) ;
double bx_1 = paraxe.get_double(16) ;
double bx_2 = paraxe.get_double(17) ;
double by_1 = paraxe.get_double(18) ;
double by_2 = paraxe.get_double(19) ;
double bz_1 = paraxe.get_double(20) ;
double bz_2 = paraxe.get_double(21) ;
double dg00_1 = paraxe.get_double(22) ;
double dg00_2 = paraxe.get_double(23) ;
double dg10_1 = paraxe.get_double(24) ;
double dg10_2 = paraxe.get_double(25) ;
double dg20_1 = paraxe.get_double(26) ;
double dg20_2 = paraxe.get_double(27) ;
double dg11_1 = paraxe.get_double(28) ;
double dg11_2 = paraxe.get_double(29) ;
double dg21_1 = paraxe.get_double(30) ;
double dg21_2 = paraxe.get_double(31) ;
double dg22_1 = paraxe.get_double(32) ;
double dg22_2 = paraxe.get_double(33) ;
double dbx_1 = paraxe.get_double(34) ;
double dbx_2 = paraxe.get_double(35) ;
double dbz_1 = paraxe.get_double(36) ;
double dbz_2 = paraxe.get_double(37) ;
double dbymo_1 = paraxe.get_double(38) ;
double dbymo_2 = paraxe.get_double(39) ;
double d1sn2_1 = paraxe.get_double(40) ;
double d1sn2_2 = paraxe.get_double(41) ;
double unsn2_1 = paraxe.get_double(42) ;
double unsn2_2 = paraxe.get_double(43) ;
double omega = paraxe.get_double(44) ;
double om2_star1 ;
double om2_star2 ;
double x1 = ori_x1 - x_rot ;
double x2 = ori_x2 - x_rot ;
double bymxo_1 = by_1-x1*omega ;
double bymxo_2 = by_2-x2*omega ;
double beta1 = g00_1*bx_1*bx_1 + 2*g10_1*bx_1*bymxo_1 + 2*g20_1*bx_1*bz_1 ;
double beta2 = g11_1*bymxo_1*bymxo_1 + 2*g21_1*bz_1*bymxo_1
+ g22_1*bz_1*bz_1 ;
double beta_1 = beta1 + beta2 ;
double delta1 = dg00_1*bx_1*bx_1 + 2*g00_1*dbx_1*bx_1
+ 2*dg10_1*bx_1*bymxo_1 ;
double delta2 = 2*g10_1*bymxo_1*dbx_1 + 2*g10_1*bx_1*dbymo_1
+ 2*dg20_1*bx_1*bz_1 ;
double delta3 = 2*g20_1*bx_1*dbz_1 +2*g20_1*bz_1*dbx_1 + dg11_1*bymxo_1*bymxo_1 ;
double delta4 = 2*g11_1*bymxo_1*dbymo_1 + 2*dg21_1*bz_1*bymxo_1;
double delta5 = 2*g21_1*bymxo_1*dbz_1 +2*g21_1*bz_1*dbymo_1
+ dg22_1*bz_1*bz_1 + 2*g22_1*bz_1*dbz_1 ;
double delta_1 = delta1 + delta2 + delta3 + delta4 + delta5 ;
// Computation of omega for star 1
//---------------------------------
om2_star1 = dnulg_1 / (beta_1/(omega*omega)*(dnulg_1*unsn2_1 + d1sn2_1/2.)
+ unsn2_1*delta_1/(omega*omega)/2.) ;
double beta3 = g00_2*bx_2*bx_2 + 2*g10_2*bx_2*bymxo_2 + 2*g20_2*bx_2*bz_2 ;
double beta4 = g11_2*bymxo_2*bymxo_2 + 2*g21_2*bz_2*bymxo_2
+ g22_2*bz_2*bz_2 ;
double beta_2 = beta3 + beta4 ;
double delta6 = dg00_2*bx_2*bx_2 + 2*g00_2*dbx_2*bx_2
+ 2*dg10_2*bx_2*bymxo_2 ;
double delta7 = 2*g10_2*bymxo_2*dbx_2 + 2*g10_2*bx_2*dbymo_2
+ 2*dg20_2*bx_2*bz_2 ;
double delta8 = 2*g20_2*bx_2*dbz_2 +2*g20_2*bz_2*dbx_2
+ dg11_2*bymxo_2*bymxo_2 ;
double delta9 = 2*g11_2*bymxo_2*dbymo_2 + 2*dg21_2*bz_2*bymxo_2;
double delta10 = 2*g21_2*bymxo_2*dbz_2 +2*g21_2*bz_2*dbymo_2
+ dg22_2*bz_2*bz_2 + 2*g22_2*bz_2*dbz_2 ;
double delta_2 = delta6 + delta7 + delta8 + delta9 + delta10 ;
// Computation of omega for star 2
//---------------------------------
om2_star2 = dnulg_2 / (beta_2/(omega*omega)*(dnulg_2*unsn2_2 + d1sn2_2/2.)
+ unsn2_2*delta_2/(omega*omega)/2.) ;
;
return om2_star1 - om2_star2 ;
}
//----------------------------------------------------------------------------
// Fonction utilisee pour la recherche de omega par la methode de la secante
//----------------------------------------------------------------------------
double fonc_binary_orbit(double om, const Param& parf) {
double xc = parf.get_double(0) ;
double dnulg = parf.get_double(1) ;
double g00 = parf.get_double(2) ;
double g10 = parf.get_double(3) ;
double g20 = parf.get_double(4) ;
double g11 = parf.get_double(5) ;
double g21 = parf.get_double(6) ;
double g22 = parf.get_double(7) ;
double bx = parf.get_double(8) ;
double by = parf.get_double(9) ;
double bz = parf.get_double(10) ;
double dg00 = parf.get_double(11) ;
double dg10 = parf.get_double(12) ;
double dg20 = parf.get_double(13) ;
double dg11 = parf.get_double(14) ;
double dg21 = parf.get_double(15) ;
double dg22 = parf.get_double(16) ;
double dbx = parf.get_double(17) ;
double dbz = parf.get_double(18) ;
double dby = parf.get_double(19) ;
double d1sn2 = parf.get_double(20) ;
double unsn2 = parf.get_double(21) ;
double x_axe = parf.get_double(22) ;
double dbymo = dby - om ;
double xx = xc - x_axe ;
double bymxo = by-xx*om ;
// bymxo = - bymxo ;
//dbymo = - dbymo ;
double beta1 = g00*bx*bx + 2*g10*bx*bymxo + 2*g20*bx*bz ;
double beta2 = g11*bymxo*bymxo + 2*g21*bz*bymxo+ g22*bz*bz ;
double beta = beta1 + beta2 ;
double alpha = 1 - unsn2*beta ;
double delta1 = dg00*bx*bx + 2*g00*dbx*bx + 2*dg10*bx*bymxo ;
double delta2 = 2*g10*bymxo*dbx + 2*g10*bx*dbymo + 2*dg20*bx*bz ;
double delta3 = 2*g20*bx*dbz +2*g20*bz*dbx + dg11*bymxo*bymxo ;
double delta4 = 2*g11*bymxo*dbymo + 2*dg21*bz*bymxo;
double delta5 = 2*g21*bymxo*dbz +2*g21*bz*dbymo + dg22*bz*bz
+ 2*g22*bz*dbz ;
double delta = delta1 + delta2 + delta3 + delta4 + delta5 ;
// Difference entre les 2 termes de l'eq.(95) de Gourgoulhon et.al (2001)
//centre de l'etoile
//-----------------------------------------------------------------------
double diff = dnulg + (1/(2.*alpha))*(-d1sn2*beta - unsn2*delta) ;
/*
double bpb = om*om*xx*xx - 2*om*by*xx + by*by ;
double dphi_cent = (g00*unsn2*( om*(by+xx*dby) - om*om*xx - by*dby )
- 0.5*bpb*(dg00 + g00*d1sn2/unsn2)*unsn2 )
/ (1 - g00*unsn2*bpb) ;
double diff = dnulg + dphi_cent ;
cout.precision(8) ;
cout << "bpb = " << bpb << endl ;
cout << "om = " << om << endl ;
cout << "by = " << by << endl ;
cout << "xx = " << xx << endl ;
cout << "dby = " << dby << endl ;
cout << "part11 = " << g00*unsn2 << endl ;
cout << "part12 = " << om*(by+xx*dby) << endl ;
cout << "part13 = " <<- om*om*xx << endl ;
cout << "part14 = " << - by*dby << endl ;
cout << "part1 = " << g00*unsn2*( om*(by+xx*dby) - om*om*xx - by*dby ) << endl ;
cout << "part2 = " << - 0.5*bpb*(dg00 + g00*d1sn2/unsn2)*unsn2 << endl ;
cout << "part3 = " << (1 - g00*unsn2*bpb) << endl ;
cout << "dnulg = " << dnulg << endl ;
cout << "dphi_cent = " << dphi_cent << endl ;
cout << "diff = " << diff << endl ;
cout << endl ;
*/
return diff ;
}
}
Lorene::Binary::et
Star_bin * et[2]
Array of the two stars (to perform loops on the stars): et[0] contains the address of star1 and et[1]...
Definition binary.h:89
Lorene::Binary::omega
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition binary.h:94
Lorene::Binary::orbit
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.
Definition binary_orbite.C:87
Lorene::Binary::x_axe
double x_axe
Absolute X coordinate of the rotation axis.
Definition binary.h:98
Lorene::Metric
Metric for tensor calculation.
Definition metric.h:90
Lorene::Metric::cov
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition metric.C:280
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::Param::add_double
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition param.C:315
Lorene::Param::get_double
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition param.C:361
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Lorene::Scalar::val_grid_point
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition scalar.h:637
Lorene::Scalar::dsdx
const Scalar & dsdx() const
Returns of *this , where .
Definition scalar_deriv.C:266
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::get_taille
int get_taille() const
Gives the total size (ie dim.taille).
Definition tbl.h:397
Lorene::Tensor
Tensor handling.
Definition tensor.h:288
Lorene::Vector
Tensor field of valence 1.
Definition vector.h:188
Lorene::Vector::change_triad
virtual void change_triad(const Base_vect &)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition vector_change_triad.C:75
Lorene::norme
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition cmp_math.C:481
Lorene::zero_list
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.
Definition zero_list.C:57
Lorene::zerosec_b
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.
Definition zerosec_borne.C:68
Lorene::zerosec
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.
Definition zerosec.C:89
Lorene::Tensor::change_triad
virtual void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition tensor_change_triad.C:77
Lorene
Lorene prototypes.
Definition app_hor.h:64
Unites
Standard units of space, time and mass.