LORENE
star_rot_isco.C
1/*
2 * Method of class Star_rot to compute the location of the ISCO
3 *
4 * (see file star_rot.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2010 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 star_rot_isco_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/star_rot_isco.C,v 1.5 2014/10/13 08:53:39 j_novak Exp $" ;
32
33/*
34 * $Id: star_rot_isco.C,v 1.5 2014/10/13 08:53:39 j_novak Exp $
35 * $Log: star_rot_isco.C,v $
36 * Revision 1.5 2014/10/13 08:53:39 j_novak
37 * Lorene classes and functions now belong to the namespace Lorene.
38 *
39 * Revision 1.4 2014/10/06 15:13:17 j_novak
40 * Modified #include directives to use c++ syntax.
41 *
42 * Revision 1.3 2011/01/07 18:20:21 m_bejger
43 * 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
44 *
45 * Revision 1.2 2010/01/25 22:33:04 e_gourgoulhon
46 * First implementation.
47 *
48 * Revision 1.1 2010/01/25 18:15:52 e_gourgoulhon
49 * First version.
50 *
51 *
52 * $Header: /cvsroot/Lorene/C++/Source/Star/star_rot_isco.C,v 1.5 2014/10/13 08:53:39 j_novak Exp $
53 *
54 */
55
56// Headers C
57#include <cmath>
58
59// Headers Lorene
60#include "star_rot.h"
61#include "param.h"
62#include "utilitaires.h"
63
64namespace Lorene {
65double funct_star_rot_isco(double, const Param& ) ;
66
67//=============================================================================
68// r_isco()
69//=============================================================================
70
72
73 if (p_r_isco == 0x0) { // a new computation is required
74
75 // First and second derivatives of metric functions
76 // ------------------------------------------------
77
78 int nzm1 = mp.get_mg()->get_nzone() - 1 ;
79 Scalar dnphi = nphi.dsdr() ;
80 dnphi.annule_domain(nzm1) ;
81 Scalar ddnphi = dnphi.dsdr() ; // d^2/dr^2 N^phi
82
83 Scalar tmp = nn.dsdr() ;
84 tmp.annule_domain(nzm1) ;
85 Scalar ddnnn = tmp.dsdr() ; // d^2/dr^2 N
86
87 tmp = bbb.dsdr() ;
88 tmp.annule_domain(nzm1) ;
89 Scalar ddbbb = tmp.dsdr() ; // d^2/dr^2 B
90
91 // Constructing the velocity of a particle corotating with the star
92 // ----------------------------------------------------------------
93
94 Scalar bdlog = bbb.dsdr() / bbb ;
95 Scalar ndlog = nn.dsdr() / nn ;
96 Scalar bsn = bbb / nn ;
97
98 Scalar r(mp) ;
99 r = mp.r ;
100
101 Scalar r2= r*r ;
102
103 bdlog.annule_domain(nzm1) ;
104 ndlog.annule_domain(nzm1) ;
105 bsn.annule_domain(nzm1) ;
106 r2.annule_domain(nzm1) ;
107
108 // ucor_plus - the velocity of corotating particle on the circular orbit
109 Scalar ucor_plus = (r2 * bsn * dnphi +
110 sqrt ( r2 * r2 * bsn *bsn * dnphi * dnphi +
111 4 * r2 * bdlog * ndlog + 4 * r * ndlog) ) /
112 2 / (1 + r * bdlog ) ;
113
114 Scalar factor_u2 = r2 * (2 * ddbbb / bbb - 2 * bdlog * bdlog +
115 4 * bdlog * ndlog ) +
116 2 * r2 * r2 * bsn * bsn * dnphi * dnphi +
117 4 * r * ( ndlog - bdlog ) - 6 ;
118
119 Scalar factor_u1 = 2 * r * r2 * bsn * ( 2 * ( ndlog - bdlog ) *
120 dnphi - ddnphi ) ;
121
122 Scalar factor_u0 = - r2 * ( 2 * ddnnn / nn - 2 * ndlog * ndlog +
123 4 * bdlog * ndlog ) ;
124
125 // Scalar field the zero of which will give the marginally stable orbit
126 Scalar orbit = factor_u2 * ucor_plus * ucor_plus +
128 orbit.std_spectral_base() ;
129
130 // Search for the zero
131 // -------------------
132
133 double r_ms, theta_ms, phi_ms, xi_ms,
134 xi_min = -1, xi_max = 1;
135 int l_ms = nzet, l ;
136 bool find_status = false ;
137
138 const Valeur& vorbit = orbit.get_spectral_va() ;
139
140 for(l = nzet; l <= nzm1; l++) {
141
142 // Preliminary location of the zero
143 // of the orbit function with an error = 0.01
144 theta_ms = M_PI / 2. ; // pi/2
145 phi_ms = 0. ;
146
147 xi_min = -1. ;
148 xi_max = 1. ;
149
150 double resloc_old ;
151 double xi_f = xi_min;
152
153 double resloc = vorbit.val_point(l, xi_min, theta_ms, phi_ms) ;
154
155 for (int iloc=0; iloc<200; iloc++) {
156 xi_f = xi_f + 0.01 ;
158 resloc = vorbit.val_point(l, xi_f, theta_ms, phi_ms) ;
159 if ( resloc * resloc_old < double(0) ) {
160 xi_min = xi_f - 0.01 ;
161 xi_max = xi_f ;
162 l_ms = l ;
163 find_status = true ;
164 break ;
165 }
166
167 }
168
169 }
170
171 Param par_ms ;
172 par_ms.add_int(l_ms) ;
173 par_ms.add_scalar(orbit) ;
174
175 if(find_status) {
176
177 double precis_ms = 1.e-12 ; // precision in the determination
178 // of xi_ms
179 int nitermax_ms = 100 ; // max number of iterations
180
181 int niter ;
182 xi_ms = zerosec(funct_star_rot_isco, par_ms, xi_min, xi_max,
184 if (ost != 0x0) {
185 * ost <<
186 " number of iterations used in zerosec to locate the ISCO : "
187 << niter << endl ;
188 *ost << " zero found at xi = " << xi_ms << endl ;
189 }
190
192
193 } else {
194
195 // assuming the ISCO is under the surface of a star
196 r_ms = ray_eq() ;
197 xi_ms = -1 ;
198 l_ms = nzet ;
199
200 }
201
202 p_r_isco = new double (
203 (bbb.get_spectral_va()).val_point(l_ms, xi_ms, theta_ms, phi_ms) * r_ms
204 ) ;
205
206 // Determination of the frequency at the marginally stable orbit
207 // -------------------------------------------------------------
208
209 ucor_plus.std_spectral_base() ;
210 double ucor_msplus = (ucor_plus.get_spectral_va()).val_point(l_ms, xi_ms, theta_ms,
211 phi_ms) ;
212 double nobrs = (bsn.get_spectral_va()).val_point(l_ms, xi_ms, theta_ms, phi_ms) ;
213 double nphirs = (nphi.get_spectral_va()).val_point(l_ms, xi_ms, theta_ms, phi_ms) ;
214
215 p_f_isco = new double ( ( ucor_msplus / nobrs / r_ms + nphirs ) /
216 (double(2) * M_PI) ) ;
217
218 // Specific angular momentum on ms orbit
219 // -------------------------------------
221 ((bbb.get_spectral_va()).val_point(l_ms, xi_ms, theta_ms, phi_ms)) * r_ms );
222
223 // Specific energy on ms orbit
224 // ---------------------------
227 ((bbb.get_spectral_va()).val_point(l_ms, xi_ms, theta_ms, phi_ms)) * r_ms ));
228
229 // Determination of the Keplerian frequency at the equator
230 // -------------------------------------------------------------
231
232
233 double ucor_eqplus = (ucor_plus.get_spectral_va()).val_point(l_ms, -1, theta_ms,phi_ms)
234 ;
235 double nobeq = (bsn.get_spectral_va()).val_point(l_ms, -1, theta_ms, phi_ms) ;
236 double nphieq = (nphi.get_spectral_va()).val_point(l_ms, -1, theta_ms, phi_ms) ;
237
238 p_f_eq = new double ( ( ucor_eqplus / nobeq / ray_eq() + nphieq ) /
239 (double(2) * M_PI) ) ;
240
241
242
243 } // End of computation
244
245 return *p_r_isco ;
246
247}
248
249
250
251//=============================================================================
252// f_isco()
253//=============================================================================
254
255double Star_rot::f_isco() const {
256
257 if (p_f_isco == 0x0) { // a new computation is required
258
259 r_isco() ; // f_isco is computed in the method r_isco()
260
261 assert(p_f_isco != 0x0) ;
262 }
263
264 return *p_f_isco ;
265
266}
267
268//=============================================================================
269// lspec_isco()
270//=============================================================================
271
272double Star_rot::lspec_isco() const {
273
274 if (p_lspec_isco == 0x0) { // a new computation is required
275
276 r_isco() ; // lspec_isco is computed in the method r_isco()
277
278 assert(p_lspec_isco != 0x0) ;
279 }
280
281 return *p_lspec_isco ;
282
283}
284
285//=============================================================================
286// espec_isco()
287//=============================================================================
288
289double Star_rot::espec_isco() const {
290
291 if (p_espec_isco == 0x0) { // a new computation is required
292
293 r_isco() ; // espec_isco is computed in the method r_isco()
294
295 assert(p_espec_isco != 0x0) ;
296 }
297
298 return *p_espec_isco ;
299
300}
301
302
303//=============================================================================
304// f_eq()
305//=============================================================================
306
307double Star_rot::f_eq() const {
308
309 if (p_f_eq == 0x0) { // a new computation is required
310
311 r_isco() ; // f_eq is computed in the method r_isco()
312
313 assert(p_f_eq != 0x0) ;
314 }
315
316 return *p_f_eq ;
317
318}
319
320
321//=============================================================================
322// Function used to locate the MS orbit
323//=============================================================================
324
325
326double funct_star_rot_isco(double xi, const Param& par){
327
328 // Retrieval of the information:
329 int l_ms = par.get_int() ;
330 const Scalar& orbit = par.get_scalar() ;
331 const Valeur& vorbit = orbit.get_spectral_va() ;
332
333 // Value et the desired point
334 double theta = M_PI / 2. ;
335 double phi = 0 ;
336 return vorbit.val_point(l_ms, xi, theta, phi) ;
337
338}
339
340
341
342
343
344
345}
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
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
Parameter storage.
Definition param.h:125
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
const Scalar & dsdr() const
Returns of *this .
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition scalar.h:601
virtual double lspec_isco() const
Angular momentum of a particle on the ISCO.
double * p_r_isco
Circumferential radius of the ISCO.
Definition star_rot.h:243
virtual double espec_isco() const
Energy of a particle on the ISCO.
virtual double r_isco(ostream *ost=0x0) const
Circumferential radius of the innermost stable circular orbit (ISCO).
double * p_lspec_isco
Specific angular momentum of a particle on the ISCO.
Definition star_rot.h:248
Scalar bbb
Metric factor B.
Definition star_rot.h:107
Scalar nphi
Metric coefficient .
Definition star_rot.h:113
double * p_f_isco
Orbital frequency of the ISCO.
Definition star_rot.h:244
virtual double f_isco() const
Orbital frequency at the innermost stable circular orbit (ISCO).
double * p_espec_isco
Specific energy of a particle on the ISCO.
Definition star_rot.h:246
double * p_f_eq
Orbital frequency at the equator.
Definition star_rot.h:249
virtual double f_eq() const
Orbital frequency at the equator.
Scalar nn
Lapse function N .
Definition star.h:225
double ray_eq() const
Coordinate radius at , [r_unit].
Map & mp
Mapping associated with the star.
Definition star.h:180
int nzet
Number of domains of *mp occupied by the star.
Definition star.h:183
Values and coefficients of a (real-value) function.
Definition valeur.h:287
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