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
31char 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//=============================================================================
74namespace Lorene {
75// r_isco()
76//=============================================================================
77
78double 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}
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
Valeur va
The numerical value of the Cmp
Definition cmp.h:464
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition cmp.C:644
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:348
const Cmp & dsdr() const
Returns of *this .
Definition cmp_deriv.C:84
virtual double f_eccentric(double ecc, double periast, ostream *ost=0x0) const
Computation of frequency of eccentric orbits.
Tenseur nphi
Metric coefficient .
Definition etoile.h:1510
Tenseur bbb
Metric factor B.
Definition etoile.h:1504
double * p_f_isco
Orbital frequency of the ISCO.
Definition etoile.h:1645
double * p_r_isco
Circumferential radius of the ISCO.
Definition etoile.h:1644
int nzet
Number of domains of *mp occupied by the star.
Definition etoile.h:432
double ray_eq() const
Coordinate radius at , [r_unit].
Tenseur nnn
Total lapse function.
Definition etoile.h:509
Map & mp
Mapping associated with the star.
Definition etoile.h:429
Coord r
r coordinate centered on the grid
Definition map.h:718
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 Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
Parameter storage.
Definition param.h:125
void add_cmp(const Cmp &ti, int position=0)
Adds the address of a new Cmp to the list.
Definition param.C:935
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:246
Values and coefficients of a (real-value) function.
Definition valeur.h:287
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
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Lorene prototypes.
Definition app_hor.h:64