30 char compobj_QI_global_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Compobj/compobj_QI_global.C,v 1.11 2014/10/13 08:52:49 j_novak Exp $" ;
79 #include "utilitaires.h"
82 double funct_compobj_QI_isco(
double,
const Param&) ;
83 double funct_compobj_QI_rmb(
double,
const Param&) ;
94 cerr <<
"Compobj_QI::angu_mom() : not implemented yet !" << endl ;
146 Scalar ucor_plus = (r2 * bsn * dnphi +
147 sqrt ( r2 * r2 * bsn *bsn * dnphi * dnphi +
148 4 * r2 * bdlog * ndlog + 4 * r * ndlog) ) /
149 2 / (1 + r * bdlog ) ;
151 Scalar factor_u2 = r2 * (2 * ddbbb /
bbb - 2 * bdlog * bdlog +
152 4 * bdlog * ndlog ) +
153 2 * r2 * r2 * bsn * bsn * dnphi * dnphi +
154 4 * r * ( ndlog - bdlog ) - 6 ;
156 Scalar factor_u1 = 2 * r * r2 * bsn * ( 2 * ( ndlog - bdlog ) *
159 Scalar factor_u0 = - r2 * ( 2 * ddnnn /
nn - 2 * ndlog * ndlog +
160 4 * bdlog * ndlog ) ;
163 Scalar orbit = factor_u2 * ucor_plus * ucor_plus +
164 factor_u1 * ucor_plus + factor_u0 ;
170 double r_ms, theta_ms, phi_ms, xi_ms,
171 xi_min = -1, xi_max = 1;
174 bool find_status = false ;
176 const Valeur& vorbit = orbit.get_spectral_va() ;
180 theta_ms = M_PI / 2. ;
183 for(l = nzm1-1; l >= lmin; l--) {
189 double xi_f = xi_min;
191 double resloc = vorbit.
val_point(l, xi_min, theta_ms, phi_ms) ;
193 for (
int iloc=0; iloc<200; iloc++) {
195 resloc_old = resloc ;
196 resloc = vorbit.
val_point(l, xi_f, theta_ms, phi_ms) ;
197 if ( resloc * resloc_old <
double(0) ) {
198 xi_min = xi_f - 0.01 ;
206 if (find_status) break ;
218 double precis_ms = 1.e-12 ;
220 int nitermax_ms = 100 ;
223 xi_ms =
zerosec(funct_compobj_QI_isco, par_ms, xi_min, xi_max,
224 precis_ms, nitermax_ms, niter) ;
226 *ost <<
"ISCO search: " << endl ;
227 *ost <<
" Domain number: " << l_ms << endl ;
228 *ost <<
" xi_min, xi_max : " << xi_min <<
" , " << xi_max << endl ;
229 *ost <<
" number of iterations used in zerosec: " << niter << endl ;
230 *ost <<
" zero found at xi = " << xi_ms << endl ;
233 r_ms =
mp.
val_r(l_ms, xi_ms, theta_ms, phi_ms) ;
254 double ucor_msplus = (ucor_plus.
get_spectral_va()).val_point(l_ms, xi_ms, theta_ms,
256 double nobrs = (bsn.get_spectral_va()).val_point(l_ms, xi_ms, theta_ms, phi_ms) ;
259 p_f_isco =
new double ( ( ucor_msplus / nobrs / r_ms + nphirs ) /
260 (
double(2) * M_PI) ) ;
269 p_espec_isco=
new double (( 1./nobrs / r_ms / ucor_msplus + nphirs) *
270 (ucor_msplus/
sqrt(1.-ucor_msplus*ucor_msplus)*
372 Scalar C1 = (A1 * B3 - A3 * B1) ;
373 Scalar C2 = (A2 * B3 - A3 * B2) ;
374 Scalar C3 = (A4 * B3 - A3 * B4) ;
377 Scalar D1 = B4 * C1 - B1 * C3 ;
378 Scalar D2 = B4 * C2 - B2 * C3 ;
379 Scalar D3 = B3 * B3 * C1 * C3 ;
380 Scalar D4 = B3 * B3 * C2 * C3 ;
399 Scalar bound_orbit = -(2.*D1*D2 + D3) -
sqrt((2.*D1*D2 + D3) * (2.*D1*D2 + D3) - 4.*D1*D1 *
400 (D2*D2 + D4)) - 2.*(D2*D2 + D4) ;
413 double zeros[2][noz] ;
417 double rmb, theta_mb, phi_mb, xi_mb;
418 double xi_min = -1, xi_max = 1 ;
428 theta_mb = M_PI / 2. ;
432 for(l = lmin; l <= nzm1; l++) {
437 double xi_f = xi_min;
439 double resloc = vorbit.
val_point(l, xi_f, theta_mb, phi_mb) ;
441 while(xi_f <= xi_max) {
445 resloc_old = resloc ;
446 resloc = vorbit.
val_point(l, xi_f, theta_mb, phi_mb) ;
448 if ( resloc * resloc_old <
double(0) ) {
462 int number_of_zeros = i ;
464 cout <<
"number of zeros: " << number_of_zeros << endl ;
467 double precis_mb = 1.e-9 ;
470 int nitermax_mb = 100 ;
473 for(
int i = 0; i < number_of_zeros; i++) {
477 int l_mb = int(zeros[1][i]) ;
478 xi_max = zeros[0][i] ;
479 xi_min = xi_max - dx ;
486 xi_mb =
zerosec(funct_compobj_QI_rmb, par_mb, xi_min, xi_max, precis_mb, nitermax_mb, niter) ;
488 *ost <<
"RMB search: " << endl ;
489 *ost <<
" Domain number: " << l_mb << endl ;
490 *ost <<
" xi_min, xi_max : " << xi_min <<
" , " << xi_max << endl ;
491 *ost <<
" number of iterations used in zerosec: " << niter << endl ;
492 *ost <<
" zero found at xi = " << xi_mb << endl ;
495 if (niter < nitermax_mb) {
496 double zero_mb =
mp.
val_r(l_mb, xi_mb, theta_mb, phi_mb) ;
499 if (zero_mb < (1 + r_ms)/2){
506 p_r_mb =
new double (rmb) ;
528 double funct_compobj_QI_isco(
double xi,
const Param& par){
536 double theta = M_PI / 2. ;
538 return vorbit.
val_point(l_ms, xi, theta, phi) ;
548 double funct_compobj_QI_rmb(
double zeros,
const Param& par){
551 int l_mb = par.get_int() ;
552 const Scalar& orbit = par.get_scalar() ;
553 const Valeur& vorbit = orbit.get_spectral_va() ;
556 double theta = M_PI / 2. ;
558 return vorbit.
val_point(l_mb, zeros, theta, phi) ;