77 #include "utilitaires.h" 80 double fonct_etoile_rot_isco(
double,
const Param&) ;
125 Cmp ucor_plus = (r2 * bsn * dnphi +
126 sqrt ( r2 * r2 * bsn *bsn * dnphi * dnphi +
127 4 * r2 * bdlog * ndlog + 4 * r * ndlog) ) /
128 2 / (1 + r * bdlog ) ;
130 Cmp factor_u2 = r2 * (2 * ddbbb /
bbb() - 2 * bdlog * bdlog +
131 4 * bdlog * ndlog ) +
132 2 * r2 * r2 * bsn * bsn * dnphi * dnphi +
133 4 * r * ( ndlog - bdlog ) - 6 ;
135 Cmp factor_u1 = 2 * r * r2 * bsn * ( 2 * ( ndlog - bdlog ) *
138 Cmp factor_u0 = - r2 * ( 2 * ddnnn /
nnn() - 2 * ndlog * ndlog +
139 4 * bdlog * ndlog ) ;
142 Cmp orbit = factor_u2 * ucor_plus * ucor_plus +
143 factor_u1 * ucor_plus + factor_u0 ;
149 double r_ms, theta_ms, phi_ms, xi_ms,
150 xi_min = -1, xi_max = 1;
152 bool find_status = false ;
154 const Valeur& vorbit = orbit.va ;
156 for(l =
nzet; l <= nzm1; l++) {
160 theta_ms = M_PI / 2. ;
167 double xi_f = xi_min;
169 double resloc = vorbit.
val_point(l, xi_min, theta_ms, phi_ms) ;
171 for (
int iloc=0; iloc<200; iloc++) {
173 resloc_old = resloc ;
174 resloc = vorbit.
val_point(l, xi_f, theta_ms, phi_ms) ;
175 if ( resloc * resloc_old <
double(0) ) {
176 xi_min = xi_f - 0.01 ;
193 double precis_ms = 1.e-13 ;
194 int nitermax_ms = 200 ;
197 xi_ms =
zerosec(fonct_etoile_rot_isco, par_ms, xi_min, xi_max,
198 precis_ms, nitermax_ms, niter,
false) ;
200 if (niter > nitermax_ms) {
201 cerr <<
"Etoile_rot::r_isco : " << endl ;
202 cerr <<
"result may be wrong ... " << endl ;
203 cerr <<
"Continuing." << endl ;
208 " number of iterations used in zerosec to locate the ISCO : " 210 *ost <<
" zero found at xi = " << xi_ms << endl ;
213 r_ms =
mp.
val_r(l_ms, xi_ms, theta_ms, phi_ms) ;
224 p_r_isco =
new double ((
bbb().va).val_point(l_ms, xi_ms, theta_ms, phi_ms)
231 double ucor_msplus = (ucor_plus.
va).val_point(l_ms, xi_ms, theta_ms,
233 double nobrs = (bsn.va).val_point(l_ms, xi_ms, theta_ms, phi_ms) ;
234 double nphirs = (
nphi().va).val_point(l_ms, xi_ms, theta_ms, phi_ms) ;
236 p_f_isco =
new double ( ( ucor_msplus / nobrs / r_ms + nphirs ) /
237 (
double(2) * M_PI) ) ;
242 ((
bbb().va).val_point(l_ms, xi_ms, theta_ms, phi_ms)) * r_ms );
246 p_espec_isco=
new double (( 1./nobrs / r_ms / ucor_msplus + nphirs) *
247 (ucor_msplus/
sqrt(1.-ucor_msplus*ucor_msplus)*
248 ((
bbb().va).val_point(l_ms, xi_ms, theta_ms, phi_ms)) * r_ms ));
253 double ucor_eqplus = (ucor_plus.
va).val_point(l_ms, -1, theta_ms, phi_ms) ;
254 double nobeq = (bsn.va).val_point(l_ms, -1, theta_ms, phi_ms) ;
255 double nphieq = (
nphi().va).val_point(l_ms, -1, theta_ms, phi_ms) ;
257 p_f_eq =
new double ( ( ucor_eqplus / nobeq /
ray_eq() + nphieq ) /
258 (
double(2) * M_PI) ) ;
341 double fonct_etoile_rot_isco(
double xi,
const Param& par){
349 double theta = M_PI / 2. ;
351 return vorbit.
val_point(l_ms, xi, theta, phi) ;
const Cmp & dsdr() const
Returns of *this .
const Cmp & get_cmp(int position=0) const
Returns the reference of a Cmp stored in the list.
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
virtual double lspec_isco() const
Angular momentum of a particle on the ISCO.
double * p_r_isco
Circumferential radius of the ISCO.
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Cmp sqrt(const Cmp &)
Square root.
void annule(int l)
Sets the Cmp to zero in a given domain.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Tenseur nnn
Total lapse function.
Tenseur nphi
Metric coefficient .
double ray_eq() const
Coordinate radius at , [r_unit].
Values and coefficients of a (real-value) function.
virtual double r_isco(ostream *ost=0x0) const
Circumferential radius of the innermost stable circular orbit (ISCO).
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.
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
void add_cmp(const Cmp &ti, int position=0)
Adds the address of a new Cmp to the list.
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 ...
Map & mp
Mapping associated with the star.
int get_nzone() const
Returns the number of domains.
virtual double espec_isco() const
Energy of a particle on the ISCO.
double * p_f_isco
Orbital frequency of the ISCO.
virtual double f_eq() const
Orbital frequency at the equator.
Tenseur bbb
Metric factor B.
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
int nzet
Number of domains of *mp occupied by the star.
double * p_lspec_isco
Specific angular momentum of a particle on the ISCO.
virtual double f_isco() const
Orbital frequency at the innermost stable circular orbit (ISCO).
Valeur va
The numerical value of the Cmp.
double * p_espec_isco
Specific energy of a particle on the ISCO.
Coord r
r coordinate centered on the grid
double * p_f_eq
Orbital frequency at the equator.