LORENE
bin_bhns_orbit.C
1 /*
2  * Methods of class Bin_bhns to compute the orbital angular velocity
3  *
4  * (see file bin_bhns.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2005-2007 Keisuke Taniguchi
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License version 2
15  * as published by the Free Software Foundation.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 char bin_bhns_orbit_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bin_bhns/bin_bhns_orbit.C,v 1.4 2014/10/13 08:52:41 j_novak Exp $" ;
29 
30 /*
31  * $Id: bin_bhns_orbit.C,v 1.4 2014/10/13 08:52:41 j_novak Exp $
32  * $Log: bin_bhns_orbit.C,v $
33  * Revision 1.4 2014/10/13 08:52:41 j_novak
34  * Lorene classes and functions now belong to the namespace Lorene.
35  *
36  * Revision 1.3 2014/10/06 15:13:00 j_novak
37  * Modified #include directives to use c++ syntax.
38  *
39  * Revision 1.2 2008/05/15 19:01:28 k_taniguchi
40  * Change of some parameters.
41  *
42  * Revision 1.1 2007/06/22 01:10:20 k_taniguchi
43  * *** empty log message ***
44  *
45  *
46  * $Header: /cvsroot/Lorene/C++/Source/Bin_bhns/bin_bhns_orbit.C,v 1.4 2014/10/13 08:52:41 j_novak Exp $
47  *
48  */
49 
50 // C++ headers
51 //#include <>
52 
53 // C headers
54 #include <cmath>
55 
56 // Lorene headers
57 #include "bin_bhns.h"
58 #include "param.h"
59 #include "utilitaires.h"
60 #include "unites.h"
61 
62 namespace Lorene {
63 double func_binbhns_orbit_ks(double , const Param& ) ;
64 double func_binbhns_orbit_is(double , const Param& ) ;
65 
66 //**********************************************************************
67 
68 void Bin_bhns::orbit_omega(double fact_omeg_min, double fact_omeg_max) {
69 
70  using namespace Unites ;
71 
72  if (hole.is_kerrschild()) {
73 
74  //--------------------------------------------------------------------
75  // Evaluation of various quantities at the center of the neutron star
76  //--------------------------------------------------------------------
77 
78  double dnulg, p6sl2, dp6sl2 ;
79  double shiftx, shifty, dshiftx, dshifty, shift2, dshift2 ;
80  double x_orb, y_orb, y_separ, xbh_orb, mhsr ;
81 
82  const Map& mp = star.get_mp() ;
83 
84  const Scalar& lapconf = star.get_lapconf_tot() ;
85  const Scalar& lapconf_auto = star.get_lapconf_auto() ;
86  const Scalar& confo = star.get_confo_tot() ;
87  const Scalar& confo_auto = star.get_confo_auto() ;
88  const Scalar& loggam = star.get_loggam() ;
89  const Vector& shift = star.get_shift_tot() ;
90  const Vector& shift_auto = star.get_shift_auto() ;
91 
92  const Vector& dlapconf_comp = star.get_d_lapconf_comp() ;
93  const Vector& dconfo_comp = star.get_d_confo_comp() ;
94  const Tensor& dshift_comp = star.get_d_shift_comp() ;
95 
96  const double& massbh = hole.get_mass_bh() ;
97  double mass = ggrav * massbh ;
98 
99  //----------------------------------------------------------
100  // Calculation of d/dX( ln(lapconf) - ln(psi) + ln(Gamma) )
101  // at the center of NS
102  //----------------------------------------------------------
103 
104  // Factor to translate x --> X
105  double factx ;
106  if (fabs(mp.get_rot_phi()) < 1.e-14) {
107  factx = 1. ;
108  }
109  else {
110  if (fabs(mp.get_rot_phi() - M_PI) < 1.e-14) {
111  factx = - 1. ;
112  }
113  else {
114  cout << "Bin_bhns::orbit_omega : unknown value of rot_phi !"
115  << endl ;
116  abort() ;
117  }
118  }
119 
120  Scalar tmp1(mp) ;
121  tmp1 = loggam ;
122  tmp1.std_spectral_base() ;
123 
124  // d/dX tmp1
125  dnulg = factx * ( ((lapconf_auto.dsdx()).val_grid_point(0,0,0,0)
126  + dlapconf_comp(1).val_grid_point(0,0,0,0))
127  / lapconf.val_grid_point(0,0,0,0)
128  - ((confo_auto.dsdx()).val_grid_point(0,0,0,0)
129  + dconfo_comp(1).val_grid_point(0,0,0,0))
130  / confo.val_grid_point(0,0,0,0)
131  + (tmp1.dsdx()).val_grid_point(0,0,0,0) ) ;
132 
133 
134  //----------------------------------------------------
135  // Calculation of psi^6/lapconf^2 at the center of NS
136  //----------------------------------------------------
137 
138  double lapconf_c = lapconf.val_grid_point(0,0,0,0) ;
139  double confo_c = confo.val_grid_point(0,0,0,0) ;
140 
141  p6sl2 = pow(confo_c,6.) / lapconf_c / lapconf_c ;
142 
143 
144  //----------------------------------------------------------
145  // Calculation of d/dX(psi^6/lapconf^2) at the center of NS
146  //----------------------------------------------------------
147 
148  double dlapconf_c = factx *
149  ( (lapconf_auto.dsdx()).val_grid_point(0,0,0,0)
150  + dlapconf_comp(1).val_grid_point(0,0,0,0) ) ;
151 
152  double dpsi6_c = 6. * factx * pow(confo_c,5.)
153  * ((confo_auto.dsdx()).val_grid_point(0,0,0,0)
154  + dconfo_comp(1).val_grid_point(0,0,0,0)) ;
155 
156  dp6sl2 = (dpsi6_c - 2.*pow(confo_c,6.)*dlapconf_c/lapconf_c)
157  / lapconf_c / lapconf_c ;
158 
159 
160  //--------------------------------------------------------
161  // Calculation of shift^X and shift^Y at the center of NS
162  //--------------------------------------------------------
163 
164  shiftx = shift(1).val_grid_point(0,0,0,0) ;
165  shifty = shift(2).val_grid_point(0,0,0,0) ;
166 
167 
168  //------------------------------------------------------------------
169  // Calculation of d shift^X/dX and d shift^Y/dX at the center of NS
170  //------------------------------------------------------------------
171 
172  dshiftx = factx * ((shift_auto(1).dsdx()).val_grid_point(0,0,0,0)
173  + dshift_comp(1,1).val_grid_point(0,0,0,0)) ;
174 
175  dshifty = factx * ((shift_auto(2).dsdx()).val_grid_point(0,0,0,0)
176  + dshift_comp(1,2).val_grid_point(0,0,0,0)) ;
177 
178 
179  //--------------------------------------------------------
180  // Calculation of (shift^X)^2 + (shift^Y)^2 + (shift^Z)^2
181  // at the center of NS
182  //--------------------------------------------------------
183 
184  Scalar tmp2(mp) ;
185  tmp2 = shift(1)%shift(1) + shift(2)%shift(2) + shift(3)%shift(3) ;
186  tmp2.std_spectral_base() ;
187 
188  shift2 = tmp2.val_grid_point(0,0,0,0) ;
189 
190 
191  //----------------------------------------------------------------
192  // Calculation of d/dX( (shift^X)^2 + (shift^Y)^2 + (shift^Z)^2 )
193  // at the center of NS
194  //----------------------------------------------------------------
195 
196  dshift2 = 2.*factx*((shift(1).val_grid_point(0,0,0,0))
197  * ((shift_auto(1).dsdx()).val_grid_point(0,0,0,0)
198  + dshift_comp(1,1).val_grid_point(0,0,0,0))
199  +(shift(2).val_grid_point(0,0,0,0))
200  * ((shift_auto(2).dsdx()).val_grid_point(0,0,0,0)
201  + dshift_comp(1,2).val_grid_point(0,0,0,0))
202  +(shift(3).val_grid_point(0,0,0,0))
203  * ((shift_auto(3).dsdx()).val_grid_point(0,0,0,0)
204  + dshift_comp(1,3).val_grid_point(0,0,0,0)) ) ;
205 
206 
207  //-------------------------
208  // Information of position
209  //-------------------------
210 
211  x_orb = (star.get_mp()).get_ori_x() - x_rot ;
212  y_orb = (star.get_mp()).get_ori_y() - y_rot ;
213  y_separ = (star.get_mp()).get_ori_y() - (hole.get_mp()).get_ori_y() ;
214 
215  xbh_orb = (hole.get_mp()).get_ori_x() - x_rot ;
216 
217  //------------------------------
218  // Calculation of H_BH / r_bh^4
219  //------------------------------
220 
221  mhsr = mass / pow( separ*separ+y_separ*y_separ, 2.5 ) ;
222 
223  // Output
224  // ------
225 
226  cout << "Bin_bhns::orbit_omega: central d(log(lap)+log(Gam))/dX : "
227  << dnulg << endl ;
228  cout << "Bin_bhns::orbit_omega: central psi^6/lapconf^2 : "
229  << p6sl2 << endl ;
230  cout << "Bin_bhns::orbit_omega: central d(psi^6/lapconf^2)/dX : "
231  << dp6sl2 << endl ;
232  cout << "Bin_bhns::orbit_omega: central shift^X : "
233  << shiftx << endl ;
234  cout << "Bin_bhns::orbit_omega: central shift^Y : "
235  << shifty << endl ;
236  cout << "Bin_bhns::orbit_omega: central d(shift^X)/dX : "
237  << dshiftx << endl ;
238  cout << "Bin_bhns::orbit_omega: central d(shift^Y)/dX : "
239  << dshifty << endl ;
240  cout << "Bin_bhns::orbit_omega: central shift^i shift_i : "
241  << shift2 << endl ;
242  cout << "Bin_bhns::orbit_omega: central d(shift^i shift_i)/dX : "
243  << dshift2 << endl ;
244 
245 
246  //---------------------------------------------------------------//
247  // Calculation of the orbital angular velocity //
248  //---------------------------------------------------------------//
249 
250  Param parorb ;
251  parorb.add_double(dnulg, 0) ;
252  parorb.add_double(p6sl2, 1) ;
253  parorb.add_double(dp6sl2, 2) ;
254  parorb.add_double(shiftx, 3) ;
255  parorb.add_double(shifty, 4) ;
256  parorb.add_double(dshiftx, 5) ;
257  parorb.add_double(dshifty, 6) ;
258  parorb.add_double(shift2, 7) ;
259  parorb.add_double(dshift2, 8) ;
260  parorb.add_double(x_orb, 9) ;
261  parorb.add_double(y_orb, 10) ;
262  parorb.add_double(separ, 11) ;
263  parorb.add_double(y_separ, 12) ;
264  parorb.add_double(xbh_orb, 13) ;
265  parorb.add_double(mhsr, 14) ;
266 
267  double omega1 = fact_omeg_min * omega ;
268  double omega2 = fact_omeg_max * omega ;
269 
270  cout << "Bin_bhns::orbit_omega: omega1, omega2 [rad/s] : "
271  << omega1 * f_unit << " " << omega2 * f_unit << endl ;
272 
273  // Search for the various zeros in the interval [omega1,omega2]
274  // ------------------------------------------------------------
275 
276  int nsub = 50 ; // total number of subdivisions of the interval
277  Tbl* azer = 0x0 ;
278  Tbl* bzer = 0x0 ;
279  zero_list(func_binbhns_orbit_ks, parorb, omega1, omega2, nsub,
280  azer, bzer) ;
281 
282  // Search for the zero closest to the previous value of omega
283  // ----------------------------------------------------------
284 
285  double omeg_min, omeg_max ;
286  int nzer = azer->get_taille() ; // number of zeros found by zero_list
287 
288  cout << "Bin_bhns::orbit_omega: " << nzer
289  << " zero(s) found in the interval [omega1, omega2]." << endl ;
290  cout << "omega, omega1, omega2 : " << omega << " " << omega1
291  << " " << omega2 << endl ;
292  cout << "azer : " << *azer << endl ;
293  cout << "bzer : " << *bzer << endl ;
294 
295  if (nzer == 0) {
296  cout <<
297  "Bin_bhns::orbit_omega: WARNING : no zero detected in the interval"
298  << endl << " [" << omega1 * f_unit << ", "
299  << omega2 * f_unit << "] rad/s !" << endl ;
300  omeg_min = omega1 ;
301  omeg_max = omega2 ;
302  }
303  else {
304  double dist_min = fabs(omega2 - omega1) ;
305  int i_dist_min = -1 ;
306  for (int i=0; i<nzer; i++) {
307  // Distance of previous value of omega from the center of the
308  // interval [azer(i), bzer(i)]
309 
310  double dist = fabs( omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
311 
312  if (dist < dist_min) {
313  dist_min = dist ;
314  i_dist_min = i ;
315  }
316  }
317  omeg_min = (*azer)(i_dist_min) ;
318  omeg_max = (*bzer)(i_dist_min) ;
319  }
320 
321  delete azer ; // Tbl allocated by zero_list
322  delete bzer ;
323 
324  cout << "Bin_bhns::orbit_omega: interval selected for the search of the zero : "
325  << endl << " [" << omeg_min << ", " << omeg_max << "] = ["
326  << omeg_min * f_unit << ", " << omeg_max * f_unit << "] rad/s " << endl ;
327 
328  // Computation of the zero in the selected interval by the secant method
329  // ---------------------------------------------------------------------
330 
331  int nitermax = 200 ;
332  int niter ;
333  double precis = 1.e-13 ;
334  omega = zerosec_b(func_binbhns_orbit_ks, parorb, omeg_min, omeg_max,
335  precis, nitermax, niter) ;
336 
337  cout << "Bin_bhns::orbit_omega: Number of iterations in zerosec for omega : "
338  << niter << endl ;
339 
340  cout << "Bin_bhns::orbit_omega: omega [rad/s] : "
341  << omega * f_unit << endl ;
342 
343 
344  }
345  else { // Isotropic coordinates with the maximal slicing
346 
347  //--------------------------------------------------------------------
348  // Evaluation of various quantities at the center of the neutron star
349  //--------------------------------------------------------------------
350 
351  double dnulg, p6sl2, dp6sl2 ;
352  double shiftx, shifty, dshiftx, dshifty, shift2, dshift2 ;
353  double x_orb, y_orb ;
354 
355  const Map& mp = star.get_mp() ;
356 
357  const Scalar& lapconf = star.get_lapconf_tot() ;
358  const Scalar& lapconf_auto = star.get_lapconf_auto() ;
359  const Scalar& confo = star.get_confo_tot() ;
360  const Scalar& confo_auto = star.get_confo_auto() ;
361  const Scalar& loggam = star.get_loggam() ;
362  const Vector& shift = star.get_shift_tot() ;
363  const Vector& shift_auto = star.get_shift_auto() ;
364 
365  const Vector& dlapconf_comp = star.get_d_lapconf_comp() ;
366  const Vector& dconfo_comp = star.get_d_confo_comp() ;
367  const Tensor& dshift_comp = star.get_d_shift_comp() ;
368 
369  //----------------------------------------------------------
370  // Calculation of d/dX( ln(lapconf) - ln(psi) + ln(Gamma) )
371  // at the center of NS
372  //----------------------------------------------------------
373 
374  // Factor to translate x --> X
375  double factx ;
376  if (fabs(mp.get_rot_phi()) < 1.e-14) {
377  factx = 1. ;
378  }
379  else {
380  if (fabs(mp.get_rot_phi() - M_PI) < 1.e-14) {
381  factx = - 1. ;
382  }
383  else {
384  cout << "Bin_bhns::orbit_omega: unknown value of rot_phi !"
385  << endl ;
386  abort() ;
387  }
388  }
389 
390  Scalar tmp1(mp) ;
391  tmp1 = loggam ;
392  tmp1.std_spectral_base() ;
393 
394  // d/dX tmp1
395  dnulg = factx * ( ((lapconf_auto.dsdx()).val_grid_point(0,0,0,0)
396  + dlapconf_comp(1).val_grid_point(0,0,0,0))
397  / lapconf.val_grid_point(0,0,0,0)
398  - ((confo_auto.dsdx()).val_grid_point(0,0,0,0)
399  + dconfo_comp(1).val_grid_point(0,0,0,0))
400  / confo.val_grid_point(0,0,0,0)
401  + (tmp1.dsdx()).val_grid_point(0,0,0,0) ) ;
402 
403 
404  //----------------------------------------------------
405  // Calculation of psi^6/lapconf^2 at the center of NS
406  //----------------------------------------------------
407 
408  double lapconf_c = lapconf.val_grid_point(0,0,0,0) ;
409  double confo_c = confo.val_grid_point(0,0,0,0) ;
410 
411  p6sl2 = pow(confo_c,6.) / lapconf_c / lapconf_c ;
412 
413 
414  //----------------------------------------------------------
415  // Calculation of d/dX(psi^6/lapconf^2) at the center of NS
416  //----------------------------------------------------------
417 
418  double dlapconf_c = factx *
419  ( (lapconf_auto.dsdx()).val_grid_point(0,0,0,0)
420  + dlapconf_comp(1).val_grid_point(0,0,0,0) ) ;
421 
422  double dpsi6_c = 6. * factx * pow(confo_c,5.)
423  * ((confo_auto.dsdx()).val_grid_point(0,0,0,0)
424  + dconfo_comp(1).val_grid_point(0,0,0,0)) ;
425 
426  dp6sl2 = (dpsi6_c - 2.*pow(confo_c,6.)*dlapconf_c/lapconf_c)
427  / lapconf_c / lapconf_c ;
428 
429 
430  //--------------------------------------------------------
431  // Calculation of shift^X and shift^Y at the center of NS
432  //--------------------------------------------------------
433 
434  shiftx = shift(1).val_grid_point(0,0,0,0) ;
435  shifty = shift(2).val_grid_point(0,0,0,0) ;
436 
437 
438  //------------------------------------------------------------------
439  // Calculation of d shift^X/dX and d shift^Y/dX at the center of NS
440  //------------------------------------------------------------------
441 
442  dshiftx = factx * ( (shift_auto(1).dsdx()).val_grid_point(0,0,0,0)
443  + dshift_comp(1,1).val_grid_point(0,0,0,0) ) ;
444 
445  dshifty = factx * ( (shift_auto(2).dsdx()).val_grid_point(0,0,0,0)
446  + dshift_comp(1,2).val_grid_point(0,0,0,0) ) ;
447 
448 
449  //--------------------------------------------------------
450  // Calculation of (shift^X)^2 + (shift^Y)^2 + (shift^Z)^2
451  // at the center of NS
452  //--------------------------------------------------------
453 
454  Scalar tmp2(mp) ;
455  tmp2 = shift(1)%shift(1) + shift(2)%shift(2) + shift(3)%shift(3) ;
456  tmp2.std_spectral_base() ;
457 
458  shift2 = tmp2.val_grid_point(0,0,0,0) ;
459 
460 
461  //----------------------------------------------------------------
462  // Calculation of d/dX( (shift^X)^2 + (shift^Y)^2 + (shift^Z)^2 )
463  // at the center of NS
464  //----------------------------------------------------------------
465 
466  dshift2 = 2.*factx*( (shift(1).val_grid_point(0,0,0,0))
467  * ((shift_auto(1).dsdx()).val_grid_point(0,0,0,0)
468  + dshift_comp(1,1).val_grid_point(0,0,0,0))
469  +(shift(2).val_grid_point(0,0,0,0))
470  * ((shift_auto(2).dsdx()).val_grid_point(0,0,0,0)
471  + dshift_comp(1,2).val_grid_point(0,0,0,0))
472  +(shift(3).val_grid_point(0,0,0,0))
473  * ((shift_auto(3).dsdx()).val_grid_point(0,0,0,0)
474  + dshift_comp(1,3).val_grid_point(0,0,0,0)) ) ;
475 
476 
477  //-------------------------
478  // Information of position
479  //-------------------------
480 
481  x_orb = (star.get_mp()).get_ori_x() - x_rot ;
482  y_orb = (star.get_mp()).get_ori_y() - y_rot ;
483 
484 
485  // Output
486  // ------
487 
488  cout << "Bin_bhns::orbit_omega: central d(log(lap)+log(Gam))/dX : "
489  << dnulg << endl ;
490  cout << "Bin_bhns::orbit_omega: central psi^6/lapconf^2 : "
491  << p6sl2 << endl ;
492  cout << "Bin_bhns::orbit_omega: central d(psi^6/lapconf^2)/dX : "
493  << dp6sl2 << endl ;
494  cout << "Bin_bhns::orbit_omega: central shift^X : "
495  << shiftx << endl ;
496  cout << "Bin_bhns::orbit_omega: central shift^Y : "
497  << shifty << endl ;
498  cout << "Bin_bhns::orbit_omega: central d(shift^X)/dX : "
499  << dshiftx << endl ;
500  cout << "Bin_bhns::orbit_omega: central d(shift^Y)/dX : "
501  << dshifty << endl ;
502  cout << "Bin_bhns::orbit_omega: central shift^i shift_i : "
503  << shift2 << endl ;
504  cout << "Bin_bhns::orbit_omega: central d(shift^i shift_i)/dX : "
505  << dshift2 << endl ;
506 
507 
508  //---------------------------------------------------------------//
509  // Calculation of the orbital angular velocity //
510  //---------------------------------------------------------------//
511 
512  Param parorb ;
513  parorb.add_double(dnulg, 0) ;
514  parorb.add_double(p6sl2, 1) ;
515  parorb.add_double(dp6sl2, 2) ;
516  parorb.add_double(shiftx, 3) ;
517  parorb.add_double(shifty, 4) ;
518  parorb.add_double(dshiftx, 5) ;
519  parorb.add_double(dshifty, 6) ;
520  parorb.add_double(shift2, 7) ;
521  parorb.add_double(dshift2, 8) ;
522  parorb.add_double(x_orb, 9) ;
523  parorb.add_double(y_orb, 10) ;
524 
525  double omega1 = fact_omeg_min * omega ;
526  double omega2 = fact_omeg_max * omega ;
527 
528  cout << "Bin_bhns::orbit_omega: omega1, omega2 [rad/s] : "
529  << omega1 * f_unit << " " << omega2 * f_unit << endl ;
530 
531  // Search for the various zeros in the interval [omega1,omega2]
532  // ------------------------------------------------------------
533 
534  int nsub = 50 ; // total number of subdivisions of the interval
535  Tbl* azer = 0x0 ;
536  Tbl* bzer = 0x0 ;
537  zero_list(func_binbhns_orbit_is, parorb, omega1, omega2, nsub,
538  azer, bzer) ;
539 
540  // Search for the zero closest to the previous value of omega
541  // ----------------------------------------------------------
542 
543  double omeg_min, omeg_max ;
544  int nzer = azer->get_taille() ; // number of zeros found by zero_list
545 
546  cout << "Bin_bhns::orbit_omega: " << nzer
547  << " zero(s) found in the interval [omega1, omega2]." << endl ;
548  cout << "omega, omega1, omega2 : " << omega << " " << omega1
549  << " " << omega2 << endl ;
550  cout << "azer : " << *azer << endl ;
551  cout << "bzer : " << *bzer << endl ;
552 
553  if (nzer == 0) {
554  cout <<
555  "Bin_bhns::orbit_omega: WARNING : no zero detected in the interval"
556  << endl << " [" << omega1 * f_unit << ", "
557  << omega2 * f_unit << "] rad/s !" << endl ;
558  omeg_min = omega1 ;
559  omeg_max = omega2 ;
560  }
561  else {
562  double dist_min = fabs(omega2 - omega1) ;
563  int i_dist_min = -1 ;
564  for (int i=0; i<nzer; i++) {
565  // Distance of previous value of omega from the center of the
566  // interval [azer(i), bzer(i)]
567 
568  double dist = fabs( omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
569 
570  if (dist < dist_min) {
571  dist_min = dist ;
572  i_dist_min = i ;
573  }
574  }
575  omeg_min = (*azer)(i_dist_min) ;
576  omeg_max = (*bzer)(i_dist_min) ;
577  }
578 
579  delete azer ; // Tbl allocated by zero_list
580  delete bzer ;
581 
582  cout << "Bin_bhns::orbit_omega: interval selected for the search of the zero : "
583  << endl << " [" << omeg_min << ", " << omeg_max << "] = ["
584  << omeg_min * f_unit << ", " << omeg_max * f_unit << "] rad/s " << endl ;
585 
586  // Computation of the zero in the selected interval by the secant method
587  // ---------------------------------------------------------------------
588 
589  int nitermax = 200 ;
590  int niter ;
591  double precis = 1.e-13 ;
592  omega = zerosec_b(func_binbhns_orbit_is, parorb, omeg_min, omeg_max,
593  precis, nitermax, niter) ;
594 
595  cout << "Bin_bhns::orbit_omega: Number of iterations in zerosec for omega : "
596  << niter << endl ;
597 
598  cout << "Bin_bhns::orbit_omega: omega [rad/s] : "
599  << omega * f_unit << endl ;
600 
601  }
602 }
603 
604 //******************************
605 // Function for searching omega
606 //******************************
607 
608 double func_binbhns_orbit_ks(double om, const Param& parorb) {
609 
610  double dnulg = parorb.get_double(0) ;
611  double p6sl2 = parorb.get_double(1) ;
612  double dp6sl2 = parorb.get_double(2) ;
613  double shiftx = parorb.get_double(3) ;
614  double shifty = parorb.get_double(4) ;
615  double dshiftx = parorb.get_double(5) ;
616  double dshifty = parorb.get_double(6) ;
617  double shift2 = parorb.get_double(7) ;
618  double dshift2 = parorb.get_double(8) ;
619  double x_orb = parorb.get_double(9) ;
620  double y_orb = parorb.get_double(10) ;
621  double x_separ = parorb.get_double(11) ;
622  double y_separ = parorb.get_double(12) ;
623  double xbh_orb = parorb.get_double(13) ;
624  double mhsr = parorb.get_double(14) ;
625 
626  double om2 = om * om ;
627  /*
628  double bpb = om2 * (x_orb * x_orb + y_orb * y_orb)
629  + 2. * om * (shifty * x_orb - shiftx * y_orb) + shift2
630  + 2. * mhsr * (x_separ*x_separ+y_separ*y_separ)
631  * (x_separ*shiftx + y_separ*shifty + om * y_separ * xbh_orb)
632  * (x_separ*shiftx + y_separ*shifty + om * y_separ * xbh_orb) ;
633  */
634 
635  double bpb = om2 * (x_orb * x_orb + y_orb * y_orb
636  + 2.*mhsr*(x_separ*x_separ+y_separ*y_separ)
637  * y_separ * y_separ * xbh_orb * xbh_orb)
638  + 2.*om*(shifty * x_orb - shiftx * y_orb
639  + 2.*mhsr*(x_separ*x_separ+y_separ*y_separ)
640  * (shiftx*x_separ + shifty*y_separ) * y_separ * xbh_orb)
641  + shift2 + 2.*mhsr*(x_separ*x_separ+y_separ*y_separ)
642  * (shiftx*x_separ + shifty*y_separ)
643  * (shiftx*x_separ + shifty*y_separ) ;
644 
645  double dlngam0 =
646  ( 0.5 * dp6sl2 * bpb
647  + p6sl2 * (0.5*dshift2
648  + om * (shifty - dshiftx*y_orb + dshifty*x_orb)
649  + om2 * x_orb
650  - mhsr * x_separ * (x_separ*shiftx+y_separ*shifty)
651  * (x_separ*shiftx+y_separ*shifty)
652  + 2. * mhsr * (x_separ*shiftx+y_separ*shifty)
653  * (y_separ*y_separ*shiftx - x_separ*y_separ*shifty
654  +(x_separ*x_separ+y_separ*y_separ)*(x_separ*dshiftx
655  +y_separ*dshifty))
656  + 2. * mhsr * om * y_separ * xbh_orb
657  * ( (y_separ*y_separ-2.*x_separ*x_separ)*shiftx
658  - 3. * x_separ * y_separ * shifty
659  +(x_separ*x_separ+y_separ*y_separ)*(x_separ*dshiftx
660  +y_separ*dshifty) )
661  - 3. * mhsr * om2 * x_separ * y_separ * y_separ * xbh_orb
662  * xbh_orb)
663  ) / (1 - p6sl2 * bpb) ;
664 
665  return dnulg - dlngam0 ;
666 
667 }
668 
669 double func_binbhns_orbit_is(double om, const Param& parorb) {
670 
671  double dnulg = parorb.get_double(0) ;
672  double p6sl2 = parorb.get_double(1) ;
673  double dp6sl2 = parorb.get_double(2) ;
674  double shiftx = parorb.get_double(3) ;
675  double shifty = parorb.get_double(4) ;
676  double dshiftx = parorb.get_double(5) ;
677  double dshifty = parorb.get_double(6) ;
678  double shift2 = parorb.get_double(7) ;
679  double dshift2 = parorb.get_double(8) ;
680  double x_orb = parorb.get_double(9) ;
681  double y_orb = parorb.get_double(10) ;
682 
683  double om2 = om * om ;
684 
685  double bpb = om2 * (x_orb * x_orb + y_orb * y_orb)
686  + 2. * om * (shifty * x_orb - shiftx * y_orb) + shift2 ;
687 
688  double dlngam0 = ( 0.5 * dp6sl2 * bpb
689  + p6sl2 * (0.5*dshift2
690  + om *
691  (shifty - dshiftx*y_orb + dshifty*x_orb)
692  + om2 * x_orb)
693  ) / (1 - p6sl2 * bpb) ;
694 
695  return dnulg - dlngam0 ;
696 
697 }
698 }
Lorene::Bin_bhns::hole
Hole_bhns hole
Black hole.
Definition: bin_bhns.h:72
Lorene::Star_bhns::get_loggam
const Scalar & get_loggam() const
Returns the logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition: star_bhns.h:321
Lorene::Black_hole::get_mp
const Map & get_mp() const
Returns the mapping.
Definition: blackhole.h:213
Lorene::Bin_bhns::y_rot
double y_rot
Absolute Y coordinate of the rotation axis.
Definition: bin_bhns.h:89
Lorene::Tensor
Tensor handling.
Definition: tensor.h:288
Lorene::Map::get_rot_phi
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition: map.h:775
Lorene::Param::add_double
void add_double(const double &x, int position=0)
Adds the the address of a new double to the list.
Definition: param.C:315
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Star_bhns::get_confo_tot
const Scalar & get_confo_tot() const
Returns the total conformal factor.
Definition: star_bhns.h:391
Lorene::Star_bhns::get_confo_auto
const Scalar & get_confo_auto() const
Returns the part of the conformal factor generated by the star.
Definition: star_bhns.h:383
Lorene::zero_list
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
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Lorene::Star_bhns::get_d_shift_comp
const Tensor & get_d_shift_comp() const
Returns the derivative of the shift vector generated by the companion black hole.
Definition: star_bhns.h:380
Lorene::Black_hole::get_mass_bh
double get_mass_bh() const
Returns the gravitational mass of BH [{\tt m_unit}].
Definition: blackhole.h:221
Lorene::Black_hole::is_kerrschild
bool is_kerrschild() const
Returns true for a Kerr-Schild background, false for a Conformally flat one.
Definition: blackhole.h:218
Lorene::Star_bhns::get_d_lapconf_comp
const Vector & get_d_lapconf_comp() const
Returns the derivative of the lapse function generated by the companion black hole.
Definition: star_bhns.h:361
Unites
Standard units of space, time and mass.
Lorene::Star_bhns::get_lapconf_tot
const Scalar & get_lapconf_tot() const
Returns the total lapconf function.
Definition: star_bhns.h:347
Lorene::zerosec_b
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
Lorene::Bin_bhns::separ
double separ
Absolute orbital separation between two centers of BH and NS.
Definition: bin_bhns.h:83
Lorene::Param::get_double
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition: param.C:361
Lorene::Star_bhns::get_shift_auto
const Vector & get_shift_auto() const
Returns the part of the shift vector generated by the star.
Definition: star_bhns.h:364
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
Lorene::Scalar::dsdx
const Scalar & dsdx() const
Returns of *this , where .
Definition: scalar_deriv.C:266
Lorene::Star_bhns::get_d_confo_comp
const Vector & get_d_confo_comp() const
Returns the derivative of the conformal factor generated by the companion black hole.
Definition: star_bhns.h:401
Lorene::Vector
Tensor field of valence 1.
Definition: vector.h:188
Lorene::Bin_bhns::orbit_omega
void orbit_omega(double fact_omeg_min, double fact_omeg_max)
Estimates the relative error on the Hamiltonian constraint.
Definition: bin_bhns_orbit.C:68
Lorene::Scalar::std_spectral_base
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition: scalar.C:784
Lorene::Param
Parameter storage.
Definition: param.h:125
Lorene::Bin_bhns::omega
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition: bin_bhns.h:80
Lorene::Star_bhns::get_shift_tot
const Vector & get_shift_tot() const
Returns the part total shift vector.
Definition: star_bhns.h:372
Lorene::Star::get_mp
const Map & get_mp() const
Returns the mapping.
Definition: star.h:355
Lorene::Bin_bhns::star
Star_bhns star
Neutron star.
Definition: bin_bhns.h:75
Lorene::Tbl::get_taille
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:397
Lorene::Scalar::val_grid_point
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition: scalar.h:637
Lorene::Bin_bhns::x_rot
double x_rot
Absolute X coordinate of the rotation axis.
Definition: bin_bhns.h:86
Lorene::Star_bhns::get_lapconf_auto
const Scalar & get_lapconf_auto() const
Returns the part of the lapconf function generated by the star.
Definition: star_bhns.h:339
Lorene::Map
Base class for coordinate mappings.
Definition: map.h:670