80 cout <<
"Etoile_rot::f_eccentric not ready yet !" << endl ;
117 Cmp ucor_plus = (r2 * bsn * dnphi +
118 sqrt ( r2 * r2 * bsn *bsn * dnphi * dnphi +
119 4 * r2 * bdlog * ndlog + 4 * r * ndlog) ) /
120 2 / (1 + r * bdlog ) ;
122 Cmp factor_u2 = r2 * (2 * ddbbb /
bbb() - 2 * bdlog * bdlog +
123 4 * bdlog * ndlog ) +
124 2 * r2 * r2 * bsn * bsn * dnphi * dnphi +
125 4 * r * ( ndlog - bdlog ) - 6 ;
127 Cmp factor_u1 = 2 * r * r2 * bsn * ( 2 * ( ndlog - bdlog ) *
130 Cmp factor_u0 = - r2 * ( 2 * ddnnn /
nnn() - 2 * ndlog * ndlog +
131 4 * bdlog * ndlog ) ;
134 Cmp orbit = factor_u2 * ucor_plus * ucor_plus +
135 factor_u1 * ucor_plus + factor_u0 ;
150 double theta_ms = M_PI / 2. ;
153 double xi_min = -1. ;
157 double xi_f = xi_min;
159 orbit.std_base_scal() ;
160 const Valeur& vorbit = orbit.va ;
162 double resloc = vorbit.
val_point(l_ms, xi_min, theta_ms, phi_ms) ;
164 for (
int iloc=0; iloc<200; iloc++) {
166 resloc_old = resloc ;
167 resloc = vorbit.
val_point(l_ms, xi_f, theta_ms, phi_ms) ;
168 if ((resloc * resloc_old) < double(0) ) {
169 xi_min = xi_f - 0.01 ;
178 "Etoile_rot::isco : preliminary location of zero of MS function :"
180 *ost <<
" xi_min = " << xi_min <<
" f(xi_min) = " <<
181 vorbit.
val_point(l_ms, xi_min, theta_ms, phi_ms) << endl ;
182 *ost <<
" xi_max = " << xi_max <<
" f(xi_max) = " <<
183 vorbit.
val_point(l_ms, xi_max, theta_ms, phi_ms) << endl ;
189 if ( vorbit.
val_point(l_ms, xi_min, theta_ms, phi_ms) *
190 vorbit.
val_point(l_ms, xi_max, theta_ms, phi_ms) < double(0) ) {
199 " number of iterations used in zerosec to locate the ISCO : "
201 *ost <<
" zero found at xi = " << xi_ms << endl ;
204 r_ms =
mp.
val_r(l_ms, xi_ms, theta_ms, phi_ms) ;
213 (
bbb().va).val_point(l_ms, xi_ms, theta_ms, phi_ms) * r_ms
220 double ucor_msplus = (ucor_plus.
va).val_point(l_ms, xi_ms, theta_ms,
222 double nobrs = (bsn.
va).val_point(l_ms, xi_ms, theta_ms, phi_ms) ;
223 double nphirs = (
nphi().va).val_point(l_ms, xi_ms, theta_ms, phi_ms) ;
225 p_f_isco =
new double ( ( ucor_msplus / nobrs / r_ms + nphirs ) /
226 (
double(2) * M_PI) ) ;
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.
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 ...