LORENE
bin_ns_bh_orbit.C
1 /*
2  * Method of class Bin_ns_bh to compute the orbital angular velocity
3  * {\tt omega}
4  *
5  * (see file bin_ns_bh.h for documentation).
6  *
7  */
8 
9 /*
10  * Copyright (c) 2003 Keisuke Taniguchi
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 version 2
16  * as published by the Free Software Foundation.
17  *
18  * LORENE is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22  *
23  * You should have received a copy of the GNU General Public License
24  * along with LORENE; if not, write to the Free Software
25  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26  *
27  */
28 
29 
30 
31 /*
32  * $Id: bin_ns_bh_orbit.C,v 1.8 2016/12/05 16:17:46 j_novak Exp $
33  * $Log: bin_ns_bh_orbit.C,v $
34  * Revision 1.8 2016/12/05 16:17:46 j_novak
35  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
36  *
37  * Revision 1.7 2014/10/24 14:10:24 j_novak
38  * Minor change to prevent weird error from g++-4.8...
39  *
40  * Revision 1.6 2014/10/13 08:52:43 j_novak
41  * Lorene classes and functions now belong to the namespace Lorene.
42  *
43  * Revision 1.5 2014/10/06 15:13:02 j_novak
44  * Modified #include directives to use c++ syntax.
45  *
46  * Revision 1.4 2004/06/09 07:26:16 k_taniguchi
47  * Minor changes.
48  *
49  * Revision 1.3 2004/06/09 06:20:11 k_taniguchi
50  * Set the standard basis for some Cmp.
51  *
52  * Revision 1.2 2004/03/25 10:28:58 j_novak
53  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
54  *
55  * Revision 1.1 2003/10/24 16:57:43 k_taniguchi
56  * Method for the calculation of the orbital angular velocity
57  *
58  *
59  * $Header: /cvsroot/Lorene/C++/Source/Bin_ns_bh/bin_ns_bh_orbit.C,v 1.8 2016/12/05 16:17:46 j_novak Exp $
60  *
61  */
62 
63 // C headers
64 #include <cmath>
65 
66 // Lorene headers
67 #include "bin_ns_bh.h"
68 #include "eos.h"
69 #include "param.h"
70 #include "utilitaires.h"
71 #include "unites.h"
72 
73 namespace Lorene {
74 double fonc_bin_ns_bh_orbit(double , const Param& ) ;
75 
76 //*************************************************************************
77 
78 void Bin_ns_bh::orbit_omega(double fact_omeg_min, double fact_omeg_max) {
79 
80  using namespace Unites ;
81 
82  //------------------------------------------------------------------
83  // Evaluation of various quantities at the center of a neutron star
84  //------------------------------------------------------------------
85 
86  double dnulg, asn2, dasn2, ny, dny, npn, dnpn ;
87 
88  const Map& mp = star.get_mp() ;
89 
90  const Cmp& loggam = star.get_loggam()() ;
91  const Cmp& nnn = star.get_nnn()() ;
92  const Cmp& confpsi = star.get_confpsi()() ;
93  const Tenseur& shift = star.get_shift() ;
94 
95  Cmp confpsi_q = pow(confpsi, 4.) ;
96  confpsi_q.std_base_scal() ;
97 
98  //-----------------------------------------------------------------
99  // Calculation of d/dX( nu + ln(Gamma) ) at the center of the star
100  // ---> dnulg
101  //-----------------------------------------------------------------
102 
103  // Factor for the coordinate transformation x --> X :
104  double factx ;
105  if (fabs(mp.get_rot_phi()) < 1.e-14) {
106  factx = 1. ;
107  }
108  else {
109  if (fabs(mp.get_rot_phi() - M_PI) < 1.e-14) {
110  factx = - 1. ;
111  }
112  else {
113  cout << "Bin_ns_bh::orbit_omega : unknown value of rot_phi !"
114  << endl ;
115  abort() ;
116  }
117  }
118 
119  Cmp tmp = log( nnn ) + loggam ;
120  tmp.std_base_scal() ;
121 
122  // ... gradient
123  dnulg = factx * tmp.dsdx()(0, 0, 0, 0) ;
124 
125  //------------------------------------------------------------
126  // Calculation of A^2/N^2 at the center of the star ---> asn2
127  //------------------------------------------------------------
128  double nc = nnn(0, 0, 0, 0) ;
129  double a2c = confpsi_q(0, 0, 0, 0) ;
130  asn2 = a2c / (nc * nc) ;
131 
132  if ( star.is_relativistic() ) {
133 
134  //------------------------------------------------------------------
135  // Calculation of d/dX(A^2/N^2) at the center of the star ---> dasn
136  //------------------------------------------------------------------
137  double da2c = factx * confpsi_q.dsdx()(0, 0, 0, 0) ;
138  double dnc = factx * nnn.dsdx()(0, 0, 0, 0) ;
139 
140  dasn2 = ( da2c - 2 * a2c / nc * dnc ) / (nc*nc) ;
141 
142  //------------------------------------------------------
143  // Calculation of N^Y at the center of the star ---> ny
144  //------------------------------------------------------
145  ny = shift(1)(0, 0, 0, 0) ;
146 
147  //-----------------------------------------------------------
148  // Calculation of dN^Y/dX at the center of the star ---> dny
149  //-----------------------------------------------------------
150  dny = factx * shift(1).dsdx()(0, 0, 0, 0) ;
151 
152  //--------------------------------------------
153  // Calculation of (N^X)^2 + (N^Y)^2 + (N^Z)^2
154  // at the center of the star ---> npn
155  //--------------------------------------------
156  tmp = flat_scalar_prod(shift, shift)() ;
157  npn = tmp(0, 0, 0, 0) ;
158 
159  //----------------------------------------------------
160  // Calculation of d/dX( (N^X)^2 + (N^Y)^2 + (N^Z)^2 )
161  // at the center of the star ---> dnpn
162  //----------------------------------------------------
163  dnpn = factx * tmp.dsdx()(0, 0, 0, 0) ;
164 
165  } // Finish of the relativistic case
166  else {
167  cout << "Bin_ns_bh::orbit_omega : "
168  << "It should be the relativistic calculation !" << endl ;
169  abort() ;
170  }
171 
172  cout << "Bin_ns_bh::orbit_omega: central d(nu+log(Gam))/dX : "
173  << dnulg << endl ;
174  cout << "Bin_ns_bh::orbit_omega: central A^2/N^2 : " << asn2 << endl ;
175  cout << "Bin_ns_bh::orbit_omega: central d(A^2/N^2)/dX : "
176  << dasn2 << endl ;
177  cout << "Bin_ns_bh::orbit_omega: central N^Y : " << ny << endl ;
178  cout << "Bin_ns_bh::orbit_omega: central dN^Y/dX : " << dny << endl ;
179  cout << "Bin_ns_bh::orbit_omega: central N.N : " << npn << endl ;
180  cout << "Bin_ns_bh::orbit_omega: central d(N.N)/dX : "
181  << dnpn << endl ;
182 
183  //------------------------------------------------------
184  // Start of calculation of the orbital angular velocity
185  //------------------------------------------------------
186  int relat = ( star.is_relativistic() ) ? 1 : 0 ;
187 
188  double ori_x = (star.get_mp()).get_ori_x() ;
189  Param parf ;
190  parf.add_int(relat) ;
191  parf.add_double( ori_x, 0) ;
192  parf.add_double( dnulg, 1) ;
193  parf.add_double( asn2, 2) ;
194  parf.add_double( dasn2, 3) ;
195  parf.add_double( ny, 4) ;
196  parf.add_double( dny, 5) ;
197  parf.add_double( npn, 6) ;
198  parf.add_double( dnpn, 7) ;
199  parf.add_double( x_axe, 8) ;
200 
201  double omega1 = fact_omeg_min * omega ;
202  double omega2 = fact_omeg_max * omega ;
203 
204  cout << "Bin_ns_bh::orbit_omega: omega1, omega2 [rad/s] : "
205  << omega1 * f_unit << " " << omega2 * f_unit << endl ;
206 
207  // Search for the various zeros in the interval [omega1,omega2]
208  // ------------------------------------------------------------
209  int nsub = 50 ; // total number of subdivisions of the interval
210  Tbl* azer = 0x0 ;
211  Tbl* bzer = 0x0 ;
212  zero_list(fonc_bin_ns_bh_orbit, parf, omega1, omega2, nsub,
213  azer, bzer) ;
214 
215  // Search for the zero closest to the previous value of omega
216  // ----------------------------------------------------------
217  double omeg_min, omeg_max ;
218  int nzer = azer->get_taille() ; // number of zeros found by zero_list
219  cout << "Bin_ns_bh:orbit_omega : " << nzer <<
220  "zero(s) found in the interval [omega1, omega2]." << endl ;
221  cout << "omega, omega1, omega2 : " << omega << " " << omega1
222  << " " << omega2 << endl ;
223  cout << "azer : " << *azer << endl ;
224  cout << "bzer : " << *bzer << endl ;
225 
226  if (nzer == 0) {
227  cout << "Bin_ns_bh::orbit_omega: WARNING : "
228  << "no zero detected in the interval" << endl
229  << " [" << omega1 * f_unit << ", "
230  << omega2 * f_unit << "] rad/s !" << endl ;
231  omeg_min = omega1 ;
232  omeg_max = omega2 ;
233  }
234  else {
235  double dist_min = fabs(omega2 - omega1) ;
236  int i_dist_min = -1 ;
237  for (int i=0; i<nzer; i++) {
238  // Distance of previous value of omega from the center of the
239  // interval [azer(i), bzer(i)]
240  double dist = fabs( omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
241 
242  if (dist < dist_min) {
243  dist_min = dist ;
244  i_dist_min = i ;
245  }
246  }
247  omeg_min = (*azer)(i_dist_min) ;
248  omeg_max = (*bzer)(i_dist_min) ;
249  }
250 
251  delete azer ; // Tbl allocated by zero_list
252  delete bzer ; //
253 
254  cout << "Bin_ns_bh:orbit_omega : "
255  << "interval selected for the search of the zero : "
256  << endl << " [" << omeg_min << ", " << omeg_max << "] = ["
257  << omeg_min * f_unit << ", " << omeg_max * f_unit << "] rad/s "
258  << endl ;
259 
260  // Computation of the zero in the selected interval by the secant method
261  // ---------------------------------------------------------------------
262 
263  int nitermax = 200 ;
264  int niter ;
265  double precis = 1.e-13 ;
266  omega = zerosec_b(fonc_bin_ns_bh_orbit, parf, omeg_min, omeg_max,
267  precis, nitermax, niter) ;
268 
269  cout << "Bin_ns_bh::orbit_omega : "
270  << "Number of iterations in zerosec for omega : "
271  << niter << endl ;
272 
273  cout << "Bin_ns_bh::orbit_omega : omega [rad/s] : "
274  << omega * f_unit << endl ;
275 
276 }
277 
278 //***********************************************************
279 // Function used for search of the orbital angular velocity
280 //***********************************************************
281 
282 double fonc_bin_ns_bh_orbit(double om, const Param& parf) {
283 
284  int relat = parf.get_int() ;
285 
286  double xc = parf.get_double(0) ;
287  double dnulg = parf.get_double(1) ;
288  double asn2 = parf.get_double(2) ;
289  double dasn2 = parf.get_double(3) ;
290  double ny = parf.get_double(4) ;
291  double dny = parf.get_double(5) ;
292  double npn = parf.get_double(6) ;
293  double dnpn = parf.get_double(7) ;
294  double x_axe = parf.get_double(8) ;
295 
296  double xx = xc - x_axe ;
297  double om2 = om*om ;
298 
299  double dphi_cent ;
300 
301  if (relat == 1) {
302  double bpb = om2 * xx*xx - 2*om * ny * xx + npn ;
303 
304  dphi_cent = ( asn2* ( om* (ny + xx*dny) - om2*xx - 0.5*dnpn )
305  - 0.5*bpb* dasn2 )
306  / ( 1 - asn2 * bpb ) ;
307  }
308  else {
309  cout << "Bin_ns_bh::orbit_omega : "
310  << "It should be the relativistic calculation !" << endl ;
311  abort() ;
312  }
313 
314  return dnulg + dphi_cent ;
315 
316 }
317 }
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
const Cmp & dsdx() const
Returns of *this , where .
Definition: cmp_deriv.C:151
const Map & get_mp() const
Returns the mapping.
Definition: etoile.h:662
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:249
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
const Tenseur & get_shift() const
Returns the total shift vector .
Definition: etoile.h:733
Base class for coordinate mappings.
Definition: map.h:688
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
double zerosec_b(double(*f)(double, const Param &), const Param &par, double a, double b, double precis, int nitermax, int &niter)
Finding the zero a function on a bounded domain.
Definition: zerosec_borne.C:71
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition: map.h:793
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:295
const Tenseur & get_loggam() const
Returns the logarithm of the Lorentz factor between the fluid and the co-orbiting observer...
Definition: etoile.h:1114
double x_axe
Absolute X coordinate of the rotation axis.
Definition: bin_ns_bh.h:143
Et_bin_nsbh star
The neutron star.
Definition: bin_ns_bh.h:131
void zero_list(double(*f)(double, const Param &), const Param &par, double xmin, double xmax, int nsub, Tbl *&az, Tbl *&bz)
Locates approximatively all the zeros of a function in a given interval.
Definition: zero_list.C:60
Parameter storage.
Definition: param.h:125
const Tenseur & get_nnn() const
Returns the total lapse function N.
Definition: etoile.h:730
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition: bin_ns_bh.h:139
const Tenseur & get_confpsi() const
Returns the part of the conformal factor $$.
Definition: et_bin_nsbh.h:257
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
void orbit_omega(double fact_omeg_min, double fact_omeg_max)
Computes the orbital angular velocity { omega}.
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:647
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:417
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:318
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition: param.C:364
Basic array class.
Definition: tbl.h:164
bool is_relativistic() const
Returns true for a relativistic star, false for a Newtonian one.
Definition: etoile.h:670
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:304