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 char 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 
64 namespace Lorene {
65 double funct_star_rot_isco(double, const Param& ) ;
66 
67 //=============================================================================
68 // r_isco()
69 //=============================================================================
70 
71 double Star_rot::r_isco(ostream* ost) const {
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 +
127  factor_u1 * ucor_plus + factor_u0 ;
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 ;
157  resloc_old = resloc ;
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,
183  precis_ms, nitermax_ms, niter) ;
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 
191  r_ms = mp.val_r(l_ms, xi_ms, theta_ms, phi_ms) ;
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  // -------------------------------------
220  p_lspec_isco=new double (ucor_msplus/sqrt(1.-ucor_msplus*ucor_msplus)*
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  // ---------------------------
225  p_espec_isco=new double (( 1./nobrs / r_ms / ucor_msplus + nphirs) *
226  (ucor_msplus/sqrt(1.-ucor_msplus*ucor_msplus)*
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 
255 double 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 
272 double 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 
289 double 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 
307 double 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 
326 double 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 }
Lorene::Star_rot::p_lspec_isco
double * p_lspec_isco
Specific angular momentum of a particle on the ISCO.
Definition: star_rot.h:248
Lorene::Star_rot::f_eq
virtual double f_eq() const
Orbital frequency at the equator.
Definition: star_rot_isco.C:307
Lorene::Star_rot::r_isco
virtual double r_isco(ostream *ost=0x0) const
Circumferential radius of the innermost stable circular orbit (ISCO).
Definition: star_rot_isco.C:71
Lorene::Valeur
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
Lorene::Valeur::val_point
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
Lorene::Star::ray_eq
double ray_eq() const
Coordinate radius at , [r_unit].
Definition: star_global.C:108
Lorene::Scalar::dsdr
const Scalar & dsdr() const
Returns of *this .
Definition: scalar_deriv.C:113
Lorene::Star_rot::p_espec_isco
double * p_espec_isco
Specific energy of a particle on the ISCO.
Definition: star_rot.h:246
Lorene::Star_rot::p_f_eq
double * p_f_eq
Orbital frequency at the equator.
Definition: star_rot.h:249
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Star::nzet
int nzet
Number of domains of *mp occupied by the star.
Definition: star.h:183
Lorene::Star_rot::p_r_isco
double * p_r_isco
Circumferential radius of the ISCO.
Definition: star_rot.h:243
Lorene::Star_rot::bbb
Scalar bbb
Metric factor B.
Definition: star_rot.h:107
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Star_rot::nphi
Scalar nphi
Metric coefficient .
Definition: star_rot.h:113
Lorene::Map::r
Coord r
r coordinate centered on the grid
Definition: map.h:718
Lorene::Star_rot::espec_isco
virtual double espec_isco() const
Energy of a particle on the ISCO.
Definition: star_rot_isco.C:289
Lorene::Map::val_r
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.
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::Star::nn
Scalar nn
Lapse function N .
Definition: star.h:225
Lorene::Tensor::annule_domain
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:666
Lorene::zerosec
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::Star_rot::lspec_isco
virtual double lspec_isco() const
Angular momentum of a particle on the ISCO.
Definition: star_rot_isco.C:272
Lorene::Scalar::get_spectral_va
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:601
Lorene::Star_rot::p_f_isco
double * p_f_isco
Orbital frequency of the ISCO.
Definition: star_rot.h:244
Lorene::Star::mp
Map & mp
Mapping associated with the star.
Definition: star.h:180
Lorene::Param::add_scalar
void add_scalar(const Scalar &ti, int position=0)
Adds the address of a new Scalar to the list.
Definition: param.C:1348
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Lorene::Star_rot::f_isco
virtual double f_isco() const
Orbital frequency at the innermost stable circular orbit (ISCO).
Definition: star_rot_isco.C:255
Lorene::Scalar::std_spectral_base
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition: scalar.C:784
Lorene::Param
Parameter storage.
Definition: param.h:125
Lorene::Param::get_int
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:292
Lorene::Param::add_int
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:246
Lorene::Param::get_scalar
const Scalar & get_scalar(int position=0) const
Returns the reference of a Scalar stored in the list.
Definition: param.C:1393