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 
31 char 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 
76 namespace Lorene {
77 double fonct_etoile_rot_isco(double, const Param&) ;
78 
79 
80 //=============================================================================
81 // r_isco()
82 //=============================================================================
83 
84 double 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 
267 double 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 
284 double Etoile_rot::lspec_isco() const {
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 
301 double Etoile_rot::espec_isco() const {
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 
319 double 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 
338 double 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 }
const Cmp & dsdr() const
Returns of *this .
Definition: cmp_deriv.C:84
const Cmp & get_cmp(int position=0) const
Returns the reference of a Cmp stored in the list.
Definition: param.C:980
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
virtual double lspec_isco() const
Angular momentum of a particle on the ISCO.
Definition: et_rot_isco.C:284
double * p_r_isco
Circumferential radius of the ISCO.
Definition: etoile.h:1644
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:246
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:348
Lorene prototypes.
Definition: app_hor.h:64
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Tenseur nnn
Total lapse function.
Definition: etoile.h:509
Tenseur nphi
Metric coefficient .
Definition: etoile.h:1510
double ray_eq() const
Coordinate radius at , [r_unit].
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
virtual double r_isco(ostream *ost=0x0) const
Circumferential radius of the innermost stable circular orbit (ISCO).
Definition: et_rot_isco.C:84
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
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:292
void add_cmp(const Cmp &ti, int position=0)
Adds the address of a new Cmp to the list.
Definition: param.C:935
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
Parameter storage.
Definition: param.h:125
Map & mp
Mapping associated with the star.
Definition: etoile.h:429
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
virtual double espec_isco() const
Energy of a particle on the ISCO.
Definition: et_rot_isco.C:301
double * p_f_isco
Orbital frequency of the ISCO.
Definition: etoile.h:1645
virtual double f_eq() const
Orbital frequency at the equator.
Definition: et_rot_isco.C:319
Tenseur bbb
Metric factor B.
Definition: etoile.h:1504
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:644
int nzet
Number of domains of *mp occupied by the star.
Definition: etoile.h:432
double * p_lspec_isco
Specific angular momentum of a particle on the ISCO.
Definition: etoile.h:1649
virtual double f_isco() const
Orbital frequency at the innermost stable circular orbit (ISCO).
Definition: et_rot_isco.C:267
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
double * p_espec_isco
Specific energy of a particle on the ISCO.
Definition: etoile.h:1647
Coord r
r coordinate centered on the grid
Definition: map.h:718
double * p_f_eq
Orbital frequency at the equator.
Definition: etoile.h:1650