122 Cmp ucor_plus = (r2 * bsn * dnphi +
123 sqrt ( r2 * r2 * bsn *bsn * dnphi * dnphi +
124 4 * r2 * bdlog * ndlog + 4 * r * ndlog) ) /
125 2 / (1 + r * bdlog ) ;
127 Cmp factor_u2 = r2 * (2 * ddbbb /
bbb() - 2 * bdlog * bdlog +
128 4 * bdlog * ndlog ) +
129 2 * r2 * r2 * bsn * bsn * dnphi * dnphi +
130 4 * r * ( ndlog - bdlog ) - 6 ;
132 Cmp factor_u1 = 2 * r * r2 * bsn * ( 2 * ( ndlog - bdlog ) *
135 Cmp factor_u0 = - r2 * ( 2 * ddnnn /
nnn() - 2 * ndlog * ndlog +
136 4 * bdlog * ndlog ) ;
139 Cmp orbit = factor_u2 * ucor_plus * ucor_plus +
140 factor_u1 * ucor_plus + factor_u0 ;
146 double r_ms, theta_ms, phi_ms, xi_ms,
147 xi_min = -1, xi_max = 1;
149 bool find_status = false ;
151 const Valeur& vorbit = orbit.va ;
153 for(l =
nzet; l <= nzm1; l++) {
157 theta_ms = M_PI / 2. ;
164 double xi_f = xi_min;
166 double resloc = vorbit.
val_point(l, xi_min, theta_ms, phi_ms) ;
168 for (
int iloc=0; iloc<200; iloc++) {
170 resloc_old = resloc ;
171 resloc = vorbit.
val_point(l, xi_f, theta_ms, phi_ms) ;
172 if ( resloc * resloc_old <
double(0) ) {
173 xi_min = xi_f - 0.01 ;
190 double precis_ms = 1.e-13 ;
191 int nitermax_ms = 200 ;
194 xi_ms =
zerosec(fonct_etoile_rot_isco, par_ms, xi_min, xi_max,
195 precis_ms, nitermax_ms, niter,
false) ;
197 if (niter > nitermax_ms) {
198 cerr <<
"Etoile_rot::r_isco : " << endl ;
199 cerr <<
"result may be wrong ... " << endl ;
200 cerr <<
"Continuing." << endl ;
205 " number of iterations used in zerosec to locate the ISCO : "
207 *ost <<
" zero found at xi = " << xi_ms << endl ;
210 r_ms =
mp.
val_r(l_ms, xi_ms, theta_ms, phi_ms) ;
221 p_r_isco =
new double ((
bbb().va).val_point(l_ms, xi_ms, theta_ms, phi_ms)
228 double ucor_msplus = (ucor_plus.
va).val_point(l_ms, xi_ms, theta_ms,
230 double nobrs = (bsn.
va).val_point(l_ms, xi_ms, theta_ms, phi_ms) ;
231 double nphirs = (
nphi().va).val_point(l_ms, xi_ms, theta_ms, phi_ms) ;
233 p_f_isco =
new double ( ( ucor_msplus / nobrs / r_ms + nphirs ) /
234 (
double(2) * M_PI) ) ;
239 ((
bbb().va).val_point(l_ms, xi_ms, theta_ms, phi_ms)) * r_ms );
243 p_espec_isco=
new double (( 1./nobrs / r_ms / ucor_msplus + nphirs) *
244 (ucor_msplus/
sqrt(1.-ucor_msplus*ucor_msplus)*
245 ((
bbb().va).val_point(l_ms, xi_ms, theta_ms, phi_ms)) * r_ms ));
250 double ucor_eqplus = (ucor_plus.
va).val_point(l_ms, -1, theta_ms, phi_ms) ;
251 double nobeq = (bsn.
va).val_point(l_ms, -1, theta_ms, phi_ms) ;
252 double nphieq = (
nphi().va).val_point(l_ms, -1, theta_ms, phi_ms) ;
254 p_f_eq =
new double ( ( ucor_eqplus / nobeq /
ray_eq() + nphieq ) /
255 (
double(2) * M_PI) ) ;
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.