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