LORENE
compobj_QI_global.C
1 /*
2  * Method of class Compobj_QI to compute the location of the ISCO
3  *
4  * (see file compobj.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2012 Odele Straub, Claire Some, Eric Gourgoulhon
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License as published by
15  * the Free Software Foundation; either version 2 of the License, or
16  * (at your option) any later version.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 
31 
32 /*
33  * $Id: compobj_QI_global.C,v 1.12 2016/12/05 16:17:49 j_novak Exp $
34  * $Log: compobj_QI_global.C,v $
35  * Revision 1.12 2016/12/05 16:17:49 j_novak
36  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
37  *
38  * Revision 1.11 2014/10/13 08:52:49 j_novak
39  * Lorene classes and functions now belong to the namespace Lorene.
40  *
41  * Revision 1.10 2014/10/06 15:13:04 j_novak
42  * Modified #include directives to use c++ syntax.
43  *
44  * Revision 1.9 2014/02/12 16:46:54 o_straub
45  * Rmb : cleaner prompt
46  *
47  * Revision 1.8 2014/01/14 20:52:53 e_gourgoulhon
48  * ISCO searched downwards
49  *
50  * Revision 1.7 2014/01/14 16:35:46 e_gourgoulhon
51  * Changed output printing in ISCO search
52  *
53  * Revision 1.6 2013/11/13 11:20:01 j_novak
54  * Minor correction to compile with older versions of g++
55  *
56  * Revision 1.5 2013/07/25 19:44:11 o_straub
57  * calculation of the marginally bound radius
58  *
59  * Revision 1.4 2013/04/04 15:32:32 e_gourgoulhon
60  * Better computation of the ISCO
61  *
62  * Revision 1.3 2012/11/22 16:04:51 c_some
63  * Minor modifications
64  *
65  * Revision 1.2 2012/11/21 14:53:45 c_some
66  * corrected mom_euler
67  *
68  * Revision 1.1 2012/11/16 16:14:11 c_some
69  * New class Compobj_QI
70  *
71  *
72  * $Header: /cvsroot/Lorene/C++/Source/Compobj/compobj_QI_global.C,v 1.12 2016/12/05 16:17:49 j_novak Exp $
73  *
74  */
75 
76 // Headers C
77 #include <cmath>
78 
79 // Headers Lorene
80 #include "compobj.h"
81 #include "param.h"
82 #include "utilitaires.h"
83 
84 namespace Lorene {
85 double funct_compobj_QI_isco(double, const Param&) ;
86 double funct_compobj_QI_rmb(double, const Param&) ;
87 
88 
89  //----------------------------//
90  // Angular momentum //
91  //----------------------------//
92 
93 double Compobj_QI::angu_mom() const {
94 
95  if (p_angu_mom == 0x0) { // a new computation is required
96 
97  cerr << "Compobj_QI::angu_mom() : not implemented yet !" << endl ; //## provisory
98  abort() ;
99  }
100 
101  return *p_angu_mom ;
102 
103 }
104 
105 
106 
107 //=============================================================================
108 // r_isco()
109 //=============================================================================
110 
111 double Compobj_QI::r_isco(int lmin, ostream* ost) const {
112 
113  if (p_r_isco == 0x0) { // a new computation is required
114 
115  // First and second derivatives of metric functions
116  // ------------------------------------------------
117 
118  int nzm1 = mp.get_mg()->get_nzone() - 1 ;
119  Scalar dnphi = nphi.dsdr() ;
120  dnphi.annule_domain(nzm1) ;
121  Scalar ddnphi = dnphi.dsdr() ; // d^2/dr^2 N^phi
122 
123  Scalar tmp = nn.dsdr() ;
124  tmp.annule_domain(nzm1) ;
125  Scalar ddnnn = tmp.dsdr() ; // d^2/dr^2 N
126 
127  tmp = bbb.dsdr() ;
128  tmp.annule_domain(nzm1) ;
129  Scalar ddbbb = tmp.dsdr() ; // d^2/dr^2 B
130 
131  // Constructing the velocity of a particle corotating with the star
132  // ----------------------------------------------------------------
133 
134  Scalar bdlog = bbb.dsdr() / bbb ;
135  Scalar ndlog = nn.dsdr() / nn ;
136  Scalar bsn = bbb / nn ;
137 
138  Scalar r(mp) ;
139  r = mp.r ;
140 
141  Scalar r2= r*r ;
142 
143  bdlog.annule_domain(nzm1) ;
144  ndlog.annule_domain(nzm1) ;
145  bsn.annule_domain(nzm1) ;
146  r2.annule_domain(nzm1) ;
147 
148  // ucor_plus - the velocity of corotating particle on the circular orbit
149  Scalar ucor_plus = (r2 * bsn * dnphi +
150  sqrt ( r2 * r2 * bsn *bsn * dnphi * dnphi +
151  4 * r2 * bdlog * ndlog + 4 * r * ndlog) ) /
152  2 / (1 + r * bdlog ) ;
153 
154  Scalar factor_u2 = r2 * (2 * ddbbb / bbb - 2 * bdlog * bdlog +
155  4 * bdlog * ndlog ) +
156  2 * r2 * r2 * bsn * bsn * dnphi * dnphi +
157  4 * r * ( ndlog - bdlog ) - 6 ;
158 
159  Scalar factor_u1 = 2 * r * r2 * bsn * ( 2 * ( ndlog - bdlog ) *
160  dnphi - ddnphi ) ;
161 
162  Scalar factor_u0 = - r2 * ( 2 * ddnnn / nn - 2 * ndlog * ndlog +
163  4 * bdlog * ndlog ) ;
164 
165  // Scalar field the zero of which will give the marginally stable orbit
166  Scalar orbit = factor_u2 * ucor_plus * ucor_plus +
167  factor_u1 * ucor_plus + factor_u0 ;
168  orbit.std_spectral_base() ;
169 
170  // Search for the zero
171  // -------------------
172 
173  double r_ms, theta_ms, phi_ms, xi_ms,
174  xi_min = -1, xi_max = 1;
175 
176  int l_ms = lmin, l ;
177  bool find_status = false ;
178 
179  const Valeur& vorbit = orbit.get_spectral_va() ;
180 
181  // Preliminary location of the zero
182  // of the orbit function with an error = 0.01
183  theta_ms = M_PI / 2. ; // pi/2
184  phi_ms = 0. ;
185 
186  for(l = nzm1-1; l >= lmin; l--) {
187 
188  xi_min = -1. ;
189  xi_max = 1. ;
190 
191  double resloc_old ;
192  double xi_f = xi_min;
193 
194  double resloc = vorbit.val_point(l, xi_min, theta_ms, phi_ms) ;
195 
196  for (int iloc=0; iloc<200; iloc++) {
197  xi_f = xi_f + 0.01 ;
198  resloc_old = resloc ;
199  resloc = vorbit.val_point(l, xi_f, theta_ms, phi_ms) ;
200  if ( resloc * resloc_old < double(0) ) {
201  xi_min = xi_f - 0.01 ;
202  xi_max = xi_f ;
203  l_ms = l ;
204  find_status = true ;
205  break ;
206  }
207 
208  }
209  if (find_status) break ;
210  }
211 
212  Param par_ms ;
213  par_ms.add_int(l_ms) ;
214  par_ms.add_scalar(orbit) ;
215 
216 
217 
218  if(find_status) {
219 
220 
221  double precis_ms = 1.e-12 ; // precision in the determination of xi_ms
222 
223  int nitermax_ms = 100 ; // max number of iterations
224 
225  int niter ;
226  xi_ms = zerosec(funct_compobj_QI_isco, par_ms, xi_min, xi_max,
227  precis_ms, nitermax_ms, niter) ;
228  if (ost != 0x0) {
229  *ost << "ISCO search: " << endl ;
230  *ost << " Domain number: " << l_ms << endl ;
231  *ost << " xi_min, xi_max : " << xi_min << " , " << xi_max << endl ;
232  *ost << " number of iterations used in zerosec: " << niter << endl ;
233  *ost << " zero found at xi = " << xi_ms << endl ;
234  }
235 
236  r_ms = mp.val_r(l_ms, xi_ms, theta_ms, phi_ms) ;
237 
238  } else {
239 
240  // ISCO not found
241  r_ms = -1 ;
242  xi_ms = -1 ;
243  l_ms = lmin ;
244 
245  }
246 
247  p_r_isco = new double (r_ms) ;
248 
249 // p_r_isco = new double (
250 // (bbb.get_spectral_va()).val_point(l_ms, xi_ms, theta_ms, phi_ms) * r_ms
251 // ) ;
252 
253  // Determination of the frequency at the marginally stable orbit
254  // -------------------------------------------------------------
255 
256  ucor_plus.std_spectral_base() ;
257  double ucor_msplus = (ucor_plus.get_spectral_va()).val_point(l_ms, xi_ms, theta_ms,
258  phi_ms) ;
259  double nobrs = (bsn.get_spectral_va()).val_point(l_ms, xi_ms, theta_ms, phi_ms) ;
260  double nphirs = (nphi.get_spectral_va()).val_point(l_ms, xi_ms, theta_ms, phi_ms) ;
261 
262  p_f_isco = new double ( ( ucor_msplus / nobrs / r_ms + nphirs ) /
263  (double(2) * M_PI) ) ;
264 
265  // Specific angular momentum on ms orbit
266  // -------------------------------------
267  p_lspec_isco=new double (ucor_msplus/sqrt(1.-ucor_msplus*ucor_msplus)*
268  ((bbb.get_spectral_va()).val_point(l_ms, xi_ms, theta_ms, phi_ms)) * r_ms );
269 
270  // Specific energy on ms orbit
271  // ---------------------------
272  p_espec_isco=new double (( 1./nobrs / r_ms / ucor_msplus + nphirs) *
273  (ucor_msplus/sqrt(1.-ucor_msplus*ucor_msplus)*
274  ((bbb.get_spectral_va()).val_point(l_ms, xi_ms, theta_ms, phi_ms)) * r_ms ));
275 
276 
277  } // End of computation
278 
279  return *p_r_isco ;
280 
281 }
282 
283 
284 
285 //=============================================================================
286 // f_isco()
287 //=============================================================================
288 
289 double Compobj_QI::f_isco(int lmin) const {
290 
291  if (p_f_isco == 0x0) { // a new computation is required
292 
293  r_isco(lmin) ; // f_isco is computed in the method r_isco()
294 
295  assert(p_f_isco != 0x0) ;
296  }
297 
298  return *p_f_isco ;
299 
300 }
301 
302 //=============================================================================
303 // lspec_isco()
304 //=============================================================================
305 
306 double Compobj_QI::lspec_isco(int lmin) const {
307 
308  if (p_lspec_isco == 0x0) { // a new computation is required
309 
310  r_isco(lmin) ; // lspec_isco is computed in the method r_isco()
311 
312  assert(p_lspec_isco != 0x0) ;
313  }
314 
315  return *p_lspec_isco ;
316 
317 }
318 
319 //=============================================================================
320 // espec_isco()
321 //=============================================================================
322 
323 double Compobj_QI::espec_isco(int lmin) const {
324 
325  if (p_espec_isco == 0x0) { // a new computation is required
326 
327  r_isco(lmin) ; // espec_isco is computed in the method r_isco()
328 
329  assert(p_espec_isco != 0x0) ;
330  }
331 
332  return *p_espec_isco ;
333 
334 }
335 
336 
337 
338 
339 
340 //=============================================================================
341 // r_mb()
342 //=============================================================================
343 
344 double Compobj_QI::r_mb(int lmin, ostream* ost) const {
345 
346  if (p_r_mb == 0x0) { // a new computation is required
347 
348  // Coefficients of the effective potential (A) and its derivative (B)
349  // ------------------------------------------------
350 
351  int nzm1 = mp.get_mg()->get_nzone() - 1 ;
352  Scalar r(mp) ;
353  r = mp.r ;
354  Scalar r2 = r*r ;
355  r2.annule_domain(nzm1) ;
356 
357  Scalar ndn = nn*nn.dsdr() ;
358  ndn.annule_domain(nzm1) ;
359 
360 
361  // Scalar V_eff = A1 + A2 * E^2 + A3 * E * L + A4 * L^2 ;
362  // Scalar dV_eff = B1 + B2 * E^2 + B3 * E * L + B4 * L^2 ;
363 
364  Scalar A1 =-bbb*bbb * nn*nn * r2 ;
365  Scalar A2 = bbb*bbb * r2 ;
366  Scalar A3 =-2. * bbb*bbb * r2 * nphi ;
367  Scalar A4 =-nn*nn + bbb*bbb * r2 * nphi*nphi ;
368 
369  Scalar B1 =-2.*r * bbb*bbb * nn*nn - 2.*r2 * bbb*bbb.dsdr() * nn*nn - 2.*r2 * bbb*bbb * nn*nn.dsdr() ;
370  Scalar B2 = 2.*r * bbb*bbb + 2.*r2 * bbb*bbb.dsdr() ;
371  Scalar B3 =-2.*nphi*B2 - 2.*r2 * bbb*bbb * nphi.dsdr() ;
372  Scalar B4 = 2.*r * bbb*bbb * nphi*nphi + 2.*r2 * bbb*bbb.dsdr() * nphi*nphi - 2.*ndn + 2.*r2 * bbb*bbb * nphi*nphi.dsdr() ;
373 
374 
375  Scalar C1 = (A1 * B3 - A3 * B1) ;
376  Scalar C2 = (A2 * B3 - A3 * B2) ;
377  Scalar C3 = (A4 * B3 - A3 * B4) ;
378 
379 
380  Scalar D1 = B4 * C1 - B1 * C3 ;
381  Scalar D2 = B4 * C2 - B2 * C3 ;
382  Scalar D3 = B3 * B3 * C1 * C3 ;
383  Scalar D4 = B3 * B3 * C2 * C3 ;
384 
385 
386 
387 
388 
389 
390  // Constructing the orbital energy of a particle corotating with the star
391  // ----------------------------------------------------------------
392 
393  /* B3 * V_eff - A3 * dV_eff = 0. ; // solve eq. for L
394 
395  Scalar L = sqrt((C1 + C2 * E2) / C3) ; // substitute into the eq. dVeff=0, then solve for E
396 
397  Scalar EE = (-(2.*D1*D2 + D3) + sqrt((2.*D1*D2 + D3) * (2.*D1*D2 + D3) - 4.*D1*D1 * (D2*D2 +
398  D4))) / (2.*(D2*D2 + D4)) ; // solve eq. EE = 1 for r
399  */
400 
401 
402  Scalar bound_orbit = -(2.*D1*D2 + D3) - sqrt((2.*D1*D2 + D3) * (2.*D1*D2 + D3) - 4.*D1*D1 *
403  (D2*D2 + D4)) - 2.*(D2*D2 + D4) ;
404 
405 
406  // cout << "bound_orbit :" << bound_orbit << endl ;
407 
408  bound_orbit.std_spectral_base() ;
409 
410 
411 
412  // Search for the zero
413  // -------------------
414 
415  const int noz(10) ; // number of zeros
416  double zeros[2][noz] ; // define array for zeros
417  int i = 0 ; // counter
418  int l ; // number of domain
419 
420  double rmb, theta_mb, phi_mb, xi_mb;
421  double xi_min = -1, xi_max = 1 ;
422 
423  const Valeur& vorbit = bound_orbit.get_spectral_va() ;
424 
425 
426  // Preliminary location of the zero
427  // of the bound_orbit function with an error = dx
428 
429  double dx = 0.001 ;
430 
431  theta_mb = M_PI / 2. ;
432  phi_mb = 0. ;
433 
434 
435  for(l = lmin; l <= nzm1; l++) {
436 
437  xi_min = -1. ;
438 
439  double resloc_old ;
440  double xi_f = xi_min;
441 
442  double resloc = vorbit.val_point(l, xi_f, theta_mb, phi_mb) ;
443 
444  while(xi_f <= xi_max) {
445 
446  xi_f = xi_f + dx ;
447 
448  resloc_old = resloc ;
449  resloc = vorbit.val_point(l, xi_f, theta_mb, phi_mb) ;
450 
451  if ( resloc * resloc_old < double(0) ) {
452 
453  zeros[0][i] = xi_f ; // xi_max
454  zeros[1][i] = l ; // domain number l
455  i++ ;
456 
457  }
458 
459  }
460 
461  }
462 
463 
464 
465  int number_of_zeros = i ;
466 
467  cout << "number of zeros: " << number_of_zeros << endl ;
468 
469 
470  double precis_mb = 1.e-9 ; // precision in the determination of xi_mb: 1.e-12
471 
472 
473  int nitermax_mb = 100 ; // max number of iterations
474 
475 
476  for(int i = 0; i < number_of_zeros; i++) {
477 
478  //cout << i << " " << zeros[0][i] << " " << zeros[1][i] << endl ;
479 
480  int l_mb = int(zeros[1][i]) ;
481  xi_max = zeros[0][i] ;
482  xi_min = xi_max - dx ;
483 
484  Param par_mb ;
485  par_mb.add_scalar(bound_orbit) ;
486  par_mb.add_int(l_mb) ;
487 
488  int niter ;
489  xi_mb = zerosec(funct_compobj_QI_rmb, par_mb, xi_min, xi_max, precis_mb, nitermax_mb, niter) ;
490  if (ost != 0x0) {
491  *ost << "RMB search: " << endl ;
492  *ost << " Domain number: " << l_mb << endl ;
493  *ost << " xi_min, xi_max : " << xi_min << " , " << xi_max << endl ;
494  *ost << " number of iterations used in zerosec: " << niter << endl ;
495  *ost << " zero found at xi = " << xi_mb << endl ;
496  }
497 
498  if (niter < nitermax_mb) {
499  double zero_mb = mp.val_r(l_mb, xi_mb, theta_mb, phi_mb) ;
500  //double r_hor = radius_hor(0); // set to 1 in the condition below
501  double r_ms = r_isco(0) ;
502  if (zero_mb < (1 + r_ms)/2){
503 
504  rmb = zero_mb ;
505  }
506  }
507  }
508 
509  p_r_mb = new double (rmb) ;
510 
511  //delete [] zeros ; not used, causes "core dump" in Code kerr_qi
512 
513  } // End of computation
514 
515 
516  return *p_r_mb ;
517 
518 }
519 
520 
521 
522 
523 
524 
525 
526 //=============================================================================
527 // Function used to locate the MS orbit
528 //=============================================================================
529 
530 
531 double funct_compobj_QI_isco(double xi, const Param& par){
532 
533  // Retrieval of the information:
534  int l_ms = par.get_int() ;
535  const Scalar& orbit = par.get_scalar() ;
536  const Valeur& vorbit = orbit.get_spectral_va() ;
537 
538  // Value et the desired point
539  double theta = M_PI / 2. ;
540  double phi = 0 ;
541  return vorbit.val_point(l_ms, xi, theta, phi) ;
542 
543 }
544 
545 
546 
547 //=============================================================================
548 // Function used to locate the MB orbit
549 //=============================================================================
550 
551 double funct_compobj_QI_rmb(double zeros, const Param& par){
552 
553  // Retrieval of the information:
554  int l_mb = par.get_int() ;
555  const Scalar& orbit = par.get_scalar() ;
556  const Valeur& vorbit = orbit.get_spectral_va() ;
557 
558  // Value et the desired point
559  double theta = M_PI / 2. ;
560  double phi = 0 ;
561  return vorbit.val_point(l_mb, zeros, theta, phi) ;
562 
563 }
564 
565 
566 
567 
568 }
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:675
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:249
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
Lorene prototypes.
Definition: app_hor.h:67
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
virtual double lspec_isco(int lmin) const
Angular momentum of a particle at the ISCO.
double * p_lspec_isco
Specific angular momentum of a particle at the ISCO.
Definition: compobj.h:330
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
virtual double espec_isco(int lmin) const
Energy of a particle at the ISCO.
virtual double angu_mom() const
Angular momentum.
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.
Definition: zerosec.C:92
virtual double f_isco(int lmin) const
Orbital frequency at the innermost stable circular orbit (ISCO).
Scalar nphi
Metric coefficient .
Definition: compobj.h:299
double * p_angu_mom
Angular momentum.
Definition: compobj.h:324
const Scalar & get_scalar(int position=0) const
Returns the reference of a Scalar stored in the list.
Definition: param.C:1396
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.
Definition: param.C:295
virtual double r_isco(int lmin, ostream *ost=0x0) const
Coordinate r of the innermost stable circular orbit (ISCO).
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 ...
Definition: valeur.C:885
Scalar bbb
Metric factor B.
Definition: compobj.h:293
double * p_espec_isco
Specific energy of a particle at the ISCO.
Definition: compobj.h:328
Parameter storage.
Definition: param.h:125
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
double * p_r_isco
Coordinate r of the ISCO.
Definition: compobj.h:325
double * p_r_mb
Coordinate r of the marginally bound orbit.
Definition: compobj.h:331
void add_scalar(const Scalar &ti, int position=0)
Adds the address of a new Scalar to the list.
Definition: param.C:1351
double * p_f_isco
Orbital frequency of the ISCO.
Definition: compobj.h:326
Scalar nn
Lapse function N .
Definition: compobj.h:138
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
virtual double r_mb(int lmin, ostream *ost=0x0) const
Coordinate r of the marginally bound circular orbit (R_mb).
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607
Map & mp
Mapping describing the coordinate system (r,theta,phi)
Definition: compobj.h:135
Coord r
r coordinate centered on the grid
Definition: map.h:736