30char 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"
82double funct_compobj_QI_isco(
double,
const Param&) ;
83double 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) ;
528double 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) ;
548double 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) ;
double * p_lspec_isco
Specific angular momentum of a particle at the ISCO.
virtual double f_isco(int lmin) const
Orbital frequency at the innermost stable circular orbit (ISCO).
double * p_r_mb
Coordinate r of the marginally bound orbit.
double * p_r_isco
Coordinate r of the ISCO.
virtual double lspec_isco(int lmin) const
Angular momentum of a particle at the ISCO.
virtual double angu_mom() const
Angular momentum.
virtual double r_isco(int lmin, ostream *ost=0x0) const
Coordinate r of the innermost stable circular orbit (ISCO).
Scalar nphi
Metric coefficient .
virtual double espec_isco(int lmin) const
Energy of a particle at the ISCO.
Scalar bbb
Metric factor B.
double * p_espec_isco
Specific energy of a particle at the ISCO.
double * p_angu_mom
Angular momentum.
double * p_f_isco
Orbital frequency of the ISCO.
virtual double r_mb(int lmin, ostream *ost=0x0) const
Coordinate r of the marginally bound circular orbit (R_mb).
Scalar nn
Lapse function N .
Map & mp
Mapping describing the coordinate system (r,theta,phi)
Coord r
r coordinate centered on the grid
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 Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
int get_nzone() const
Returns the number of domains.
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
const Scalar & get_scalar(int position=0) const
Returns the reference of a Scalar stored in the list.
void add_scalar(const Scalar &ti, int position=0)
Adds the address of a new Scalar to the list.
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Tensor field of valence 0 (or component of a tensorial field).
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
const Scalar & dsdr() const
Returns of *this .
const Valeur & get_spectral_va() const
Returns va (read only version)
Values and coefficients of a (real-value) function.
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 ...
Cmp sqrt(const Cmp &)
Square root.
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.
void annule_domain(int l)
Sets the Tensor to zero in a given domain.