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 
32 
33 /*
34  * $Id: et_rot_isco.C,v 1.7 2016/12/05 16:17:54 j_novak Exp $
35  * $Log: et_rot_isco.C,v $
36  * Revision 1.7 2016/12/05 16:17:54 j_novak
37  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
38  *
39  * Revision 1.6 2014/10/13 08:52:58 j_novak
40  * Lorene classes and functions now belong to the namespace Lorene.
41  *
42  * Revision 1.5 2014/10/06 15:13:09 j_novak
43  * Modified #include directives to use c++ syntax.
44  *
45  * Revision 1.4 2014/07/04 12:09:06 j_novak
46  * New argument in zerosec(): a boolean (false by default) for aborting if the number of iteration is greater than the max.
47  *
48  * Revision 1.3 2011/01/07 18:20:08 m_bejger
49  * 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
50  *
51  * Revision 1.2 2001/12/06 15:11:43 jl_zdunik
52  * Introduction of the new function f_eq() in the class Etoile_rot
53  *
54  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
55  * LORENE
56  *
57  * Revision 2.2 2001/03/26 09:31:13 jlz
58  * New functions: espec_isco() and lspec_isco().
59  *
60  * Revision 2.1 2000/11/18 23:18:08 eric
61  * Ajout de l'argument ost a Etoile_rot::r_isco. Ajout de l'argument ost a Etoile_rot::r_isco.
62  *
63  * Revision 2.0 2000/11/18 21:10:41 eric
64  * *** empty log message ***
65  *
66  *
67  * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_isco.C,v 1.7 2016/12/05 16:17:54 j_novak Exp $
68  *
69  */
70 
71 // Headers C
72 #include <cmath>
73 
74 // Headers Lorene
75 #include "etoile.h"
76 #include "param.h"
77 #include "utilitaires.h"
78 
79 namespace Lorene {
80 double fonct_etoile_rot_isco(double, const Param&) ;
81 
82 
83 //=============================================================================
84 // r_isco()
85 //=============================================================================
86 
87 double Etoile_rot::r_isco(ostream* ost) const {
88 
89  if (p_r_isco == 0x0) { // a new computation is required
90 
91  // First and second derivatives of metric functions
92  // ------------------------------------------------
93 
94  int nzm1 = mp.get_mg()->get_nzone() - 1 ;
95  Cmp dnphi = nphi().dsdr() ;
96  dnphi.annule(nzm1) ;
97  Cmp ddnphi = dnphi.dsdr() ; // d^2/dr^2 N^phi
98 
99  Cmp tmp = nnn().dsdr() ;
100  tmp.annule(nzm1) ;
101  Cmp ddnnn = tmp.dsdr() ; // d^2/dr^2 N
102 
103  tmp = bbb().dsdr() ;
104  tmp.annule(nzm1) ;
105  Cmp ddbbb = tmp.dsdr() ; // d^2/dr^2 B
106 
107  // Constructing the velocity of a particle corotating with the star
108  // ----------------------------------------------------------------
109 
110  Cmp bdlog = bbb().dsdr() / bbb() ;
111  Cmp ndlog = nnn().dsdr() / nnn() ;
112  Cmp bsn = bbb() / nnn() ;
113 
114  Cmp r(mp) ;
115  r = mp.r ;
116 
117  Cmp r2= r*r ;
118 
119  bdlog.annule(nzm1) ;
120  ndlog.annule(nzm1) ;
121  bsn.annule(nzm1) ;
122  r2.annule(nzm1) ;
123 
124  // ucor_plus - the velocity of corotating particle on the circular orbit
125  Cmp ucor_plus = (r2 * bsn * dnphi +
126  sqrt ( r2 * r2 * bsn *bsn * dnphi * dnphi +
127  4 * r2 * bdlog * ndlog + 4 * r * ndlog) ) /
128  2 / (1 + r * bdlog ) ;
129 
130  Cmp factor_u2 = r2 * (2 * ddbbb / bbb() - 2 * bdlog * bdlog +
131  4 * bdlog * ndlog ) +
132  2 * r2 * r2 * bsn * bsn * dnphi * dnphi +
133  4 * r * ( ndlog - bdlog ) - 6 ;
134 
135  Cmp factor_u1 = 2 * r * r2 * bsn * ( 2 * ( ndlog - bdlog ) *
136  dnphi - ddnphi ) ;
137 
138  Cmp factor_u0 = - r2 * ( 2 * ddnnn / nnn() - 2 * ndlog * ndlog +
139  4 * bdlog * ndlog ) ;
140 
141  // Scalar field the zero of which will give the marginally stable orbit
142  Cmp orbit = factor_u2 * ucor_plus * ucor_plus +
143  factor_u1 * ucor_plus + factor_u0 ;
144  orbit.std_base_scal() ;
145 
146  // Search for the zero
147  // -------------------
148 
149  double r_ms, theta_ms, phi_ms, xi_ms,
150  xi_min = -1, xi_max = 1;
151  int l_ms = nzet, l ;
152  bool find_status = false ;
153 
154  const Valeur& vorbit = orbit.va ;
155 
156  for(l = nzet; l <= nzm1; l++) {
157 
158  // Preliminary location of the zero
159  // of the orbit function with an error = 0.01
160  theta_ms = M_PI / 2. ; // pi/2
161  phi_ms = 0. ;
162 
163  xi_min = -1. ;
164  xi_max = 1. ;
165 
166  double resloc_old ;
167  double xi_f = xi_min;
168 
169  double resloc = vorbit.val_point(l, xi_min, theta_ms, phi_ms) ;
170 
171  for (int iloc=0; iloc<200; iloc++) {
172  xi_f = xi_f + 0.01 ;
173  resloc_old = resloc ;
174  resloc = vorbit.val_point(l, xi_f, theta_ms, phi_ms) ;
175  if ( resloc * resloc_old < double(0) ) {
176  xi_min = xi_f - 0.01 ;
177  xi_max = xi_f ;
178  l_ms = l ;
179  find_status = true ;
180  break ;
181  }
182 
183  }
184 
185  }
186 
187  Param par_ms ;
188  par_ms.add_int(l_ms) ;
189  par_ms.add_cmp(orbit) ;
190 
191  if(find_status) {
192 
193  double precis_ms = 1.e-13 ; // precision in the determination of xi_ms
194  int nitermax_ms = 200 ; // max number of iterations
195 
196  int niter ;
197  xi_ms = zerosec(fonct_etoile_rot_isco, par_ms, xi_min, xi_max,
198  precis_ms, nitermax_ms, niter, false) ;
199 
200  if (niter > nitermax_ms) {
201  cerr << "Etoile_rot::r_isco : " << endl ;
202  cerr << "result may be wrong ... " << endl ;
203  cerr << "Continuing." << endl ;
204  }
205 
206  if (ost != 0x0) {
207  * ost <<
208  " number of iterations used in zerosec to locate the ISCO : "
209  << niter << endl ;
210  *ost << " zero found at xi = " << xi_ms << endl ;
211  }
212 
213  r_ms = mp.val_r(l_ms, xi_ms, theta_ms, phi_ms) ;
214 
215  } else {
216 
217  // assuming the ISCO is under the surface of a star
218  r_ms = ray_eq() ;
219  xi_ms = -1 ;
220  l_ms = nzet ;
221 
222  }
223 
224  p_r_isco = new double ((bbb().va).val_point(l_ms, xi_ms, theta_ms, phi_ms)
225  * r_ms ) ;
226 
227  // Determination of the frequency at the marginally stable orbit
228  // -------------------------------------------------------------
229 
230  ucor_plus.std_base_scal() ;
231  double ucor_msplus = (ucor_plus.va).val_point(l_ms, xi_ms, theta_ms,
232  phi_ms) ;
233  double nobrs = (bsn.va).val_point(l_ms, xi_ms, theta_ms, phi_ms) ;
234  double nphirs = (nphi().va).val_point(l_ms, xi_ms, theta_ms, phi_ms) ;
235 
236  p_f_isco = new double ( ( ucor_msplus / nobrs / r_ms + nphirs ) /
237  (double(2) * M_PI) ) ;
238 
239  // Specific angular momentum on ms orbit
240  // -------------------------------------
241  p_lspec_isco=new double (ucor_msplus/sqrt(1.-ucor_msplus*ucor_msplus)*
242  ((bbb().va).val_point(l_ms, xi_ms, theta_ms, phi_ms)) * r_ms );
243 
244  // Specific energy on ms orbit
245  // ---------------------------
246  p_espec_isco=new double (( 1./nobrs / r_ms / ucor_msplus + nphirs) *
247  (ucor_msplus/sqrt(1.-ucor_msplus*ucor_msplus)*
248  ((bbb().va).val_point(l_ms, xi_ms, theta_ms, phi_ms)) * r_ms ));
249 
250  // Determination of the Keplerian frequency at the equator
251  // -------------------------------------------------------------
252 
253  double ucor_eqplus = (ucor_plus.va).val_point(l_ms, -1, theta_ms, phi_ms) ;
254  double nobeq = (bsn.va).val_point(l_ms, -1, theta_ms, phi_ms) ;
255  double nphieq = (nphi().va).val_point(l_ms, -1, theta_ms, phi_ms) ;
256 
257  p_f_eq = new double ( ( ucor_eqplus / nobeq / ray_eq() + nphieq ) /
258  (double(2) * M_PI) ) ;
259 
260  } // End of computation
261 
262  return *p_r_isco ;
263 
264 }
265 
266 //=============================================================================
267 // f_isco()
268 //=============================================================================
269 
270 double Etoile_rot::f_isco() const {
271 
272  if (p_f_isco == 0x0) { // a new computation is required
273 
274  r_isco() ; // f_isco is computed in the method r_isco()
275 
276  assert(p_f_isco != 0x0) ;
277  }
278 
279  return *p_f_isco ;
280 
281 }
282 
283 //=============================================================================
284 // lspec_isco()
285 //=============================================================================
286 
287 double Etoile_rot::lspec_isco() const {
288 
289  if (p_lspec_isco == 0x0) { // a new computation is required
290 
291  r_isco() ; // lspec_isco is computed in the method r_isco()
292 
293  assert(p_lspec_isco != 0x0) ;
294  }
295 
296  return *p_lspec_isco ;
297 
298 }
299 
300 //=============================================================================
301 // espec_isco()
302 //=============================================================================
303 
304 double Etoile_rot::espec_isco() const {
305 
306  if (p_espec_isco == 0x0) { // a new computation is required
307 
308  r_isco() ; // espec_isco is computed in the method r_isco()
309 
310  assert(p_espec_isco != 0x0) ;
311  }
312 
313  return *p_espec_isco ;
314 
315 }
316 
317 
318 //=============================================================================
319 // f_eq()
320 //=============================================================================
321 
322 double Etoile_rot::f_eq() const {
323 
324  if (p_f_eq == 0x0) { // a new computation is required
325 
326  r_isco() ; // f_eq is computed in the method r_isco()
327 
328  assert(p_f_eq != 0x0) ;
329  }
330 
331  return *p_f_eq ;
332 
333 }
334 
335 
336 //=============================================================================
337 // Function used to locate the MS orbit
338 //=============================================================================
339 
340 
341 double fonct_etoile_rot_isco(double xi, const Param& par){
342 
343  // Retrieval of the information:
344  int l_ms = par.get_int() ;
345  const Cmp& orbit = par.get_cmp() ;
346  const Valeur& vorbit = orbit.va ;
347 
348  // Value et the desired point
349  double theta = M_PI / 2. ;
350  double phi = 0 ;
351  return vorbit.val_point(l_ms, xi, theta, phi) ;
352 
353 }
354 
355 }
const Cmp & dsdr() const
Returns of *this .
Definition: cmp_deriv.C:87
const Cmp & get_cmp(int position=0) const
Returns the reference of a Cmp stored in the list.
Definition: param.C:983
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:287
double * p_r_isco
Circumferential radius of the ISCO.
Definition: etoile.h:1648
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
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:351
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
Tenseur nnn
Total lapse function.
Definition: etoile.h:512
Tenseur nphi
Metric coefficient .
Definition: etoile.h:1513
double ray_eq() const
Coordinate radius at , [r_unit].
Values and coefficients of a (real-value) function.
Definition: valeur.h:297
virtual double r_isco(ostream *ost=0x0) const
Circumferential radius of the innermost stable circular orbit (ISCO).
Definition: et_rot_isco.C:87
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 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
void add_cmp(const Cmp &ti, int position=0)
Adds the address of a new Cmp to the list.
Definition: param.C:938
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
Parameter storage.
Definition: param.h:125
Map & mp
Mapping associated with the star.
Definition: etoile.h:432
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
virtual double espec_isco() const
Energy of a particle on the ISCO.
Definition: et_rot_isco.C:304
double * p_f_isco
Orbital frequency of the ISCO.
Definition: etoile.h:1649
virtual double f_eq() const
Orbital frequency at the equator.
Definition: et_rot_isco.C:322
Tenseur bbb
Metric factor B.
Definition: etoile.h:1507
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:647
int nzet
Number of domains of *mp occupied by the star.
Definition: etoile.h:435
double * p_lspec_isco
Specific angular momentum of a particle on the ISCO.
Definition: etoile.h:1653
virtual double f_isco() const
Orbital frequency at the innermost stable circular orbit (ISCO).
Definition: et_rot_isco.C:270
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:1651
Coord r
r coordinate centered on the grid
Definition: map.h:736
double * p_f_eq
Orbital frequency at the equator.
Definition: etoile.h:1654