LORENE
binaire_orbite.C
1  /*
2  * Method of class Binaire to compute the orbital angular velocity {\tt omega}
3  * and the position of the rotation axis {\tt x_axe}.
4  *
5  * (See file binaire.h for documentation)
6  *
7  */
8 
9 /*
10  * Copyright (c) 2000-2003 Eric Gourgoulhon
11  * Copyright (c) 2000-2001 Keisuke Taniguchi
12  *
13  * This file is part of LORENE.
14  *
15  * LORENE is free software; you can redistribute it and/or modify
16  * it under the terms of the GNU General Public License as published by
17  * the Free Software Foundation; either version 2 of the License, or
18  * (at your option) any later version.
19  *
20  * LORENE is distributed in the hope that it will be useful,
21  * but WITHOUT ANY WARRANTY; without even the implied warranty of
22  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23  * GNU General Public License for more details.
24  *
25  * You should have received a copy of the GNU General Public License
26  * along with LORENE; if not, write to the Free Software
27  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
28  *
29  */
30 
31 
32 char binaire_orbite_C[] = "$Header: /cvsroot/Lorene/C++/Source/Binaire/binaire_orbite.C,v 1.8 2014/10/24 11:27:49 j_novak Exp $" ;
33 
34 /*
35  * $Id: binaire_orbite.C,v 1.8 2014/10/24 11:27:49 j_novak Exp $
36  * $Log: binaire_orbite.C,v $
37  * Revision 1.8 2014/10/24 11:27:49 j_novak
38  * Minor change in the setting of parameters for zero_list, to avoid problems with g++-4.8. To be further explored...
39  *
40  * Revision 1.7 2014/10/13 08:52:44 j_novak
41  * Lorene classes and functions now belong to the namespace Lorene.
42  *
43  * Revision 1.6 2014/10/06 15:12:59 j_novak
44  * Modified #include directives to use c++ syntax.
45  *
46  * Revision 1.5 2011/03/27 16:42:21 e_gourgoulhon
47  * Added output via new function save_profile for graphics.
48  *
49  * Revision 1.4 2009/06/18 18:40:57 k_taniguchi
50  * Added a slightly modified code to determine
51  * the orbital angular velocity.
52  *
53  * Revision 1.3 2004/03/25 10:28:59 j_novak
54  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
55  *
56  * Revision 1.2 2003/09/16 13:41:27 e_gourgoulhon
57  * Search of sub-intervals containing the zero(s) via zero_list
58  * Selection of the sub-interval as the closest one to previous value of omega.
59  * Replaced zerosec by zerosec_b.
60  *
61  * Revision 1.1.1.1 2001/11/20 15:19:30 e_gourgoulhon
62  * LORENE
63  *
64  * Revision 2.9 2001/08/14 12:36:45 keisuke
65  * Change of the minimum and maximum values for searching a new position
66  * of the rotation axis.
67  *
68  * Revision 2.8 2001/03/20 16:06:55 keisuke
69  * Addition of the method directly to calculate
70  * the orbital angular velocity (we do not use this method!).
71  *
72  * Revision 2.7 2001/03/19 17:10:28 keisuke
73  * Change of the definition of the canter of mass.
74  *
75  * Revision 2.6 2001/03/19 13:37:04 keisuke
76  * Set x_axe to be zero for the case of identical stars.
77  *
78  * Revision 2.5 2001/03/17 16:17:20 keisuke
79  * Input a subroutine to determine the position of the rotation axis
80  * in the case of binary systems composed of different stars.
81  *
82  * Revision 2.4 2000/10/02 08:29:34 keisuke
83  * Change the method to construct d_logn_auto in the calculation of dnulg[i].
84  *
85  * Revision 2.3 2000/09/22 15:56:32 keisuke
86  * *** empty log message ***
87  *
88  * Revision 2.2 2000/09/22 15:52:12 keisuke
89  * Changement calcul de dlogn (prise en compte d'une partie divergente
90  * dans logn_auto par l'appel a d_logn_auto).
91  *
92  * Revision 2.1 2000/02/12 18:36:09 eric
93  * *** empty log message ***
94  *
95  * Revision 2.0 2000/02/12 17:09:21 eric
96  * *** empty log message ***
97  *
98  *
99  * $Header: /cvsroot/Lorene/C++/Source/Binaire/binaire_orbite.C,v 1.8 2014/10/24 11:27:49 j_novak Exp $
100  *
101  */
102 
103 // Headers C
104 #include <cmath>
105 
106 // Headers Lorene
107 #include "binaire.h"
108 #include "eos.h"
109 #include "param.h"
110 #include "utilitaires.h"
111 #include "unites.h"
112 
113 #include "scalar.h"
114 #include "graphique.h"
115 
116 namespace Lorene {
117 double fonc_binaire_axe(double , const Param& ) ;
118 double fonc_binaire_orbit(double , const Param& ) ;
119 
120 //******************************************************************************
121 
122 void Binaire::orbit(double fact_omeg_min, double fact_omeg_max, double& xgg1,
123  double& xgg2) {
124 
125  using namespace Unites ;
126 
127  //-------------------------------------------------------------
128  // Evaluation of various quantities at the center of each star
129  //-------------------------------------------------------------
130 
131  double dnulg[2], asn2[2], dasn2[2], ny[2], dny[2], npn[2], dnpn[2], xgg[2] ;
132  double nyso[2], dnyso[2], npnso2[2], dnpnso2[2], ori_x[2] ;
133 
134  for (int i=0; i<2; i++) {
135 
136  const Map& mp = et[i]->get_mp() ;
137 
138  const Cmp& logn_auto_regu = et[i]->get_logn_auto_regu()() ;
139  const Cmp& logn_comp = et[i]->get_logn_comp()() ;
140  const Cmp& loggam = et[i]->get_loggam()() ;
141  const Cmp& nnn = et[i]->get_nnn()() ;
142  const Cmp& a_car = et[i]->get_a_car()() ;
143  const Tenseur& shift = et[i]->get_shift() ;
144  const Tenseur& d_logn_auto_div = et[i]->get_d_logn_auto_div() ;
145 
146  Tenseur dln_auto_div = d_logn_auto_div ;
147 
148  if ( *(dln_auto_div.get_triad()) != ref_triad ) {
149 
150  // Change the basis from spherical coordinate to Cartesian one
151  dln_auto_div.change_triad( mp.get_bvect_cart() ) ;
152 
153  // Change the basis from mapping coordinate to absolute one
154  dln_auto_div.change_triad( ref_triad ) ;
155 
156  }
157 
158  //----------------------------------
159  // Calcul de d/dX( nu + ln(Gamma) ) au centre de l'etoile ---> dnulg[i]
160  //----------------------------------
161 
162  // Facteur de passage x --> X :
163  double factx ;
164  if (fabs(mp.get_rot_phi()) < 1.e-14) {
165  factx = 1. ;
166  }
167  else {
168  if (fabs(mp.get_rot_phi() - M_PI) < 1.e-14) {
169  factx = - 1. ;
170  }
171  else {
172  cout << "Binaire::orbit : unknown value of rot_phi !" << endl ;
173  abort() ;
174  }
175  }
176 
177  Cmp tmp = logn_auto_regu + logn_comp + loggam ;
178 
179  // ... gradient suivant X :
180  dnulg[i] = dln_auto_div(0)(0, 0, 0, 0)
181  + factx * tmp.dsdx()(0, 0, 0, 0) ;
182 
183  tmp = logn_auto_regu + logn_comp ;
184  cout << "dlnndx_div_c : " << dln_auto_div(0)(0, 0, 0, 0) << endl ;
185  cout << "dlnndx_c : " << dln_auto_div(0)(0, 0, 0, 0) + factx*tmp.dsdx()(0, 0, 0, 0) << endl ;
186 
187  cout << "dloggamdx_c : " << factx*loggam.dsdx()(0, 0, 0, 0) << endl ;
188  Scalar stmp(logn_comp) ;
189  save_profile(stmp, 0., 10., 0.5*M_PI, 0., "prof_logn.d") ;
190  stmp = loggam ;
191  save_profile(stmp, 0., 1.8, 0.5*M_PI, 0., "prof_loggam.d") ;
192 
193  //----------------------------------
194  // Calcul de A^2/N^2 au centre de l'etoile ---> asn2[i]
195  //----------------------------------
196 
197  double nc = nnn(0, 0, 0, 0) ;
198  double a2c = a_car(0, 0, 0, 0) ;
199  asn2[i] = a2c / (nc * nc) ;
200 
201  if ( et[i]->is_relativistic() ) {
202 
203  //----------------------------------
204  // Calcul de d/dX(A^2/N^2) au centre de l'etoile ---> dasn2[i]
205  //----------------------------------
206 
207  double da2c = factx * a_car.dsdx()(0, 0, 0, 0) ;
208  double dnc = factx * nnn.dsdx()(0, 0, 0, 0) ;
209 
210  dasn2[i] = ( da2c - 2 * a2c / nc * dnc ) / (nc*nc) ;
211 
212  //----------------------------------
213  // Calcul de N^Y au centre de l'etoile ---> ny[i]
214  //----------------------------------
215 
216  ny[i] = shift(1)(0, 0, 0, 0) ;
217  nyso[i] = ny[i] / omega ;
218 
219  //----------------------------------
220  // Calcul de dN^Y/dX au centre de l'etoile ---> dny[i]
221  //----------------------------------
222 
223  dny[i] = factx * shift(1).dsdx()(0, 0, 0, 0) ;
224  dnyso[i] = dny[i] / omega ;
225 
226  //----------------------------------
227  // Calcul de (N^X)^2 + (N^Y)^2 + (N^Z)^2
228  // au centre de l'etoile ---> npn[i]
229  //----------------------------------
230 
231  tmp = flat_scalar_prod(shift, shift)() ;
232 
233  npn[i] = tmp(0, 0, 0, 0) ;
234  npnso2[i] = npn[i] / omega / omega ;
235 
236  //----------------------------------
237  // Calcul de d/dX( (N^X)^2 + (N^Y)^2 + (N^Z)^2 )
238  // au centre de l'etoile ---> dnpn[i]
239  //----------------------------------
240 
241  dnpn[i] = factx * tmp.dsdx()(0, 0, 0, 0) ;
242  dnpnso2[i] = dnpn[i] / omega / omega ;
243 
244  } // fin du cas relativiste
245  else {
246  dasn2[i] = 0 ;
247  ny[i] = 0 ;
248  nyso[i] = 0 ;
249  dny[i] = 0 ;
250  dnyso[i] = 0 ;
251  npn[i] = 0 ;
252  npnso2[i] = 0 ;
253  dnpn[i] = 0 ;
254  dnpnso2[i] = 0 ;
255  }
256 
257  cout << "Binaire::orbit: central d(nu+log(Gam))/dX : "
258  << dnulg[i] << endl ;
259  cout << "Binaire::orbit: central A^2/N^2 : " << asn2[i] << endl ;
260  cout << "Binaire::orbit: central d(A^2/N^2)/dX : " << dasn2[i] << endl ;
261  cout << "Binaire::orbit: central N^Y : " << ny[i] << endl ;
262  cout << "Binaire::orbit: central dN^Y/dX : " << dny[i] << endl ;
263  cout << "Binaire::orbit: central N.N : " << npn[i] << endl ;
264  cout << "Binaire::orbit: central d(N.N)/dX : " << dnpn[i] << endl ;
265 
266  //----------------------
267  // Pour information seulement : 1/ calcul des positions des "centres de
268  // de masse"
269  // 2/ calcul de dH/dX en r=0
270  //-----------------------
271 
272  ori_x[i] = (et[i]->get_mp()).get_ori_x() ;
273 
274  xgg[i] = factx * (et[i]->xa_barycenter() - ori_x[i]) ;
275 
276  } // fin de la boucle sur les etoiles
277 
278  xgg1 = xgg[0] ;
279  xgg2 = xgg[1] ;
280 
281 //---------------------------------
282 // Position de l'axe de rotation
283 //---------------------------------
284 
285  int relat = ( et[0]->is_relativistic() ) ? 1 : 0 ;
286  double ori_x1 = ori_x[0] ;
287  double ori_x2 = ori_x[1] ;
288 
289  if ( et[0]->get_eos() == et[1]->get_eos() &&
290  et[0]->get_ent()()(0,0,0,0) == et[1]->get_ent()()(0,0,0,0) ) {
291 
292  x_axe = 0. ;
293 
294  }
295  else {
296 
297  Param paraxe ;
298  paraxe.add_int(relat) ;
299  paraxe.add_double( ori_x1, 0) ;
300  paraxe.add_double( ori_x2, 1) ;
301  paraxe.add_double( dnulg[0], 2) ;
302  paraxe.add_double( dnulg[1], 3) ;
303  paraxe.add_double( asn2[0], 4) ;
304  paraxe.add_double( asn2[1], 5) ;
305  paraxe.add_double( dasn2[0], 6) ;
306  paraxe.add_double( dasn2[1], 7) ;
307  paraxe.add_double( nyso[0], 8) ;
308  paraxe.add_double( nyso[1], 9) ;
309  paraxe.add_double( dnyso[0], 10) ;
310  paraxe.add_double( dnyso[1], 11) ;
311  paraxe.add_double( npnso2[0], 12) ;
312  paraxe.add_double( npnso2[1], 13) ;
313  paraxe.add_double( dnpnso2[0], 14) ;
314  paraxe.add_double( dnpnso2[1], 15) ;
315 
316  int nitmax_axe = 200 ;
317  int nit_axe ;
318  double precis_axe = 1.e-13 ;
319 
320  x_axe = zerosec(fonc_binaire_axe, paraxe, 0.9*ori_x1, 0.9*ori_x2,
321  precis_axe, nitmax_axe, nit_axe) ;
322 
323  cout << "Binaire::orbit : Number of iterations in zerosec for x_axe : "
324  << nit_axe << endl ;
325  }
326 
327  cout << "Binaire::orbit : x_axe [km] : " << x_axe / km << endl ;
328 
329 //-------------------------------------
330 // Calcul de la vitesse orbitale
331 //-------------------------------------
332 
333  Param parf ;
334  parf.add_int(relat) ;
335  parf.add_double( ori_x1, 0) ;
336  parf.add_double( dnulg[0], 1) ;
337  parf.add_double( asn2[0], 2) ;
338  parf.add_double( dasn2[0], 3) ;
339  parf.add_double( ny[0], 4) ;
340  parf.add_double( dny[0], 5) ;
341  parf.add_double( npn[0], 6) ;
342  parf.add_double( dnpn[0], 7) ;
343  parf.add_double( x_axe, 8) ;
344 
345  double omega1 = fact_omeg_min * omega ;
346  double omega2 = fact_omeg_max * omega ;
347  cout << "Binaire::orbit: omega1, omega2 [rad/s] : "
348  << omega1 * f_unit << " " << omega2 * f_unit << endl ;
349 
350  // Search for the various zeros in the interval [omega1,omega2]
351  // ------------------------------------------------------------
352  int nsub = 50 ; // total number of subdivisions of the interval
353  Tbl* azer = 0x0 ;
354  Tbl* bzer = 0x0 ;
355  zero_list(fonc_binaire_orbit, parf, omega1, omega2, nsub,
356  azer, bzer) ;
357 
358  // Search for the zero closest to the previous value of omega
359  // ----------------------------------------------------------
360  double omeg_min, omeg_max ;
361  int nzer = azer->get_taille() ; // number of zeros found by zero_list
362  cout << "Binaire:orbit : " << nzer <<
363  " zero(s) found in the interval [omega1, omega2]." << endl ;
364  cout << "omega, omega1, omega2 : " << omega << " " << omega1
365  << " " << omega2 << endl ;
366  cout << "azer : " << *azer << endl ;
367  cout << "bzer : " << *bzer << endl ;
368 
369  if (nzer == 0) {
370  cout <<
371  "Binaire::orbit: WARNING : no zero detected in the interval"
372  << endl << " [" << omega1 * f_unit << ", "
373  << omega2 * f_unit << "] rad/s !" << endl ;
374  omeg_min = omega1 ;
375  omeg_max = omega2 ;
376  }
377  else {
378  double dist_min = fabs(omega2 - omega1) ;
379  int i_dist_min = -1 ;
380  for (int i=0; i<nzer; i++) {
381  // Distance of previous value of omega from the center of the
382  // interval [azer(i), bzer(i)]
383  double dist = fabs( omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
384  if (dist < dist_min) {
385  dist_min = dist ;
386  i_dist_min = i ;
387  }
388  }
389  omeg_min = (*azer)(i_dist_min) ;
390  omeg_max = (*bzer)(i_dist_min) ;
391  }
392 
393  delete azer ; // Tbl allocated by zero_list
394  delete bzer ; //
395 
396  cout << "Binaire:orbit : interval selected for the search of the zero : "
397  << endl << " [" << omeg_min << ", " << omeg_max << "] = ["
398  << omeg_min * f_unit << ", " << omeg_max * f_unit << "] rad/s " << endl ;
399 
400  // Computation of the zero in the selected interval by the secant method
401  // ---------------------------------------------------------------------
402  int nitermax = 200 ;
403  int niter ;
404  double precis = 1.e-13 ;
405  omega = zerosec_b(fonc_binaire_orbit, parf, omeg_min, omeg_max,
406  precis, nitermax, niter) ;
407 
408  cout << "Binaire::orbit : Number of iterations in zerosec for omega : "
409  << niter << endl ;
410 
411  cout << "Binaire::orbit : omega [rad/s] : "
412  << omega * f_unit << endl ;
413 
414 
415  /* We do not use the method below.
416  // Direct calculation
417  //--------------------
418 
419  double om2_et1 ;
420  double om2_et2 ;
421 
422  double x_et1 = ori_x1 - x_axe ;
423  double x_et2 = ori_x2 - x_axe ;
424 
425  if (relat == 1) {
426 
427  double andan_et1 = 0.5 * dasn2[0] + asn2[0] * dnulg[0] ;
428  double andan_et2 = 0.5 * dasn2[1] + asn2[1] * dnulg[1] ;
429 
430  double bpb_et1 = x_et1 * x_et1 - 2. * nyso[0] * x_et1 + npnso2[0] ;
431  double bpb_et2 = x_et2 * x_et2 - 2. * nyso[1] * x_et2 + npnso2[1] ;
432 
433  double cpc_et1 = 0.5 * dnpnso2[0] + x_et1 * (1. - dnyso[0]) - nyso[0] ;
434  double cpc_et2 = 0.5 * dnpnso2[1] + x_et2 * (1. - dnyso[1]) - nyso[1] ;
435 
436  om2_et1 = dnulg[0] / (andan_et1 * bpb_et1 + asn2[0] * cpc_et1) ;
437  om2_et2 = dnulg[1] / (andan_et2 * bpb_et2 + asn2[1] * cpc_et2) ;
438 
439  }
440  else {
441 
442  om2_et1 = dnulg[0] / x_et1 ;
443  om2_et2 = dnulg[1] / x_et2 ;
444 
445  }
446 
447  double ome_et1 = sqrt( om2_et1 ) ;
448  double ome_et2 = sqrt( om2_et2 ) ;
449 
450  double diff_om1 = 1. - ome_et1 / ome_et2 ;
451  double diff_om2 = 1. - ome_et1 / omega ;
452 
453  cout << "Binaire::orbit : omega (direct) [rad/s] : "
454  << ome_et1 * f_unit << endl ;
455 
456  cout << "Binaire::orbit : relative difference between " << endl
457  << " omega_star1 and omega_star2 (direct) : " << diff_om1 << endl
458  << " omega (direct) and omega (iretation) : " << diff_om2 << endl ;
459  */
460 
461 }
462 
463 
464 void Binaire::orbit_eqmass(double fact_omeg_min, double fact_omeg_max,
465  double mass1, double mass2,
466  double& xgg1, double& xgg2) {
467 
468  using namespace Unites ;
469 
470  //-------------------------------------------------------------
471  // Evaluation of various quantities at the center of each star
472  //-------------------------------------------------------------
473 
474  double dnulg[2], asn2[2], dasn2[2], ny[2], dny[2], npn[2], dnpn[2], xgg[2] ;
475  double nyso[2], dnyso[2], npnso2[2], dnpnso2[2], ori_x[2] ;
476 
477  for (int i=0; i<2; i++) {
478 
479  const Map& mp = et[i]->get_mp() ;
480 
481  const Cmp& logn_auto_regu = et[i]->get_logn_auto_regu()() ;
482  const Cmp& logn_comp = et[i]->get_logn_comp()() ;
483  const Cmp& loggam = et[i]->get_loggam()() ;
484  const Cmp& nnn = et[i]->get_nnn()() ;
485  const Cmp& a_car = et[i]->get_a_car()() ;
486  const Tenseur& shift = et[i]->get_shift() ;
487  const Tenseur& d_logn_auto_div = et[i]->get_d_logn_auto_div() ;
488 
489  Tenseur dln_auto_div = d_logn_auto_div ;
490 
491  if ( *(dln_auto_div.get_triad()) != ref_triad ) {
492 
493  // Change the basis from spherical coordinate to Cartesian one
494  dln_auto_div.change_triad( mp.get_bvect_cart() ) ;
495 
496  // Change the basis from mapping coordinate to absolute one
497  dln_auto_div.change_triad( ref_triad ) ;
498 
499  }
500 
501  //----------------------------------
502  // Calcul de d/dX( nu + ln(Gamma) ) au centre de l'etoile ---> dnulg[i]
503  //----------------------------------
504 
505  // Facteur de passage x --> X :
506  double factx ;
507  if (fabs(mp.get_rot_phi()) < 1.e-14) {
508  factx = 1. ;
509  }
510  else {
511  if (fabs(mp.get_rot_phi() - M_PI) < 1.e-14) {
512  factx = - 1. ;
513  }
514  else {
515  cout << "Binaire::orbit : unknown value of rot_phi !" << endl ;
516  abort() ;
517  }
518  }
519 
520  Cmp tmp = logn_auto_regu + logn_comp + loggam ;
521 
522  // ... gradient suivant X :
523  dnulg[i] = dln_auto_div(0)(0, 0, 0, 0)
524  + factx * tmp.dsdx()(0, 0, 0, 0) ;
525 
526  //----------------------------------
527  // Calcul de A^2/N^2 au centre de l'etoile ---> asn2[i]
528  //----------------------------------
529 
530  double nc = nnn(0, 0, 0, 0) ;
531  double a2c = a_car(0, 0, 0, 0) ;
532  asn2[i] = a2c / (nc * nc) ;
533 
534  if ( et[i]->is_relativistic() ) {
535 
536  //----------------------------------
537  // Calcul de d/dX(A^2/N^2) au centre de l'etoile ---> dasn2[i]
538  //----------------------------------
539 
540  double da2c = factx * a_car.dsdx()(0, 0, 0, 0) ;
541  double dnc = factx * nnn.dsdx()(0, 0, 0, 0) ;
542 
543  dasn2[i] = ( da2c - 2 * a2c / nc * dnc ) / (nc*nc) ;
544 
545  //----------------------------------
546  // Calcul de N^Y au centre de l'etoile ---> ny[i]
547  //----------------------------------
548 
549  ny[i] = shift(1)(0, 0, 0, 0) ;
550  nyso[i] = ny[i] / omega ;
551 
552  //----------------------------------
553  // Calcul de dN^Y/dX au centre de l'etoile ---> dny[i]
554  //----------------------------------
555 
556  dny[i] = factx * shift(1).dsdx()(0, 0, 0, 0) ;
557  dnyso[i] = dny[i] / omega ;
558 
559  //----------------------------------
560  // Calcul de (N^X)^2 + (N^Y)^2 + (N^Z)^2
561  // au centre de l'etoile ---> npn[i]
562  //----------------------------------
563 
564  tmp = flat_scalar_prod(shift, shift)() ;
565 
566  npn[i] = tmp(0, 0, 0, 0) ;
567  npnso2[i] = npn[i] / omega / omega ;
568 
569  //----------------------------------
570  // Calcul de d/dX( (N^X)^2 + (N^Y)^2 + (N^Z)^2 )
571  // au centre de l'etoile ---> dnpn[i]
572  //----------------------------------
573 
574  dnpn[i] = factx * tmp.dsdx()(0, 0, 0, 0) ;
575  dnpnso2[i] = dnpn[i] / omega / omega ;
576 
577  } // fin du cas relativiste
578  else {
579  dasn2[i] = 0 ;
580  ny[i] = 0 ;
581  nyso[i] = 0 ;
582  dny[i] = 0 ;
583  dnyso[i] = 0 ;
584  npn[i] = 0 ;
585  npnso2[i] = 0 ;
586  dnpn[i] = 0 ;
587  dnpnso2[i] = 0 ;
588  }
589 
590  cout << "Binaire::orbit: central d(nu+log(Gam))/dX : "
591  << dnulg[i] << endl ;
592  cout << "Binaire::orbit: central A^2/N^2 : " << asn2[i] << endl ;
593  cout << "Binaire::orbit: central d(A^2/N^2)/dX : " << dasn2[i] << endl ;
594  cout << "Binaire::orbit: central N^Y : " << ny[i] << endl ;
595  cout << "Binaire::orbit: central dN^Y/dX : " << dny[i] << endl ;
596  cout << "Binaire::orbit: central N.N : " << npn[i] << endl ;
597  cout << "Binaire::orbit: central d(N.N)/dX : " << dnpn[i] << endl ;
598 
599  //----------------------
600  // Pour information seulement : 1/ calcul des positions des "centres de
601  // de masse"
602  // 2/ calcul de dH/dX en r=0
603  //-----------------------
604 
605  ori_x[i] = (et[i]->get_mp()).get_ori_x() ;
606 
607  xgg[i] = factx * (et[i]->xa_barycenter() - ori_x[i]) ;
608 
609  } // fin de la boucle sur les etoiles
610 
611  xgg1 = xgg[0] ;
612  xgg2 = xgg[1] ;
613 
614 //---------------------------------
615 // Position de l'axe de rotation
616 //---------------------------------
617 
618  int relat = ( et[0]->is_relativistic() ) ? 1 : 0 ;
619  double ori_x1 = ori_x[0] ;
620  double ori_x2 = ori_x[1] ;
621 
622  if ( et[0]->get_eos() == et[1]->get_eos() && mass1 == mass2 ) {
623 
624  x_axe = 0. ;
625 
626  }
627  else {
628 
629  Param paraxe ;
630  paraxe.add_int(relat) ;
631  paraxe.add_double( ori_x1, 0) ;
632  paraxe.add_double( ori_x2, 1) ;
633  paraxe.add_double( dnulg[0], 2) ;
634  paraxe.add_double( dnulg[1], 3) ;
635  paraxe.add_double( asn2[0], 4) ;
636  paraxe.add_double( asn2[1], 5) ;
637  paraxe.add_double( dasn2[0], 6) ;
638  paraxe.add_double( dasn2[1], 7) ;
639  paraxe.add_double( nyso[0], 8) ;
640  paraxe.add_double( nyso[1], 9) ;
641  paraxe.add_double( dnyso[0], 10) ;
642  paraxe.add_double( dnyso[1], 11) ;
643  paraxe.add_double( npnso2[0], 12) ;
644  paraxe.add_double( npnso2[1], 13) ;
645  paraxe.add_double( dnpnso2[0], 14) ;
646  paraxe.add_double( dnpnso2[1], 15) ;
647 
648  int nitmax_axe = 200 ;
649  int nit_axe ;
650  double precis_axe = 1.e-13 ;
651 
652  x_axe = zerosec(fonc_binaire_axe, paraxe, 0.9*ori_x1, 0.9*ori_x2,
653  precis_axe, nitmax_axe, nit_axe) ;
654 
655  cout << "Binaire::orbit : Number of iterations in zerosec for x_axe : "
656  << nit_axe << endl ;
657  }
658 
659  cout << "Binaire::orbit : x_axe [km] : " << x_axe / km << endl ;
660 
661 //-------------------------------------
662 // Calcul de la vitesse orbitale
663 //-------------------------------------
664 
665  Param parf ;
666  parf.add_int(relat) ;
667  parf.add_double( ori_x1, 0) ;
668  parf.add_double( dnulg[0], 1) ;
669  parf.add_double( asn2[0], 2) ;
670  parf.add_double( dasn2[0], 3) ;
671  parf.add_double( ny[0], 4) ;
672  parf.add_double( dny[0], 5) ;
673  parf.add_double( npn[0], 6) ;
674  parf.add_double( dnpn[0], 7) ;
675  parf.add_double( x_axe, 8) ;
676 
677  double omega1 = fact_omeg_min * omega ;
678  double omega2 = fact_omeg_max * omega ;
679  cout << "Binaire::orbit: omega1, omega2 [rad/s] : "
680  << omega1 * f_unit << " " << omega2 * f_unit << endl ;
681 
682  // Search for the various zeros in the interval [omega1,omega2]
683  // ------------------------------------------------------------
684  int nsub = 50 ; // total number of subdivisions of the interval
685  Tbl* azer = 0x0 ;
686  Tbl* bzer = 0x0 ;
687  zero_list(fonc_binaire_orbit, parf, omega1, omega2, nsub,
688  azer, bzer) ;
689 
690  // Search for the zero closest to the previous value of omega
691  // ----------------------------------------------------------
692  double omeg_min, omeg_max ;
693  int nzer = azer->get_taille() ; // number of zeros found by zero_list
694  cout << "Binaire:orbit : " << nzer <<
695  " zero(s) found in the interval [omega1, omega2]." << endl ;
696  cout << "omega, omega1, omega2 : " << omega << " " << omega1
697  << " " << omega2 << endl ;
698  cout << "azer : " << *azer << endl ;
699  cout << "bzer : " << *bzer << endl ;
700 
701  if (nzer == 0) {
702  cout <<
703  "Binaire::orbit: WARNING : no zero detected in the interval"
704  << endl << " [" << omega1 * f_unit << ", "
705  << omega2 * f_unit << "] rad/s !" << endl ;
706  omeg_min = omega1 ;
707  omeg_max = omega2 ;
708  }
709  else {
710  double dist_min = fabs(omega2 - omega1) ;
711  int i_dist_min = -1 ;
712  for (int i=0; i<nzer; i++) {
713  // Distance of previous value of omega from the center of the
714  // interval [azer(i), bzer(i)]
715  double dist = fabs( omega - 0.5 * ( (*azer)(i) + (*bzer)(i) ) ) ;
716  if (dist < dist_min) {
717  dist_min = dist ;
718  i_dist_min = i ;
719  }
720  }
721  omeg_min = (*azer)(i_dist_min) ;
722  omeg_max = (*bzer)(i_dist_min) ;
723  }
724 
725  delete azer ; // Tbl allocated by zero_list
726  delete bzer ; //
727 
728  cout << "Binaire:orbit : interval selected for the search of the zero : "
729  << endl << " [" << omeg_min << ", " << omeg_max << "] = ["
730  << omeg_min * f_unit << ", " << omeg_max * f_unit << "] rad/s " << endl ;
731 
732  // Computation of the zero in the selected interval by the secant method
733  // ---------------------------------------------------------------------
734  int nitermax = 200 ;
735  int niter ;
736  double precis = 1.e-13 ;
737  omega = zerosec_b(fonc_binaire_orbit, parf, omeg_min, omeg_max,
738  precis, nitermax, niter) ;
739 
740  cout << "Binaire::orbit : Number of iterations in zerosec for omega : "
741  << niter << endl ;
742 
743  cout << "Binaire::orbit : omega [rad/s] : "
744  << omega * f_unit << endl ;
745 
746 
747  /* We do not use the method below.
748  // Direct calculation
749  //--------------------
750 
751  double om2_et1 ;
752  double om2_et2 ;
753 
754  double x_et1 = ori_x1 - x_axe ;
755  double x_et2 = ori_x2 - x_axe ;
756 
757  if (relat == 1) {
758 
759  double andan_et1 = 0.5 * dasn2[0] + asn2[0] * dnulg[0] ;
760  double andan_et2 = 0.5 * dasn2[1] + asn2[1] * dnulg[1] ;
761 
762  double bpb_et1 = x_et1 * x_et1 - 2. * nyso[0] * x_et1 + npnso2[0] ;
763  double bpb_et2 = x_et2 * x_et2 - 2. * nyso[1] * x_et2 + npnso2[1] ;
764 
765  double cpc_et1 = 0.5 * dnpnso2[0] + x_et1 * (1. - dnyso[0]) - nyso[0] ;
766  double cpc_et2 = 0.5 * dnpnso2[1] + x_et2 * (1. - dnyso[1]) - nyso[1] ;
767 
768  om2_et1 = dnulg[0] / (andan_et1 * bpb_et1 + asn2[0] * cpc_et1) ;
769  om2_et2 = dnulg[1] / (andan_et2 * bpb_et2 + asn2[1] * cpc_et2) ;
770 
771  }
772  else {
773 
774  om2_et1 = dnulg[0] / x_et1 ;
775  om2_et2 = dnulg[1] / x_et2 ;
776 
777  }
778 
779  double ome_et1 = sqrt( om2_et1 ) ;
780  double ome_et2 = sqrt( om2_et2 ) ;
781 
782  double diff_om1 = 1. - ome_et1 / ome_et2 ;
783  double diff_om2 = 1. - ome_et1 / omega ;
784 
785  cout << "Binaire::orbit : omega (direct) [rad/s] : "
786  << ome_et1 * f_unit << endl ;
787 
788  cout << "Binaire::orbit : relative difference between " << endl
789  << " omega_star1 and omega_star2 (direct) : " << diff_om1 << endl
790  << " omega (direct) and omega (iretation) : " << diff_om2 << endl ;
791  */
792 
793 }
794 
795 
796 //*************************************************
797 // Function used for search of the rotation axis
798 //*************************************************
799 
800 double fonc_binaire_axe(double x_rot, const Param& paraxe) {
801 
802  int relat = paraxe.get_int() ;
803 
804  double ori_x1 = paraxe.get_double(0) ;
805  double ori_x2 = paraxe.get_double(1) ;
806  double dnulg_1 = paraxe.get_double(2) ;
807  double dnulg_2 = paraxe.get_double(3) ;
808  double asn2_1 = paraxe.get_double(4) ;
809  double asn2_2 = paraxe.get_double(5) ;
810  double dasn2_1 = paraxe.get_double(6) ;
811  double dasn2_2 = paraxe.get_double(7) ;
812  double nyso_1 = paraxe.get_double(8) ;
813  double nyso_2 = paraxe.get_double(9) ;
814  double dnyso_1 = paraxe.get_double(10) ;
815  double dnyso_2 = paraxe.get_double(11) ;
816  double npnso2_1 = paraxe.get_double(12) ;
817  double npnso2_2 = paraxe.get_double(13) ;
818  double dnpnso2_1 = paraxe.get_double(14) ;
819  double dnpnso2_2 = paraxe.get_double(15) ;
820 
821  double om2_star1 ;
822  double om2_star2 ;
823 
824  double x1 = ori_x1 - x_rot ;
825  double x2 = ori_x2 - x_rot ;
826 
827  if (relat == 1) {
828 
829  double andan_1 = 0.5 * dasn2_1 + asn2_1 * dnulg_1 ;
830  double andan_2 = 0.5 * dasn2_2 + asn2_2 * dnulg_2 ;
831 
832  double bpb_1 = x1 * x1 - 2. * nyso_1 * x1 + npnso2_1 ;
833  double bpb_2 = x2 * x2 - 2. * nyso_2 * x2 + npnso2_2 ;
834 
835  double cpc_1 = 0.5 * dnpnso2_1 + x1 * (1. - dnyso_1) - nyso_1 ;
836  double cpc_2 = 0.5 * dnpnso2_2 + x2 * (1. - dnyso_2) - nyso_2 ;
837 
838  om2_star1 = dnulg_1 / (andan_1 * bpb_1 + asn2_1 * cpc_1) ;
839  om2_star2 = dnulg_2 / (andan_2 * bpb_2 + asn2_2 * cpc_2) ;
840 
841  }
842  else {
843 
844  om2_star1 = dnulg_1 / x1 ;
845  om2_star2 = dnulg_2 / x2 ;
846 
847  }
848 
849  return om2_star1 - om2_star2 ;
850 
851 }
852 
853 //*****************************************************************************
854 // Fonction utilisee pour la recherche de omega par la methode de la secante
855 //*****************************************************************************
856 
857 double fonc_binaire_orbit(double om, const Param& parf) {
858 
859  int relat = parf.get_int() ;
860 
861  double xc = parf.get_double(0) ;
862  double dnulg = parf.get_double(1) ;
863  double asn2 = parf.get_double(2) ;
864  double dasn2 = parf.get_double(3) ;
865  double ny = parf.get_double(4) ;
866  double dny = parf.get_double(5) ;
867  double npn = parf.get_double(6) ;
868  double dnpn = parf.get_double(7) ;
869  double x_axe = parf.get_double(8) ;
870 
871  double xx = xc - x_axe ;
872  double om2 = om*om ;
873 
874  double dphi_cent ;
875 
876  if (relat == 1) {
877  double bpb = om2 * xx*xx - 2*om * ny * xx + npn ;
878 
879  dphi_cent = ( asn2* ( om* (ny + xx*dny) - om2*xx - 0.5*dnpn )
880  - 0.5*bpb* dasn2 )
881  / ( 1 - asn2 * bpb ) ;
882  }
883  else {
884  dphi_cent = - om2*xx ;
885  }
886 
887  return dnulg + dphi_cent ;
888 
889 }
890 
891 
892 }
Lorene::Tenseur
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:301
Lorene::Binaire::et
Etoile_bin * et[2]
Array of the two stars (to perform loops on the stars): {\tt et[0]} contains the address of {\tt star...
Definition: binaire.h:124
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::Binaire::omega
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition: binaire.h:129
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::Etoile::get_eos
const Eos & get_eos() const
Returns the equation of state.
Definition: etoile.h:670
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Cmp::dsdx
const Cmp & dsdx() const
Returns of *this , where .
Definition: cmp_deriv.C:148
Lorene::Etoile::is_relativistic
bool is_relativistic() const
Returns true for a relativistic star, false for a Newtonian one.
Definition: etoile.h:667
Lorene::Tenseur::get_triad
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tenseur.h:704
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::Etoile::get_logn_auto_regu
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
Lorene::Binaire::ref_triad
const Base_vect_cart ref_triad
Cartesian triad of the absolute reference frame.
Definition: binaire.h:112
Lorene::zerosec
double zerosec(double(*f)(double, const Param &), const Param &par, double a, double b, double precis, int nitermax, int &niter, bool abort=true)
Finding the zero a function.
Definition: zerosec.C:89
Lorene::Etoile::get_shift
const Tenseur & get_shift() const
Returns the total shift vector .
Definition: etoile.h:730
Lorene::Etoile_bin::get_logn_comp
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
Unites
Standard units of space, time and mass.
Lorene::Binaire::orbit_eqmass
void orbit_eqmass(double fact_omeg_min, double fact_omeg_max, double mass1, double mass2, double &xgg1, double &xgg2)
Computes the orbital angular velocity {\tt omega} and the position of the rotation axis {\tt x_axe}.
Definition: binaire_orbite.C:464
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::Cmp
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
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::Tenseur::change_triad
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
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
Lorene::flat_scalar_prod
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...
Definition: tenseur_operateur.C:653
Lorene::Map::get_bvect_cart
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping,...
Definition: map.h:791
save_profile
void save_profile(const Scalar &uu, double r_min, double r_max, double theta, double phi, const char *filename)
Saves in a file the profile of a Scalar along some radial axis determined by a fixed value of .
Lorene::Etoile::get_d_logn_auto_div
const Tenseur & get_d_logn_auto_div() const
Returns the gradient of logn_auto_div.
Definition: etoile.h:719
Lorene::Etoile::get_ent
const Tenseur & get_ent() const
Returns the enthalpy field.
Definition: etoile.h:673
Lorene::Etoile::get_a_car
const Tenseur & get_a_car() const
Returns the total conformal factor .
Definition: etoile.h:733
Lorene::Etoile::get_mp
const Map & get_mp() const
Returns the mapping.
Definition: etoile.h:659
Lorene::Etoile_bin::get_loggam
const Tenseur & get_loggam() const
Returns the logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition: etoile.h:1111
Lorene::Etoile_bin::xa_barycenter
virtual double xa_barycenter() const
Absolute coordinate X of the barycenter of the baryon density, defined according to the formula.
Definition: et_bin_global.C:170
Lorene::Param
Parameter storage.
Definition: param.h:125
Lorene::Param::get_int
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition: param.C:292
Lorene::Etoile::get_nnn
const Tenseur & get_nnn() const
Returns the total lapse function N.
Definition: etoile.h:727
Lorene::Tbl::get_taille
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:397
Lorene::Map
Base class for coordinate mappings.
Definition: map.h:670
Lorene::Param::add_int
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition: param.C:246
Lorene::Binaire::orbit
void orbit(double fact_omeg_min, double fact_omeg_max, double &xgg1, double &xgg2)
Computes the orbital angular velocity {\tt omega} and the position of the rotation axis {\tt x_axe}.
Definition: binaire_orbite.C:122
Lorene::Binaire::x_axe
double x_axe
Absolute X coordinate of the rotation axis.
Definition: binaire.h:133