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