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