LORENE
et_rot_f_eccentric.C
1 /*
2  * Method of class Etoile_rot to compute eccentric orbits
3  *
4  * (see file etoile.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2001 Dorota Gondek-Rosinska
10  * Copyright (c) 2001 Eric Gourgoulhon
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_f_eccentric_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_f_eccentric.C,v 1.7 2014/10/13 08:52:57 j_novak Exp $" ;
32 
33 /*
34  * $Id: et_rot_f_eccentric.C,v 1.7 2014/10/13 08:52:57 j_novak Exp $
35  * $Log: et_rot_f_eccentric.C,v $
36  * Revision 1.7 2014/10/13 08:52:57 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.6 2014/10/06 15:13:09 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.5 2003/12/19 16:31:52 j_novak
43  * Still warnings...
44  *
45  * Revision 1.4 2003/12/19 16:21:42 j_novak
46  * Shadow hunt
47  *
48  * Revision 1.3 2003/12/05 14:50:26 j_novak
49  * To suppress some warnings...
50  *
51  * Revision 1.2 2003/10/03 15:58:47 j_novak
52  * Cleaning of some headers
53  *
54  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
55  * LORENE
56  *
57  * Revision 2.0 2001/02/08 15:13:24 eric
58  * *** empty log message ***
59  *
60  *
61  * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_rot_f_eccentric.C,v 1.7 2014/10/13 08:52:57 j_novak Exp $
62  *
63  */
64 
65 
66 // Headers C
67 #include <cmath>
68 
69 // Headers Lorene
70 #include "etoile.h"
71 #include "param.h"
72 
73 //=============================================================================
74 namespace Lorene {
75 // r_isco()
76 //=============================================================================
77 
78 double Etoile_rot::f_eccentric(double, double, ostream* ost) const {
79 
80  cout << "Etoile_rot::f_eccentric not ready yet !" << endl ;
81  abort() ;
82 
83  // First and second derivatives of metric functions
84  // ------------------------------------------------
85 
86  int nzm1 = mp.get_mg()->get_nzone() - 1 ;
87  Cmp dnphi = nphi().dsdr() ;
88  dnphi.annule(nzm1) ;
89  Cmp ddnphi = dnphi.dsdr() ; // d^2/dr^2 N^phi
90 
91  Cmp tmp = nnn().dsdr() ;
92  tmp.annule(nzm1) ;
93  Cmp ddnnn = tmp.dsdr() ; // d^2/dr^2 N
94 
95  tmp = bbb().dsdr() ;
96  tmp.annule(nzm1) ;
97  Cmp ddbbb = tmp.dsdr() ; // d^2/dr^2 B
98 
99  // Constructing the velocity of a particle corotating with the star
100  // ----------------------------------------------------------------
101 
102  Cmp bdlog = bbb().dsdr() / bbb() ;
103  Cmp ndlog = nnn().dsdr() / nnn() ;
104  Cmp bsn = bbb() / nnn() ;
105 
106  Cmp r(mp) ;
107  r = mp.r ;
108 
109  Cmp r2= r*r ;
110 
111  bdlog.annule(nzm1) ;
112  ndlog.annule(nzm1) ;
113  bsn.annule(nzm1) ;
114  r2.annule(nzm1) ;
115 
116  // ucor_plus - the velocity of corotating particle on the circular orbit
117  Cmp ucor_plus = (r2 * bsn * dnphi +
118  sqrt ( r2 * r2 * bsn *bsn * dnphi * dnphi +
119  4 * r2 * bdlog * ndlog + 4 * r * ndlog) ) /
120  2 / (1 + r * bdlog ) ;
121 
122  Cmp factor_u2 = r2 * (2 * ddbbb / bbb() - 2 * bdlog * bdlog +
123  4 * bdlog * ndlog ) +
124  2 * r2 * r2 * bsn * bsn * dnphi * dnphi +
125  4 * r * ( ndlog - bdlog ) - 6 ;
126 
127  Cmp factor_u1 = 2 * r * r2 * bsn * ( 2 * ( ndlog - bdlog ) *
128  dnphi - ddnphi ) ;
129 
130  Cmp factor_u0 = - r2 * ( 2 * ddnnn / nnn() - 2 * ndlog * ndlog +
131  4 * bdlog * ndlog ) ;
132 
133  // Scalar field the zero of which will give the marginally stable orbit
134  Cmp orbit = factor_u2 * ucor_plus * ucor_plus +
135  factor_u1 * ucor_plus + factor_u0 ;
136 
137  // Search for the zero
138  // -------------------
139 
140  int l_ms = nzet ; // index of the domain containing the MS orbit
141 
142 
143  Param par_ms ;
144  par_ms.add_int(l_ms) ;
145  par_ms.add_cmp(orbit) ;
146 
147  // Preliminary location of the zero
148  // of the orbit function with an error = 0.01
149  // The orbit closest to the star
150  double theta_ms = M_PI / 2. ; // pi/2
151  double phi_ms = 0. ;
152 
153  double xi_min = -1. ;
154  double xi_max = 1. ;
155 
156  double resloc_old ;
157  double xi_f = xi_min;
158 
159  orbit.std_base_scal() ;
160  const Valeur& vorbit = orbit.va ;
161 
162  double resloc = vorbit.val_point(l_ms, xi_min, theta_ms, phi_ms) ;
163 
164  for (int iloc=0; iloc<200; iloc++) {
165  xi_f = xi_f + 0.01 ;
166  resloc_old = resloc ;
167  resloc = vorbit.val_point(l_ms, xi_f, theta_ms, phi_ms) ;
168  if ((resloc * resloc_old) < double(0) ) {
169  xi_min = xi_f - 0.01 ;
170  xi_max = xi_f ;
171  break ;
172  }
173  }
174 
175 
176  if (ost != 0x0) {
177  *ost <<
178  "Etoile_rot::isco : preliminary location of zero of MS function :"
179  << endl ;
180  *ost << " xi_min = " << xi_min << " f(xi_min) = " <<
181  vorbit.val_point(l_ms, xi_min, theta_ms, phi_ms) << endl ;
182  *ost << " xi_max = " << xi_max << " f(xi_max) = " <<
183  vorbit.val_point(l_ms, xi_max, theta_ms, phi_ms) << endl ;
184  }
185 
186  double xi_ms = 0 ;
187  double r_ms = 0 ;
188 
189  if ( vorbit.val_point(l_ms, xi_min, theta_ms, phi_ms) *
190  vorbit.val_point(l_ms, xi_max, theta_ms, phi_ms) < double(0) ) {
191 
192 //## double precis_ms = 1.e-12 ; // precision in the determination of xi_ms
193 //## int nitermax_ms = 100 ; // max number of iterations
194 
195  int niter = 0 ;
196 
197  if (ost != 0x0) {
198  * ost <<
199  " number of iterations used in zerosec to locate the ISCO : "
200  << niter << endl ;
201  *ost << " zero found at xi = " << xi_ms << endl ;
202  }
203 
204  r_ms = mp.val_r(l_ms, xi_ms, theta_ms, phi_ms) ;
205 
206  }
207  else {
208  xi_ms = -1 ;
209  r_ms = ray_eq() ;
210  }
211 
212  p_r_isco = new double (
213  (bbb().va).val_point(l_ms, xi_ms, theta_ms, phi_ms) * r_ms
214  ) ;
215 
216  // Determination of the frequency at the marginally stable orbit
217  // -------------------------------------------------------------
218 
219  ucor_plus.std_base_scal() ;
220  double ucor_msplus = (ucor_plus.va).val_point(l_ms, xi_ms, theta_ms,
221  phi_ms) ;
222  double nobrs = (bsn.va).val_point(l_ms, xi_ms, theta_ms, phi_ms) ;
223  double nphirs = (nphi().va).val_point(l_ms, xi_ms, theta_ms, phi_ms) ;
224 
225  p_f_isco = new double ( ( ucor_msplus / nobrs / r_ms + nphirs ) /
226  (double(2) * M_PI) ) ;
227 
228  return 0 ;
229 
230 }
231 
232 
233 
234 
235 
236 
237 
238 
239 
240 }
const Cmp & dsdr() const
Returns of *this .
Definition: cmp_deriv.C:84
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:443
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 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.
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
double * p_f_isco
Orbital frequency of the ISCO.
Definition: etoile.h:1645
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
virtual double f_eccentric(double ecc, double periast, ostream *ost=0x0) const
Computation of frequency of eccentric orbits.
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:461
Coord r
r coordinate centered on the grid
Definition: map.h:718