et_rot_f_eccentric.C
/*
* Method of class Etoile_rot to compute eccentric orbits
*
* (see file etoile.h for documentation).
*
*/
/*
* Copyright (c) 2001 Dorota Gondek-Rosinska
* Copyright (c) 2001 Eric Gourgoulhon
*
* 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 et_rot_f_eccentric_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_f_eccentric.C,v 1.7 2014/10/13 08:52:57 j_novak Exp $" ;
/*
* $Id: et_rot_f_eccentric.C,v 1.7 2014/10/13 08:52:57 j_novak Exp $
* $Log: et_rot_f_eccentric.C,v $
* Revision 1.7 2014/10/13 08:52:57 j_novak
* Lorene classes and functions now belong to the namespace Lorene.
*
* Revision 1.6 2014/10/06 15:13:09 j_novak
* Modified #include directives to use c++ syntax.
*
* Revision 1.5 2003/12/19 16:31:52 j_novak
* Still warnings...
*
* Revision 1.4 2003/12/19 16:21:42 j_novak
* Shadow hunt
*
* Revision 1.3 2003/12/05 14:50:26 j_novak
* To suppress some warnings...
*
* Revision 1.2 2003/10/03 15:58:47 j_novak
* Cleaning of some headers
*
* Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
* LORENE
*
* Revision 2.0 2001/02/08 15:13:24 eric
* *** empty log message ***
*
*
* $Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_f_eccentric.C,v 1.7 2014/10/13 08:52:57 j_novak Exp $
*
*/
// Headers C
#include <cmath>
// Headers Lorene
#include "etoile.h"
#include "param.h"
//=============================================================================
namespace Lorene {
// r_isco()
//=============================================================================
double Etoile_rot::f_eccentric(double, double, ostream* ost) const {
cout << "Etoile_rot::f_eccentric not ready yet !" << endl ;
abort() ;
// First and second derivatives of metric functions
// ------------------------------------------------
int nzm1 = mp.get_mg()->get_nzone() - 1 ;
Cmp dnphi = nphi().dsdr() ;
dnphi.annule(nzm1) ;
Cmp ddnphi = dnphi.dsdr() ; // d^2/dr^2 N^phi
Cmp tmp = nnn().dsdr() ;
tmp.annule(nzm1) ;
Cmp ddnnn = tmp.dsdr() ; // d^2/dr^2 N
tmp = bbb().dsdr() ;
tmp.annule(nzm1) ;
Cmp ddbbb = tmp.dsdr() ; // d^2/dr^2 B
// Constructing the velocity of a particle corotating with the star
// ----------------------------------------------------------------
Cmp bdlog = bbb().dsdr() / bbb() ;
Cmp ndlog = nnn().dsdr() / nnn() ;
Cmp bsn = bbb() / nnn() ;
Cmp r(mp) ;
r = mp.r ;
Cmp r2= r*r ;
bdlog.annule(nzm1) ;
ndlog.annule(nzm1) ;
bsn.annule(nzm1) ;
r2.annule(nzm1) ;
// ucor_plus - the velocity of corotating particle on the circular orbit
Cmp ucor_plus = (r2 * bsn * dnphi +
sqrt ( r2 * r2 * bsn *bsn * dnphi * dnphi +
4 * r2 * bdlog * ndlog + 4 * r * ndlog) ) /
2 / (1 + r * bdlog ) ;
Cmp factor_u2 = r2 * (2 * ddbbb / bbb() - 2 * bdlog * bdlog +
4 * bdlog * ndlog ) +
2 * r2 * r2 * bsn * bsn * dnphi * dnphi +
4 * r * ( ndlog - bdlog ) - 6 ;
Cmp factor_u1 = 2 * r * r2 * bsn * ( 2 * ( ndlog - bdlog ) *
dnphi - ddnphi ) ;
Cmp factor_u0 = - r2 * ( 2 * ddnnn / nnn() - 2 * ndlog * ndlog +
4 * bdlog * ndlog ) ;
// Scalar field the zero of which will give the marginally stable orbit
Cmp orbit = factor_u2 * ucor_plus * ucor_plus +
factor_u1 * ucor_plus + factor_u0 ;
// Search for the zero
// -------------------
int l_ms = nzet ; // index of the domain containing the MS orbit
Param par_ms ;
par_ms.add_int(l_ms) ;
par_ms.add_cmp(orbit) ;
// Preliminary location of the zero
// of the orbit function with an error = 0.01
// The orbit closest to the star
double theta_ms = M_PI / 2. ; // pi/2
double phi_ms = 0. ;
double xi_min = -1. ;
double xi_max = 1. ;
double resloc_old ;
double xi_f = xi_min;
orbit.std_base_scal() ;
const Valeur& vorbit = orbit.va ;
double resloc = vorbit.val_point(l_ms, xi_min, theta_ms, phi_ms) ;
for (int iloc=0; iloc<200; iloc++) {
xi_f = xi_f + 0.01 ;
resloc_old = resloc ;
resloc = vorbit.val_point(l_ms, xi_f, theta_ms, phi_ms) ;
if ((resloc * resloc_old) < double(0) ) {
xi_min = xi_f - 0.01 ;
xi_max = xi_f ;
break ;
}
}
if (ost != 0x0) {
*ost <<
"Etoile_rot::isco : preliminary location of zero of MS function :"
<< endl ;
*ost << " xi_min = " << xi_min << " f(xi_min) = " <<
vorbit.val_point(l_ms, xi_min, theta_ms, phi_ms) << endl ;
*ost << " xi_max = " << xi_max << " f(xi_max) = " <<
vorbit.val_point(l_ms, xi_max, theta_ms, phi_ms) << endl ;
}
double xi_ms = 0 ;
double r_ms = 0 ;
if ( vorbit.val_point(l_ms, xi_min, theta_ms, phi_ms) *
vorbit.val_point(l_ms, xi_max, theta_ms, phi_ms) < double(0) ) {
//## double precis_ms = 1.e-12 ; // precision in the determination of xi_ms
//## int nitermax_ms = 100 ; // max number of iterations
int niter = 0 ;
if (ost != 0x0) {
* ost <<
" number of iterations used in zerosec to locate the ISCO : "
<< niter << endl ;
*ost << " zero found at xi = " << xi_ms << endl ;
}
r_ms = mp.val_r(l_ms, xi_ms, theta_ms, phi_ms) ;
}
else {
xi_ms = -1 ;
r_ms = ray_eq() ;
}
p_r_isco = new double (
(bbb().va).val_point(l_ms, xi_ms, theta_ms, phi_ms) * r_ms
) ;
// Determination of the frequency at the marginally stable orbit
// -------------------------------------------------------------
ucor_plus.std_base_scal() ;
double ucor_msplus = (ucor_plus.va).val_point(l_ms, xi_ms, theta_ms,
phi_ms) ;
double nobrs = (bsn.va).val_point(l_ms, xi_ms, theta_ms, phi_ms) ;
double nphirs = (nphi().va).val_point(l_ms, xi_ms, theta_ms, phi_ms) ;
p_f_isco = new double ( ( ucor_msplus / nobrs / r_ms + nphirs ) /
(double(2) * M_PI) ) ;
return 0 ;
}
}
Lorene::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Lorene::Cmp::va
Valeur va
The numerical value of the Cmp.
Definition cmp.h:464
Lorene::Cmp::std_base_scal
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition cmp.C:644
Lorene::Cmp::annule
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:348
Lorene::Cmp::dsdr
const Cmp & dsdr() const
Returns of *this .
Definition cmp_deriv.C:84
Lorene::Etoile_rot::f_eccentric
virtual double f_eccentric(double ecc, double periast, ostream *ost=0x0) const
Computation of frequency of eccentric orbits.
Definition et_rot_f_eccentric.C:78
Lorene::Etoile_rot::nphi
Tenseur nphi
Metric coefficient .
Definition etoile.h:1510
Lorene::Etoile_rot::bbb
Tenseur bbb
Metric factor B.
Definition etoile.h:1504
Lorene::Etoile_rot::p_f_isco
double * p_f_isco
Orbital frequency of the ISCO.
Definition etoile.h:1645
Lorene::Etoile_rot::p_r_isco
double * p_r_isco
Circumferential radius of the ISCO.
Definition etoile.h:1644
Lorene::Etoile::nzet
int nzet
Number of domains of *mp occupied by the star.
Definition etoile.h:432
Lorene::Etoile::ray_eq
double ray_eq() const
Coordinate radius at , [r_unit].
Definition etoile_global.C:120
Lorene::Etoile::nnn
Tenseur nnn
Total lapse function.
Definition etoile.h:509
Lorene::Etoile::mp
Map & mp
Mapping associated with the star.
Definition etoile.h:429
Lorene::Param
Parameter storage.
Definition param.h:125
Lorene::Param::add_cmp
void add_cmp(const Cmp &ti, int position=0)
Adds the address of a new Cmp to the list.
Definition param.C:935
Lorene::Param::add_int
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:246
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Lorene::Valeur::val_point
double val_point(int l, double x, double theta, double phi) const
Computes the value of the field represented by *this at an arbitrary point, by means of the spectral ...
Definition valeur.C:882
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Lorene
Lorene prototypes.
Definition app_hor.h:64
Lorene::r
Coord r
r coordinate centered on the grid
Definition map.h:718