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
32char 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
116namespace Lorene {
117double fonc_binaire_axe(double , const Param& ) ;
118double fonc_binaire_orbit(double , const Param& ) ;
119
120//******************************************************************************
121
122void 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
464void 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
800double 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
857double 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}
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
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}.
const Base_vect_cart ref_triad
Cartesian triad of the absolute reference frame.
Definition binaire.h:112
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition binaire.h:129
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}.
double x_axe
Absolute X coordinate of the rotation axis.
Definition binaire.h:133
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:148
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 Tenseur & get_loggam() const
Returns the logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition etoile.h:1111
virtual double xa_barycenter() const
Absolute coordinate X of the barycenter of the baryon density, defined according to the formula.
const Tenseur & get_nnn() const
Returns the total lapse function N.
Definition etoile.h:727
const Tenseur & get_shift() const
Returns the total shift vector .
Definition etoile.h:730
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
const Map & get_mp() const
Returns the mapping.
Definition etoile.h:659
bool is_relativistic() const
Returns true for a relativistic star, false for a Newtonian one.
Definition etoile.h:667
const Tenseur & get_d_logn_auto_div() const
Returns the gradient of logn_auto_div.
Definition etoile.h:719
const Tenseur & get_a_car() const
Returns the total conformal factor .
Definition etoile.h:733
Base class for coordinate mappings.
Definition map.h:670
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
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition map.h:775
Parameter storage.
Definition param.h:125
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 int & get_int(int position=0) const
Returns the reference of a int stored in the list.
Definition param.C:292
const double & get_double(int position=0) const
Returns the reference of a double stored in the list.
Definition param.C:361
void add_int(const int &n, int position=0)
Adds the address of a new int to the list.
Definition param.C:246
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Basic array class.
Definition tbl.h:161
int get_taille() const
Gives the total size (ie dim.taille)
Definition tbl.h:397
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tenseur.h:704
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 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 .
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
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.
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
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...
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.