LORENE
et_rot_isco.C
1/*
2 * Method of class Etoile_rot to compute the location of the ISCO
3 *
4 * (see file etoile.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2000-2001 Eric Gourgoulhon
10 * Copyright (c) 2000-2001 J. Leszek Zdunik
11 *
12 * This file is part of LORENE.
13 *
14 * LORENE is free software; you can redistribute it and/or modify
15 * it under the terms of the GNU General Public License as published by
16 * the Free Software Foundation; either version 2 of the License, or
17 * (at your option) any later version.
18 *
19 * LORENE is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU General Public License for more details.
23 *
24 * You should have received a copy of the GNU General Public License
25 * along with LORENE; if not, write to the Free Software
26 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
27 *
28 */
29
30
31char et_rot_isco_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_isco.C,v 1.6 2014/10/13 08:52:58 j_novak Exp $" ;
32
33/*
34 * $Id: et_rot_isco.C,v 1.6 2014/10/13 08:52:58 j_novak Exp $
35 * $Log: et_rot_isco.C,v $
36 * Revision 1.6 2014/10/13 08:52:58 j_novak
37 * Lorene classes and functions now belong to the namespace Lorene.
38 *
39 * Revision 1.5 2014/10/06 15:13:09 j_novak
40 * Modified #include directives to use c++ syntax.
41 *
42 * Revision 1.4 2014/07/04 12:09:06 j_novak
43 * New argument in zerosec(): a boolean (false by default) for aborting if the number of iteration is greater than the max.
44 *
45 * Revision 1.3 2011/01/07 18:20:08 m_bejger
46 * Correcting for the case of stiff EOS, in which ISCO may be farther than the first domain outside the star - now searching all non-compactified domains
47 *
48 * Revision 1.2 2001/12/06 15:11:43 jl_zdunik
49 * Introduction of the new function f_eq() in the class Etoile_rot
50 *
51 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
52 * LORENE
53 *
54 * Revision 2.2 2001/03/26 09:31:13 jlz
55 * New functions: espec_isco() and lspec_isco().
56 *
57 * Revision 2.1 2000/11/18 23:18:08 eric
58 * Ajout de l'argument ost a Etoile_rot::r_isco. Ajout de l'argument ost a Etoile_rot::r_isco.
59 *
60 * Revision 2.0 2000/11/18 21:10:41 eric
61 * *** empty log message ***
62 *
63 *
64 * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_isco.C,v 1.6 2014/10/13 08:52:58 j_novak Exp $
65 *
66 */
67
68// Headers C
69#include <cmath>
70
71// Headers Lorene
72#include "etoile.h"
73#include "param.h"
74#include "utilitaires.h"
75
76namespace Lorene {
77double fonct_etoile_rot_isco(double, const Param&) ;
78
79
80//=============================================================================
81// r_isco()
82//=============================================================================
83
84double Etoile_rot::r_isco(ostream* ost) const {
85
86 if (p_r_isco == 0x0) { // a new computation is required
87
88 // First and second derivatives of metric functions
89 // ------------------------------------------------
90
91 int nzm1 = mp.get_mg()->get_nzone() - 1 ;
92 Cmp dnphi = nphi().dsdr() ;
93 dnphi.annule(nzm1) ;
94 Cmp ddnphi = dnphi.dsdr() ; // d^2/dr^2 N^phi
95
96 Cmp tmp = nnn().dsdr() ;
97 tmp.annule(nzm1) ;
98 Cmp ddnnn = tmp.dsdr() ; // d^2/dr^2 N
99
100 tmp = bbb().dsdr() ;
101 tmp.annule(nzm1) ;
102 Cmp ddbbb = tmp.dsdr() ; // d^2/dr^2 B
103
104 // Constructing the velocity of a particle corotating with the star
105 // ----------------------------------------------------------------
106
107 Cmp bdlog = bbb().dsdr() / bbb() ;
108 Cmp ndlog = nnn().dsdr() / nnn() ;
109 Cmp bsn = bbb() / nnn() ;
110
111 Cmp r(mp) ;
112 r = mp.r ;
113
114 Cmp r2= r*r ;
115
116 bdlog.annule(nzm1) ;
117 ndlog.annule(nzm1) ;
118 bsn.annule(nzm1) ;
119 r2.annule(nzm1) ;
120
121 // ucor_plus - the velocity of corotating particle on the circular orbit
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 ) ;
126
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 ;
131
132 Cmp factor_u1 = 2 * r * r2 * bsn * ( 2 * ( ndlog - bdlog ) *
133 dnphi - ddnphi ) ;
134
135 Cmp factor_u0 = - r2 * ( 2 * ddnnn / nnn() - 2 * ndlog * ndlog +
136 4 * bdlog * ndlog ) ;
137
138 // Scalar field the zero of which will give the marginally stable orbit
139 Cmp orbit = factor_u2 * ucor_plus * ucor_plus +
140 factor_u1 * ucor_plus + factor_u0 ;
141 orbit.std_base_scal() ;
142
143 // Search for the zero
144 // -------------------
145
146 double r_ms, theta_ms, phi_ms, xi_ms,
147 xi_min = -1, xi_max = 1;
148 int l_ms = nzet, l ;
149 bool find_status = false ;
150
151 const Valeur& vorbit = orbit.va ;
152
153 for(l = nzet; l <= nzm1; l++) {
154
155 // Preliminary location of the zero
156 // of the orbit function with an error = 0.01
157 theta_ms = M_PI / 2. ; // pi/2
158 phi_ms = 0. ;
159
160 xi_min = -1. ;
161 xi_max = 1. ;
162
163 double resloc_old ;
164 double xi_f = xi_min;
165
166 double resloc = vorbit.val_point(l, xi_min, theta_ms, phi_ms) ;
167
168 for (int iloc=0; iloc<200; iloc++) {
169 xi_f = xi_f + 0.01 ;
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 ;
174 xi_max = xi_f ;
175 l_ms = l ;
176 find_status = true ;
177 break ;
178 }
179
180 }
181
182 }
183
184 Param par_ms ;
185 par_ms.add_int(l_ms) ;
186 par_ms.add_cmp(orbit) ;
187
188 if(find_status) {
189
190 double precis_ms = 1.e-13 ; // precision in the determination of xi_ms
191 int nitermax_ms = 200 ; // max number of iterations
192
193 int niter ;
194 xi_ms = zerosec(fonct_etoile_rot_isco, par_ms, xi_min, xi_max,
195 precis_ms, nitermax_ms, niter, false) ;
196
197 if (niter > nitermax_ms) {
198 cerr << "Etoile_rot::r_isco : " << endl ;
199 cerr << "result may be wrong ... " << endl ;
200 cerr << "Continuing." << endl ;
201 }
202
203 if (ost != 0x0) {
204 * ost <<
205 " number of iterations used in zerosec to locate the ISCO : "
206 << niter << endl ;
207 *ost << " zero found at xi = " << xi_ms << endl ;
208 }
209
210 r_ms = mp.val_r(l_ms, xi_ms, theta_ms, phi_ms) ;
211
212 } else {
213
214 // assuming the ISCO is under the surface of a star
215 r_ms = ray_eq() ;
216 xi_ms = -1 ;
217 l_ms = nzet ;
218
219 }
220
221 p_r_isco = new double ((bbb().va).val_point(l_ms, xi_ms, theta_ms, phi_ms)
222 * r_ms ) ;
223
224 // Determination of the frequency at the marginally stable orbit
225 // -------------------------------------------------------------
226
227 ucor_plus.std_base_scal() ;
228 double ucor_msplus = (ucor_plus.va).val_point(l_ms, xi_ms, theta_ms,
229 phi_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) ;
232
233 p_f_isco = new double ( ( ucor_msplus / nobrs / r_ms + nphirs ) /
234 (double(2) * M_PI) ) ;
235
236 // Specific angular momentum on ms orbit
237 // -------------------------------------
238 p_lspec_isco=new double (ucor_msplus/sqrt(1.-ucor_msplus*ucor_msplus)*
239 ((bbb().va).val_point(l_ms, xi_ms, theta_ms, phi_ms)) * r_ms );
240
241 // Specific energy on ms orbit
242 // ---------------------------
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 ));
246
247 // Determination of the Keplerian frequency at the equator
248 // -------------------------------------------------------------
249
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) ;
253
254 p_f_eq = new double ( ( ucor_eqplus / nobeq / ray_eq() + nphieq ) /
255 (double(2) * M_PI) ) ;
256
257 } // End of computation
258
259 return *p_r_isco ;
260
261}
262
263//=============================================================================
264// f_isco()
265//=============================================================================
266
267double Etoile_rot::f_isco() const {
268
269 if (p_f_isco == 0x0) { // a new computation is required
270
271 r_isco() ; // f_isco is computed in the method r_isco()
272
273 assert(p_f_isco != 0x0) ;
274 }
275
276 return *p_f_isco ;
277
278}
279
280//=============================================================================
281// lspec_isco()
282//=============================================================================
283
285
286 if (p_lspec_isco == 0x0) { // a new computation is required
287
288 r_isco() ; // lspec_isco is computed in the method r_isco()
289
290 assert(p_lspec_isco != 0x0) ;
291 }
292
293 return *p_lspec_isco ;
294
295}
296
297//=============================================================================
298// espec_isco()
299//=============================================================================
300
302
303 if (p_espec_isco == 0x0) { // a new computation is required
304
305 r_isco() ; // espec_isco is computed in the method r_isco()
306
307 assert(p_espec_isco != 0x0) ;
308 }
309
310 return *p_espec_isco ;
311
312}
313
314
315//=============================================================================
316// f_eq()
317//=============================================================================
318
319double Etoile_rot::f_eq() const {
320
321 if (p_f_eq == 0x0) { // a new computation is required
322
323 r_isco() ; // f_eq is computed in the method r_isco()
324
325 assert(p_f_eq != 0x0) ;
326 }
327
328 return *p_f_eq ;
329
330}
331
332
333//=============================================================================
334// Function used to locate the MS orbit
335//=============================================================================
336
337
338double fonct_etoile_rot_isco(double xi, const Param& par){
339
340 // Retrieval of the information:
341 int l_ms = par.get_int() ;
342 const Cmp& orbit = par.get_cmp() ;
343 const Valeur& vorbit = orbit.va ;
344
345 // Value et the desired point
346 double theta = M_PI / 2. ;
347 double phi = 0 ;
348 return vorbit.val_point(l_ms, xi, theta, phi) ;
349
350}
351
352}
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Valeur va
The numerical value of the Cmp
Definition cmp.h:464
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition cmp.C:644
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:348
const Cmp & dsdr() const
Returns of *this .
Definition cmp_deriv.C:84
virtual double espec_isco() const
Energy of a particle on the ISCO.
virtual double f_isco() const
Orbital frequency at the innermost stable circular orbit (ISCO).
Tenseur nphi
Metric coefficient .
Definition etoile.h:1510
virtual double r_isco(ostream *ost=0x0) const
Circumferential radius of the innermost stable circular orbit (ISCO).
Definition et_rot_isco.C:84
Tenseur bbb
Metric factor B.
Definition etoile.h:1504
virtual double lspec_isco() const
Angular momentum of a particle on the ISCO.
double * p_f_eq
Orbital frequency at the equator.
Definition etoile.h:1650
double * p_f_isco
Orbital frequency of the ISCO.
Definition etoile.h:1645
double * p_lspec_isco
Specific angular momentum of a particle on the ISCO.
Definition etoile.h:1649
double * p_espec_isco
Specific energy of a particle on the ISCO.
Definition etoile.h:1647
double * p_r_isco
Circumferential radius of the ISCO.
Definition etoile.h:1644
virtual double f_eq() const
Orbital frequency at the equator.
int nzet
Number of domains of *mp occupied by the star.
Definition etoile.h:432
double ray_eq() const
Coordinate radius at , [r_unit].
Tenseur nnn
Total lapse function.
Definition etoile.h:509
Map & mp
Mapping associated with the star.
Definition etoile.h:429
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
void add_cmp(const Cmp &ti, int position=0)
Adds the address of a new Cmp to the list.
Definition param.C:935
const Cmp & get_cmp(int position=0) const
Returns the reference of a Cmp stored in the list.
Definition param.C:980
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:246
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
Lorene prototypes.
Definition app_hor.h:64