LORENE
bin_bhns_extr_orbit.C
1 /*
2  * Method of class Bin_bhns_extr to compute the orbital angular velocity
3  * {\tt omega}
4  *
5  * (see file bin_bhns_extr.h for documentation).
6  *
7  */
8 
9 /*
10  * Copyright (c) 2004-2005 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_bhns_extr_orbit.C,v 1.6 2016/12/05 16:17:46 j_novak Exp $
33  * $Log: bin_bhns_extr_orbit.C,v $
34  * Revision 1.6 2016/12/05 16:17:46 j_novak
35  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
36  *
37  * Revision 1.5 2014/10/24 14:10:23 j_novak
38  * Minor change to prevent weird error from g++-4.8...
39  *
40  * Revision 1.4 2014/10/13 08:52:42 j_novak
41  * Lorene classes and functions now belong to the namespace Lorene.
42  *
43  * Revision 1.3 2014/10/06 15:13:00 j_novak
44  * Modified #include directives to use c++ syntax.
45  *
46  * Revision 1.2 2005/02/28 23:07:47 k_taniguchi
47  * Addition of the code for the conformally flat case
48  *
49  * Revision 1.1 2004/11/30 20:46:57 k_taniguchi
50  * *** empty log message ***
51  *
52  *
53  * $Header: /cvsroot/Lorene/C++/Source/Bin_bhns_extr/bin_bhns_extr_orbit.C,v 1.6 2016/12/05 16:17:46 j_novak Exp $
54  *
55  */
56 
57 // C headers
58 #include <cmath>
59 
60 // Lorene headers
61 #include "bin_bhns_extr.h"
62 #include "eos.h"
63 #include "param.h"
64 #include "utilitaires.h"
65 #include "unites.h"
66 
67 namespace Lorene {
68 double fonc_bhns_orbit_ks(double, const Param& ) ;
69 double fonc_bhns_orbit_cf(double, const Param& ) ;
70 
71 //************************************************************************
72 
73 void Bin_bhns_extr::orbit_omega_ks(double fact_omeg_min,
74  double fact_omeg_max) {
75 
76  using namespace Unites ;
77 
78  //------------------------------------------------------------------
79  // Evaluation of various quantities at the center of a neutron star
80  //------------------------------------------------------------------
81 
82  double dnulg, asn2, dasn2, nx, dnx, ny, dny, npn, dnpn ;
83  double msr ;
84 
85  const Map& mp = star.get_mp() ;
86 
87  const Cmp& logn_auto_regu = star.get_logn_auto_regu()() ;
88  const Cmp& loggam = star.get_loggam()() ;
89  const Cmp& nnn = star.get_nnn()() ;
90  const Cmp& a_car = star.get_a_car()() ;
91  const Tenseur& shift = star.get_shift() ;
92  const Tenseur& d_logn_auto_div = star.get_d_logn_auto_div() ;
93 
94  Tenseur dln_auto_div = d_logn_auto_div ;
95 
96  if ( *(dln_auto_div.get_triad()) != ref_triad ) {
97 
98  // Change the basis from spherical coordinate to Cartesian one
99  dln_auto_div.change_triad( mp.get_bvect_cart() ) ;
100 
101  // Change the basis from mapping coordinate to absolute one
102  dln_auto_div.change_triad( ref_triad ) ;
103 
104  }
105 
106  const Tenseur& d_logn_comp = star.get_d_logn_comp() ;
107 
108  Tenseur dln_comp = d_logn_comp ;
109 
110  if ( *(dln_comp.get_triad()) != ref_triad ) {
111 
112  // Change the basis from spherical coordinate to Cartesian one
113  dln_comp.change_triad( mp.get_bvect_cart() ) ;
114 
115  // Change the basis from mapping coordinate to absolute one
116  dln_comp.change_triad( ref_triad ) ;
117 
118  }
119 
120  //-------------------------------
121  // Parameters with respect to BH
122  //-------------------------------
123 
124  msr = ggrav * mass_bh / separ ;
125 
126  //-----------------------------------------------------------------
127  // Calculation of d/dX( nu + ln(Gamma) ) at the center of the star
128  // ---> dnulg
129  //-----------------------------------------------------------------
130 
131  // Factor for the coordinate transformation x --> X :
132  double factx ;
133  if (fabs(mp.get_rot_phi()) < 1.e-14) {
134  factx = 1. ;
135  }
136  else {
137  if (fabs(mp.get_rot_phi() - M_PI) < 1.e-14) {
138  factx = - 1. ;
139  }
140  else {
141  cout << "Bin_bhns_extr::orbit_omega : unknown value of rot_phi !"
142  << endl ;
143  abort() ;
144  }
145  }
146 
147  Cmp tmp = logn_auto_regu + loggam ;
148 
149  // ... gradient
150  dnulg = dln_auto_div(0)(0, 0, 0, 0) + dln_comp(0)(0, 0, 0, 0)
151  + factx * tmp.dsdx()(0, 0, 0, 0) ;
152 
153  //------------------------------------------------------------
154  // Calculation of A^2/N^2 at the center of the star ---> asn2
155  //------------------------------------------------------------
156  double nc = nnn(0, 0, 0, 0) ;
157  double a2c = a_car(0, 0, 0, 0) ;
158  asn2 = a2c / (nc * nc) ;
159 
160  if ( star.is_relativistic() ) {
161 
162  //------------------------------------------------------------------
163  // Calculation of d/dX(A^2/N^2) at the center of the star ---> dasn
164  //------------------------------------------------------------------
165 
166  double da2c = factx * a_car.dsdx()(0, 0, 0, 0) ;
167  double dnc = factx * nnn.dsdx()(0, 0, 0, 0) ;
168 
169  dasn2 = ( da2c - 2. * a2c / nc * dnc ) / (nc*nc) ;
170 
171  //------------------------------------------------------
172  // Calculation of N^X at the center of the star ---> nx
173  //------------------------------------------------------
174  nx = shift(0)(0, 0, 0, 0) ;
175 
176  //-----------------------------------------------------------
177  // Calculation of dN^X/dX at the center of the star ---> dnx
178  //-----------------------------------------------------------
179  dnx = factx * shift(0).dsdx()(0, 0, 0, 0) ;
180 
181  //------------------------------------------------------
182  // Calculation of N^Y at the center of the star ---> ny
183  //------------------------------------------------------
184  ny = shift(1)(0, 0, 0, 0) ;
185 
186  //-----------------------------------------------------------
187  // Calculation of dN^Y/dX at the center of the star ---> dny
188  //-----------------------------------------------------------
189  dny = factx * shift(1).dsdx()(0, 0, 0, 0) ;
190 
191  //--------------------------------------------
192  // Calculation of (N^X)^2 + (N^Y)^2 + (N^Z)^2
193  // at the center of the star ---> npn
194  //--------------------------------------------
195  tmp = flat_scalar_prod(shift, shift)() ; // For the flat part
196  npn = tmp(0, 0, 0, 0) ;
197 
198  //----------------------------------------------------
199  // Calculation of d/dX( (N^X)^2 + (N^Y)^2 + (N^Z)^2 )
200  // at the center of the star ---> dnpn
201  //----------------------------------------------------
202  dnpn = factx * tmp.dsdx()(0, 0, 0, 0) ;
203 
204  } // Finish of the relativistic case
205  else {
206  cout << "Bin_bhns_extr::orbit_omega : "
207  << "It should be the relativistic calculation !" << endl ;
208  abort() ;
209  }
210 
211  cout << "Bin_bhns_extr::orbit_omega: coord. ori. d(nu+log(Gam))/dX : "
212  << dnulg << endl ;
213  cout << "Bin_bhns_extr::orbit_omega: coord. ori. A^2/N^2 : "
214  << asn2 << endl ;
215  cout << "Bin_bhns_extr::orbit_omega: coord. ori. d(A^2/N^2)/dX : "
216  << dasn2 << endl ;
217  cout << "Bin_bhns_extr::orbit_omega: coord. ori. N^X : " << nx << endl ;
218  cout << "Bin_bhns_extr::orbit_omega: coord. ori. dN^X/dX : "
219  << dnx << endl ;
220  cout << "Bin_bhns_extr::orbit_omega: coord. ori. N^Y : " << ny << endl ;
221  cout << "Bin_bhns_extr::orbit_omega: coord. ori. dN^Y/dX : "
222  << dny << endl ;
223  cout << "Bin_bhns_extr::orbit_omega: coord. ori. N.N : " << npn << endl ;
224  cout << "Bin_bhns_extr::orbit_omega: coord. ori. d(N.N)/dX : "
225  << dnpn << endl ;
226 
227  //------------------------------------------------------
228  // Start of calculation of the orbital angular velocity
229  //------------------------------------------------------
230  int relat = ( star.is_relativistic() ) ? 1 : 0 ;
231 
232  double ori_x = (star.get_mp()).get_ori_x() ;
233  Param parf ;
234  parf.add_int(relat) ;
235  parf.add_double( ori_x, 0) ;
236  parf.add_double( dnulg, 1) ;
237  parf.add_double( asn2, 2) ;
238  parf.add_double( dasn2, 3) ;
239  parf.add_double( nx, 4) ;
240  parf.add_double( dnx, 5) ;
241  parf.add_double( ny, 6) ;
242  parf.add_double( dny, 7) ;
243  parf.add_double( npn, 8) ;
244  parf.add_double( dnpn, 9) ;
245  parf.add_double( msr, 10) ;
246 
247  double omega1 = fact_omeg_min * omega ;
248  double omega2 = fact_omeg_max * omega ;
249 
250  cout << "Bin_bhns_extr::orbit_omega: omega1, omega2 [rad/s] : "
251  << omega1 * f_unit << " " << omega2 * f_unit << endl ;
252 
253  // Search for the various zeros in the interval [omega1,omega2]
254  // ------------------------------------------------------------
255  int nsub = 50 ; // total number of subdivisions of the interval
256  Tbl* azer = 0x0 ;
257  Tbl* bzer = 0x0 ;
258  zero_list(fonc_bhns_orbit_ks, parf, omega1, omega2, nsub,
259  azer, bzer) ;
260 
261  // Search for the zero closest to the previous value of omega
262  // ----------------------------------------------------------
263  double omeg_min, omeg_max ;
264  int nzer = azer->get_taille() ; // number of zeros found by zero_list
265  cout << "Bin_bhns_extr::orbit_omega : " << nzer <<
266  "zero(s) found in the interval [omega1, omega2]." << endl ;
267  cout << "omega, omega1, omega2 : " << omega << " " << omega1
268  << " " << omega2 << endl ;
269  cout << "azer : " << *azer << endl ;
270  cout << "bzer : " << *bzer << endl ;
271 
272  if (nzer == 0) {
273  cout << "Bin_bhns_extr::orbit_omega: WARNING : "
274  << "no zero detected in the interval" << endl
275  << " [" << omega1 * f_unit << ", "
276  << omega2 * f_unit << "] rad/s !" << endl ;
277  omeg_min = omega1 ;
278  omeg_max = omega2 ;
279  }
280  else {
281  double dist_min = fabs(omega2 - omega1) ;
282  int i_dist_min = -1 ;
283  for (int i=0; i<nzer; i++) {
284  // Distance of previous value of omega from the center of the
285  // interval [azer(i), bzer(i)]
286  double dist = fabs( omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
287 
288  if (dist < dist_min) {
289  dist_min = dist ;
290  i_dist_min = i ;
291  }
292  }
293  omeg_min = (*azer)(i_dist_min) ;
294  omeg_max = (*bzer)(i_dist_min) ;
295  }
296 
297  delete azer ; // Tbl allocated by zero_list
298  delete bzer ; //
299 
300  cout << "Bin_bhns_extr:orbit_omega : "
301  << "interval selected for the search of the zero : "
302  << endl << " [" << omeg_min << ", " << omeg_max << "] = ["
303  << omeg_min * f_unit << ", " << omeg_max * f_unit << "] rad/s "
304  << endl ;
305 
306  // Computation of the zero in the selected interval by the secant method
307  // ---------------------------------------------------------------------
308 
309  int nitermax = 200 ;
310  int niter ;
311  double precis = 1.e-13 ;
312  omega = zerosec_b(fonc_bhns_orbit_ks, parf, omeg_min, omeg_max,
313  precis, nitermax, niter) ;
314 
315  cout << "Bin_bhns_extr::orbit_omega : "
316  << "Number of iterations in zerosec for omega : "
317  << niter << endl ;
318 
319  cout << "Bin_bhns_extr::orbit_omega : omega [rad/s] : "
320  << omega * f_unit << endl ;
321 
322 }
323 
324 
325 void Bin_bhns_extr::orbit_omega_cf(double fact_omeg_min,
326  double fact_omeg_max) {
327 
328  using namespace Unites ;
329 
330  //------------------------------------------------------------------
331  // Evaluation of various quantities at the center of a neutron star
332  //------------------------------------------------------------------
333 
334  double dnulg, asn2, dasn2, ny, dny, npn, dnpn ;
335 
336  const Map& mp = star.get_mp() ;
337 
338  const Cmp& logn_auto_regu = star.get_logn_auto_regu()() ;
339  const Cmp& logn_comp = star.get_logn_comp()() ;
340  const Cmp& loggam = star.get_loggam()() ;
341  const Cmp& nnn = star.get_nnn()() ;
342  const Cmp& a_car = star.get_a_car()() ;
343  const Tenseur& shift = star.get_shift() ;
344  const Tenseur& d_logn_auto_div = star.get_d_logn_auto_div() ;
345 
346  Tenseur dln_auto_div = d_logn_auto_div ;
347 
348  if ( *(dln_auto_div.get_triad()) != ref_triad ) {
349 
350  // Change the basis from spherical coordinate to Cartesian one
351  dln_auto_div.change_triad( mp.get_bvect_cart() ) ;
352 
353  // Change the basis from mapping coordinate to absolute one
354  dln_auto_div.change_triad( ref_triad ) ;
355 
356  }
357 
358  //-----------------------------------------------------------------
359  // Calculation of d/dX( nu + ln(Gamma) ) at the center of the star
360  // ---> dnulg
361  //-----------------------------------------------------------------
362 
363  // Factor for the coordinate transformation x --> X :
364  double factx ;
365  if (fabs(mp.get_rot_phi()) < 1.e-14) {
366  factx = 1. ;
367  }
368  else {
369  if (fabs(mp.get_rot_phi() - M_PI) < 1.e-14) {
370  factx = - 1. ;
371  }
372  else {
373  cout <<
374  "Bin_bhns_extr::orbit_omega_cfc : unknown value of rot_phi !"
375  << endl ;
376  abort() ;
377  }
378  }
379 
380  Cmp tmp = logn_auto_regu + logn_comp + loggam ;
381 
382  // ... gradient
383  dnulg = dln_auto_div(0)(0, 0, 0, 0)
384  + factx * tmp.dsdx()(0, 0, 0, 0) ;
385 
386  //------------------------------------------------------------
387  // Calculation of A^2/N^2 at the center of the star ---> asn2
388  //------------------------------------------------------------
389  double nc = nnn(0, 0, 0, 0) ;
390  double a2c = a_car(0, 0, 0, 0) ;
391  asn2 = a2c / (nc * nc) ;
392 
393  if ( star.is_relativistic() ) {
394 
395  //------------------------------------------------------------------
396  // Calculation of d/dX(A^2/N^2) at the center of the star ---> dasn
397  //------------------------------------------------------------------
398 
399  double da2c = factx * a_car.dsdx()(0, 0, 0, 0) ;
400  double dnc = factx * nnn.dsdx()(0, 0, 0, 0) ;
401 
402  dasn2 = ( da2c - 2. * a2c / nc * dnc ) / (nc*nc) ;
403 
404  //------------------------------------------------------
405  // Calculation of N^Y at the center of the star ---> ny
406  //------------------------------------------------------
407  ny = shift(1)(0, 0, 0, 0) ;
408 
409  //-----------------------------------------------------------
410  // Calculation of dN^Y/dX at the center of the star ---> dny
411  //-----------------------------------------------------------
412  dny = factx * shift(1).dsdx()(0, 0, 0, 0) ;
413 
414  //--------------------------------------------
415  // Calculation of (N^X)^2 + (N^Y)^2 + (N^Z)^2
416  // at the center of the star ---> npn
417  //--------------------------------------------
418  tmp = flat_scalar_prod(shift, shift)() ; // For the flat part
419  npn = tmp(0, 0, 0, 0) ;
420 
421  //----------------------------------------------------
422  // Calculation of d/dX( (N^X)^2 + (N^Y)^2 + (N^Z)^2 )
423  // at the center of the star ---> dnpn
424  //----------------------------------------------------
425  dnpn = factx * tmp.dsdx()(0, 0, 0, 0) ;
426 
427  } // Finish of the relativistic case
428  else {
429  cout << "Bin_bhns_extr::orbit_omega_cfc : "
430  << "It should be the relativistic calculation !" << endl ;
431  abort() ;
432  }
433 
434  cout << "Bin_bhns_extr::orbit_omega_cfc: coord. ori. d(nu+log(Gam))/dX : "
435  << dnulg << endl ;
436  cout << "Bin_bhns_extr::orbit_omega_cfc: coord. ori. A^2/N^2 : "
437  << asn2 << endl ;
438  cout << "Bin_bhns_extr::orbit_omega_cfc: coord. ori. d(A^2/N^2)/dX : "
439  << dasn2 << endl ;
440  cout << "Bin_bhns_extr::orbit_omega_cfc: coord. ori. N^Y : "
441  << ny << endl ;
442  cout << "Bin_bhns_extr::orbit_omega_cfc: coord. ori. dN^Y/dX : "
443  << dny << endl ;
444  cout << "Bin_bhns_extr::orbit_omega_cfc: coord. ori. N.N : "
445  << npn << endl ;
446  cout << "Bin_bhns_extr::orbit_omega_cfc: coord. ori. d(N.N)/dX : "
447  << dnpn << endl ;
448 
449  //------------------------------------------------------
450  // Start of calculation of the orbital angular velocity
451  //------------------------------------------------------
452  int relat = ( star.is_relativistic() ) ? 1 : 0 ;
453  double ori_x = (star.get_mp()).get_ori_x() ;
454 
455  Param parf ;
456  parf.add_int(relat) ;
457  parf.add_double( ori_x, 0) ;
458  parf.add_double( dnulg, 1) ;
459  parf.add_double( asn2, 2) ;
460  parf.add_double( dasn2, 3) ;
461  parf.add_double( ny, 4) ;
462  parf.add_double( dny, 5) ;
463  parf.add_double( npn, 6) ;
464  parf.add_double( dnpn, 7) ;
465 
466  double omega1 = fact_omeg_min * omega ;
467  double omega2 = fact_omeg_max * omega ;
468 
469  cout << "Bin_bhns_extr::orbit_omega_cfc: omega1, omega2 [rad/s] : "
470  << omega1 * f_unit << " " << omega2 * f_unit << endl ;
471 
472  // Search for the various zeros in the interval [omega1,omega2]
473  // ------------------------------------------------------------
474  int nsub = 50 ; // total number of subdivisions of the interval
475  Tbl* azer = 0x0 ;
476  Tbl* bzer = 0x0 ;
477  zero_list(fonc_bhns_orbit_cf, parf, omega1, omega2, nsub,
478  azer, bzer) ;
479 
480  // Search for the zero closest to the previous value of omega
481  // ----------------------------------------------------------
482  double omeg_min, omeg_max ;
483  int nzer = azer->get_taille() ; // number of zeros found by zero_list
484  cout << "Bin_bhns_extr::orbit_omega_cfc : " << nzer <<
485  "zero(s) found in the interval [omega1, omega2]." << endl ;
486  cout << "omega, omega1, omega2 : " << omega << " " << omega1
487  << " " << omega2 << endl ;
488  cout << "azer : " << *azer << endl ;
489  cout << "bzer : " << *bzer << endl ;
490 
491  if (nzer == 0) {
492  cout << "Bin_bhns_extr::orbit_omega_cfc: WARNING : "
493  << "no zero detected in the interval" << endl
494  << " [" << omega1 * f_unit << ", "
495  << omega2 * f_unit << "] rad/s !" << endl ;
496  omeg_min = omega1 ;
497  omeg_max = omega2 ;
498  }
499  else {
500  double dist_min = fabs(omega2 - omega1) ;
501  int i_dist_min = -1 ;
502  for (int i=0; i<nzer; i++) {
503  // Distance of previous value of omega from the center of the
504  // interval [azer(i), bzer(i)]
505  double dist = fabs( omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
506 
507  if (dist < dist_min) {
508  dist_min = dist ;
509  i_dist_min = i ;
510  }
511  }
512  omeg_min = (*azer)(i_dist_min) ;
513  omeg_max = (*bzer)(i_dist_min) ;
514  }
515 
516  delete azer ; // Tbl allocated by zero_list
517  delete bzer ; //
518 
519  cout << "Bin_bhns_extr:orbit_omega_cfc : "
520  << "interval selected for the search of the zero : "
521  << endl << " [" << omeg_min << ", " << omeg_max << "] = ["
522  << omeg_min * f_unit << ", " << omeg_max * f_unit << "] rad/s "
523  << endl ;
524 
525  // Computation of the zero in the selected interval by the secant method
526  // ---------------------------------------------------------------------
527 
528  int nitermax = 200 ;
529  int niter ;
530  double precis = 1.e-13 ;
531  omega = zerosec_b(fonc_bhns_orbit_cf, parf, omeg_min, omeg_max,
532  precis, nitermax, niter) ;
533 
534  cout << "Bin_bhns_extr::orbit_omega_cfc : "
535  << "Number of iterations in zerosec for omega : "
536  << niter << endl ;
537 
538  cout << "Bin_bhns_extr::orbit_omega_cfc : omega [rad/s] : "
539  << omega * f_unit << endl ;
540 
541 }
542 
543 
544 //***********************************************************
545 // Function used for search of the orbital angular velocity
546 //***********************************************************
547 
548 double fonc_bhns_orbit_ks(double om, const Param& parf) {
549 
550  int relat = parf.get_int() ;
551 
552  double xc = parf.get_double(0) ;
553  double dnulg = parf.get_double(1) ;
554  double asn2 = parf.get_double(2) ;
555  double dasn2 = parf.get_double(3) ;
556  double nx = parf.get_double(4) ;
557  double dnx = parf.get_double(5) ;
558  double ny = parf.get_double(6) ;
559  double dny = parf.get_double(7) ;
560  double npn = parf.get_double(8) ;
561  double dnpn = parf.get_double(9) ;
562  double msr = parf.get_double(10) ;
563 
564  double om2 = om*om ;
565 
566  double dphi_cent ;
567 
568  if (relat == 1) {
569  double bpb = om2 * xc*xc - 2.*om * ny * xc + npn + 2.*msr * nx*nx;
570 
571  dphi_cent = ( asn2* ( om* (ny + xc*dny) - om2*xc - 0.5*dnpn
572  - 2.*msr * nx * dnx + msr * nx*nx / xc )
573  - 0.5*bpb* dasn2 ) / ( 1 - asn2 * bpb ) ;
574  }
575  else {
576  cout << "Bin_bhns_extr::orbit_omega_ks : "
577  << "It should be the relativistic calculation !" << endl ;
578  abort() ;
579  }
580 
581  return dnulg + dphi_cent ;
582 
583 }
584 
585 double fonc_bhns_orbit_cf(double om, const Param& parf) {
586 
587  int relat = parf.get_int() ;
588 
589  double xc = parf.get_double(0) ;
590  double dnulg = parf.get_double(1) ;
591  double asn2 = parf.get_double(2) ;
592  double dasn2 = parf.get_double(3) ;
593  double ny = parf.get_double(4) ;
594  double dny = parf.get_double(5) ;
595  double npn = parf.get_double(6) ;
596  double dnpn = parf.get_double(7) ;
597 
598  double om2 = om*om ;
599 
600  double dphi_cent ;
601 
602  if (relat == 1) {
603  double bpb = om2 * xc*xc - 2.*om * ny * xc + npn ;
604 
605  dphi_cent = ( asn2* ( om* (ny + xc*dny) - om2*xc - 0.5*dnpn )
606  - 0.5*bpb* dasn2 ) / ( 1 - asn2 * bpb ) ;
607  }
608  else {
609  cout << "Bin_bhns_extr::orbit_omega_cf : "
610  << "It should be the relativistic calculation !" << endl ;
611  abort() ;
612  }
613 
614  return dnulg + dphi_cent ;
615 
616 }
617 }
const Tenseur & get_logn_auto_regu() const
Returns the regular part of the logarithm of the part of the lapse N generated principaly by the star...
Definition: etoile.h:711
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
const Tenseur & get_a_car() const
Returns the total conformal factor .
Definition: etoile.h:736
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
Et_bin_bhns_extr star
Neutron star.
Definition: bin_bhns_extr.h:66
const Tenseur & get_d_logn_auto_div() const
Returns the gradient of logn_auto_div.
Definition: etoile.h:722
const Tenseur & get_logn_comp() const
Returns the part of the lapse logarithm (gravitational potential at the Newtonian limit) generated pr...
Definition: etoile.h:1119
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
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition: tenseur.C:684
void orbit_omega_cf(double fact_omeg_min, double fact_omeg_max)
Computes the orbital angular velocity { omega} in the conformally flat background metric...
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tenseur.h:707
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
const Base_vect_cart ref_triad
Cartesian triad of the absolute reference frame.
Definition: bin_bhns_extr.h:63
Parameter storage.
Definition: param.h:125
const Tenseur & get_nnn() const
Returns the total lapse function N.
Definition: etoile.h:730
const Tenseur & get_d_logn_comp() const
Returns the gradient of logn_comp (Cartesian components with respect to ref_triad ) ...
Definition: etoile.h:1134
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:809
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
double separ
Absolute orbital separation between two centers of BH and NS.
Definition: bin_bhns_extr.h:74
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition: bin_bhns_extr.h:71
bool is_relativistic() const
Returns true for a relativistic star, false for a Newtonian one.
Definition: etoile.h:670
double mass_bh
Gravitational mass of BH.
Definition: bin_bhns_extr.h:77
void orbit_omega_ks(double fact_omeg_min, double fact_omeg_max)
Computes the orbital angular velocity { omega} in the Kerr-Schild background metric.
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:304