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