LORENE
etoile_bin.C
1/*
2 * Methods for the class Etoile_bin
3 *
4 * (see file etoile.h for documentation)
5 */
6
7/*
8 * Copyright (c) 2000-2001 Eric Gourgoulhon
9 * Copyright (c) 2000-2001 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 as published by
15 * the Free Software Foundation; either version 2 of the License, or
16 * (at your option) any later version.
17 *
18 * LORENE is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU General Public License for more details.
22 *
23 * You should have received a copy of the GNU General Public License
24 * along with LORENE; if not, write to the Free Software
25 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26 *
27 */
28
29
30char etoile_bin_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/etoile_bin.C,v 1.13 2014/10/13 08:52:58 j_novak Exp $" ;
31
32/*
33 * $Id: etoile_bin.C,v 1.13 2014/10/13 08:52:58 j_novak Exp $
34 * $Log: etoile_bin.C,v $
35 * Revision 1.13 2014/10/13 08:52:58 j_novak
36 * Lorene classes and functions now belong to the namespace Lorene.
37 *
38 * Revision 1.12 2004/11/25 09:53:55 m_bejger
39 * Corrected error in fait_d_psi() in the case where nzet > 1.
40 *
41 * Revision 1.11 2004/05/07 12:35:16 f_limousin
42 * Add new member ssjm1_psi
43 *
44 * Revision 1.10 2004/03/25 10:29:06 j_novak
45 * All LORENE's units are now defined in the namespace Unites (in file unites.h).
46 *
47 * Revision 1.9 2003/10/21 11:47:56 k_taniguchi
48 * Delete various things for the Bin_ns_bh project.
49 * They are moved to et_bin_nsbh.C.
50 *
51 * Revision 1.8 2003/10/20 15:08:22 k_taniguchi
52 * Minor changes.
53 *
54 * Revision 1.7 2003/10/20 14:51:26 k_taniguchi
55 * Addition of various things for the Bin_ns_bh project
56 * which are related with the part of the neutron star.
57 *
58 * Revision 1.6 2003/02/06 16:08:36 f_limousin
59 * Modified Etoile_bin::sprod in order to avoid a warning from the compiler
60 *
61 * Revision 1.5 2003/02/03 12:52:15 f_limousin
62 * *** empty log message ***
63 *
64 * Revision 1.4 2003/01/31 16:57:12 p_grandclement
65 * addition of the member Cmp decouple used to compute the K_ij auto, once
66 * the K_ij total is known
67 *
68 * Revision 1.3 2003/01/17 13:39:51 f_limousin
69 * Modification of sprod routine
70 *
71 * Revision 1.2 2002/12/17 21:20:29 e_gourgoulhon
72 * Suppression of the member p_companion,
73 * as well as the associated function set_companion.
74 *
75 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
76 * LORENE
77 *
78 * Revision 2.31 2001/06/25 12:53:15 eric
79 * Ajout du membre p_companion et des fonctions associees set_companion() et get_companion().
80 *
81 * Revision 2.30 2000/12/22 13:09:09 eric
82 * Modif fait_d_psi : prolongement C^1 de dpsi en dehors de l'etoile.
83 *
84 * Revision 2.29 2000/09/29 11:54:40 keisuke
85 * Add the relaxations for logn_auto_div and d_logn_auto_div.
86 *
87 * Revision 2.28 2000/09/29 09:57:26 keisuke
88 * Add the relaxation for logn_auto_regu.
89 *
90 * Revision 2.27 2000/09/22 15:51:19 keisuke
91 * d_logn_auto_div devient desormais un membre de la classe Etoile
92 * et non plus de la classe derivee Etoile_bin.
93 *
94 * Revision 2.26 2000/09/07 14:35:31 keisuke
95 * Ajout du membre d_logn_auto_regu.
96 *
97 * Revision 2.25 2000/08/29 11:38:03 eric
98 * Ajout du membre d_logn_auto_div.
99 *
100 * Revision 2.24 2000/07/06 09:53:22 eric
101 * Ajout de xa_barycenter dans l'affichage.
102 *
103 * Revision 2.23 2000/07/06 09:40:11 eric
104 * Ajout du membre derive p_xa_barycenter.
105 *
106 * Revision 2.22 2000/06/15 15:50:02 eric
107 * Ajout du calcul de d_tilde dans l'affichage.
108 *
109 * Revision 2.21 2000/05/23 12:28:00 phil
110 * changement apres modification de skxk
111 *
112 * Revision 2.20 2000/03/15 11:04:58 eric
113 * Ajout des fonctions Etoile_bin::set_w_shift() et Etoile_bin::set_khi_shift()
114 *
115 * Revision 2.19 2000/02/24 09:12:56 keisuke
116 * Add an output for the velocity field in the corotating frame.
117 *
118 * Revision 2.18 2000/02/21 16:31:58 eric
119 * Modif affichage.
120 *
121 * Revision 2.17 2000/02/21 14:54:13 eric
122 * fait_d_psi appel d_psi.set_etat_nondef() et return dans le cas
123 * corotation.
124 *
125 * Revision 2.16 2000/02/21 14:31:43 eric
126 * gam_euler est desormais sauve dans le fichier (cas irrotationnel seulement)
127 * psi0 n'est sauve que dans le cas irrotationnel.
128 *
129 * Revision 2.15 2000/02/21 13:58:22 eric
130 * Suppression du membre psi: remplacement par psi0.
131 *
132 * Revision 2.14 2000/02/17 15:30:04 eric
133 * Ajout de la fonction Etoile_bin::relaxation.
134 *
135 * Revision 2.13 2000/02/17 14:42:09 eric
136 * Modif affichage.
137 *
138 * Revision 2.12 2000/02/16 17:12:25 eric
139 * Modif initialisation de w_shift, khi_shift et ssjm1_wshift dans le
140 * constructeur standard.
141 *
142 * Revision 2.11 2000/02/16 16:08:10 eric
143 * w_shift et ssjm1_wshift sont desormais definis sur la triade du mapping.
144 *
145 * Revision 2.10 2000/02/16 15:06:11 eric
146 * Ajout des membres w_shift et khi_shift.
147 * (sauves dans les fichiers a la place de shift_auto).
148 * Ajout de la fonction Etoile_bin::fait_shift_auto.
149 *
150 * Revision 2.9 2000/02/16 13:47:22 eric
151 * Ajout des membres ssjm1_khi et ssjm1_wshift.
152 * (sauvegardes dans les fichiers).
153 *
154 * Revision 2.8 2000/02/16 11:54:40 eric
155 * Ajout des membres ssjm1_logn et ssjm1_beta (sauvegarde dans les fichiers).
156 *
157 * Revision 2.7 2000/02/10 18:56:52 eric
158 * Modifs routine d'affichage.
159 *
160 * Revision 2.6 2000/02/10 16:12:05 eric
161 * Ajout de la fonction fait_d_psi
162 *
163 * Revision 2.5 2000/02/09 20:24:03 eric
164 * Appel de set_triad(ref_triad) sur u_euler et shift dans les
165 * constructeurs.
166 * ,.
167 *
168 * Revision 2.4 2000/02/09 19:31:33 eric
169 * La triade de decomposition doit desormais figurer en argument des
170 * constructeurs de Tenseur.
171 *
172 * Revision 2.3 2000/02/08 19:29:50 eric
173 * La fonction Etoile_bin::scal_prod est rebaptisee Etoile_bin::sprod
174 *
175 * Revision 2.2 2000/02/04 17:15:36 eric
176 * Ajout du membre ref_triad.
177 *
178 * Revision 2.1 2000/02/01 16:00:00 eric
179 * Ajout de la fonction Etoile_bin::scal_prod.
180 *
181 * Revision 2.0 2000/01/31 15:57:12 eric
182 * *** empty log message ***
183 *
184 *
185 * $Header: /cvsroot/Lorene/C++/Source/Etoile/etoile_bin.C,v 1.13 2014/10/13 08:52:58 j_novak Exp $
186 *
187 */
188
189// Headers C
190#include "math.h"
191
192// Headers Lorene
193#include "etoile.h"
194#include "eos.h"
195#include "unites.h"
196
197// Local prototype
198namespace Lorene {
199Cmp raccord_c1(const Cmp& uu, int l1) ;
200
201 //--------------//
202 // Constructors //
203 //--------------//
204
205// Standard constructor
206// --------------------
207Etoile_bin::Etoile_bin(Map& mpi, int nzet_i, bool relat, const Eos& eos_i,
208 bool irrot, const Base_vect& ref_triad_i)
209 : Etoile(mpi, nzet_i, relat, eos_i),
210 irrotational(irrot),
211 ref_triad(ref_triad_i),
212 psi0(mpi),
213 d_psi(mpi, 1, COV, ref_triad),
214 wit_w(mpi, 1, CON, ref_triad),
215 loggam(mpi),
216 logn_comp(mpi),
217 d_logn_auto(mpi, 1, COV, ref_triad),
218 d_logn_auto_regu(mpi, 1, COV, ref_triad),
219 d_logn_comp(mpi, 1, COV, ref_triad),
220 beta_comp(mpi),
221 d_beta_auto(mpi, 1, COV, ref_triad),
222 d_beta_comp(mpi, 1, COV, ref_triad),
223 shift_auto(mpi, 1, CON, ref_triad),
224 shift_comp(mpi, 1, CON, ref_triad),
225 w_shift(mpi, 1, CON, mp.get_bvect_cart()),
226 khi_shift(mpi),
227 tkij_auto(mpi, 2, CON, ref_triad),
228 tkij_comp(mpi, 2, CON, ref_triad),
229 akcar_auto(mpi),
230 akcar_comp(mpi),
231 bsn(mpi, 1, CON, ref_triad),
232 pot_centri(mpi),
233 ssjm1_logn(mpi),
234 ssjm1_beta(mpi),
235 ssjm1_khi(mpi),
236 ssjm1_wshift(mpi, 1, CON, mp.get_bvect_cart()),
237 ssjm1_psi(mpi),
238 decouple(mpi)
239{
240
241 // Pointers of derived quantities initialized to zero :
242 set_der_0x0() ;
243
244 // The reference triad is assigned to the vectors u_euler and shift :
247
248 // All quantities are initialized to zero :
249 psi0 = 0 ;
250 d_psi = 0 ;
251 wit_w = 0 ;
252 loggam = 0 ;
253 logn_comp = 0 ;
254 d_logn_auto = 0 ;
255 d_logn_auto_regu = 0 ;
256 d_logn_comp = 0 ;
257 beta_comp = 0 ;
258 d_beta_auto = 0 ;
259 d_beta_comp = 0 ;
260 shift_auto = 0 ;
261 shift_comp = 0 ;
262
264 for (int i=0; i<3; i++) {
265 w_shift.set(i) = 0 ;
266 }
267
269 khi_shift.set() = 0 ;
270
273 akcar_auto = 0 ;
274 akcar_comp = 0 ;
275 bsn = 0 ;
276 pot_centri = 0 ;
277 ssjm1_logn = 0 ;
278 ssjm1_beta = 0 ;
279 ssjm1_khi = 0 ;
280 ssjm1_psi = 0 ;
281
283 for (int i=0; i<3; i++) {
284 ssjm1_wshift.set(i) = 0 ;
285 }
286
287}
288
289// Copy constructor
290// ----------------
292 : Etoile(et),
293 irrotational(et.irrotational),
294 ref_triad(et.ref_triad),
295 psi0(et.psi0),
296 d_psi(et.d_psi),
297 wit_w(et.wit_w),
298 loggam(et.loggam),
299 logn_comp(et.logn_comp),
300 d_logn_auto(et.d_logn_auto),
301 d_logn_auto_regu(et.d_logn_auto_regu),
302 d_logn_comp(et.d_logn_comp),
303 beta_comp(et.beta_comp),
304 d_beta_auto(et.d_beta_auto),
305 d_beta_comp(et.d_beta_comp),
306 shift_auto(et.shift_auto),
307 shift_comp(et.shift_comp),
308 w_shift(et.w_shift),
309 khi_shift(et.khi_shift),
310 tkij_auto(et.tkij_auto),
311 tkij_comp(et.tkij_comp),
312 akcar_auto(et.akcar_auto),
313 akcar_comp(et.akcar_comp),
314 bsn(et.bsn),
315 pot_centri(et.pot_centri),
316 ssjm1_logn(et.ssjm1_logn),
317 ssjm1_beta(et.ssjm1_beta),
318 ssjm1_khi(et.ssjm1_khi),
319 ssjm1_wshift(et.ssjm1_wshift),
320 ssjm1_psi(et.ssjm1_khi),
321 decouple(et.decouple)
322{
323 set_der_0x0() ;
324
325}
326
327// Constructor from a file
328// -----------------------
329Etoile_bin::Etoile_bin(Map& mpi, const Eos& eos_i, const Base_vect& ref_triad_i,
330 FILE* fich)
331 : Etoile(mpi, eos_i, fich),
332 ref_triad(ref_triad_i),
333 psi0(mpi),
334 d_psi(mpi, 1, COV, ref_triad),
335 wit_w(mpi, 1, CON, ref_triad),
336 loggam(mpi),
337 logn_comp(mpi),
338 d_logn_auto(mpi, 1, COV, ref_triad),
339 d_logn_auto_regu(mpi, 1, COV, ref_triad),
340 d_logn_comp(mpi, 1, COV, ref_triad),
341 beta_comp(mpi),
342 d_beta_auto(mpi, 1, COV, ref_triad),
343 d_beta_comp(mpi, 1, COV, ref_triad),
344 shift_auto(mpi, 1, CON, ref_triad),
345 shift_comp(mpi, 1, CON, ref_triad),
346 w_shift(mpi, 1, CON, mp.get_bvect_cart()),
347 khi_shift(mpi),
348 tkij_auto(mpi, 2, CON, ref_triad),
349 tkij_comp(mpi, 2, CON, ref_triad),
350 akcar_auto(mpi),
351 akcar_comp(mpi),
352 bsn(mpi, 1, CON, ref_triad),
353 pot_centri(mpi),
354 ssjm1_logn(mpi),
355 ssjm1_beta(mpi),
356 ssjm1_khi(mpi),
357 ssjm1_wshift(mpi, 1, CON, mp.get_bvect_cart()),
358 ssjm1_psi(mpi),
359 decouple(mpi)
360{
361
362 // The reference triad is assigned to the vectors u_euler and shift :
365
366 // Etoile parameters
367 // -----------------
368
369 // irrotational is read in the file:
370 fread(&irrotational, sizeof(bool), 1, fich) ;
371
372
373 // Read of the saved fields:
374 // ------------------------
375
376 if (irrotational) {
377 Tenseur gam_euler_file(mp, fich) ;
378 gam_euler = gam_euler_file ;
379
380 Tenseur psi0_file(mp, fich) ;
381 psi0 = psi0_file ;
382 }
383
384
385 Tenseur w_shift_file(mp, mp.get_bvect_cart(), fich) ;
386 w_shift = w_shift_file ;
387
388 Tenseur khi_shift_file(mp, fich) ;
389 khi_shift = khi_shift_file ;
390
391 fait_shift_auto() ; // constructs shift_auto from w_shift and khi_shift
392
393 Cmp ssjm1_logn_file(mp, *(mp.get_mg()), fich) ;
394 ssjm1_logn = ssjm1_logn_file ;
395
396 Cmp ssjm1_beta_file(mp, *(mp.get_mg()), fich) ;
397 ssjm1_beta = ssjm1_beta_file ;
398
399 Cmp ssjm1_khi_file(mp, *(mp.get_mg()), fich) ;
400 ssjm1_khi = ssjm1_khi_file ;
401
402 Tenseur ssjm1_wshift_file(mp, mp.get_bvect_cart(), fich) ;
403 ssjm1_wshift = ssjm1_wshift_file ;
404
405 ssjm1_psi = 0 ;
406
407 // All other fields are initialized to zero :
408 // ----------------------------------------
409 d_psi = 0 ;
410 wit_w = 0 ;
411 loggam = 0 ;
412 logn_comp = 0 ;
413 d_logn_auto = 0 ;
414 d_logn_auto_regu = 0 ;
415 d_logn_comp = 0 ;
416 beta_comp = 0 ;
417 d_beta_auto = 0 ;
418 d_beta_comp = 0 ;
419 shift_comp = 0 ;
422 akcar_auto = 0 ;
423 akcar_comp = 0 ;
424 bsn = 0 ;
425 pot_centri = 0 ;
426
427 // Pointers of derived quantities initialized to zero
428 // --------------------------------------------------
429 set_der_0x0() ;
430
431}
432
433 //------------//
434 // Destructor //
435 //------------//
436
438
439 del_deriv() ;
440
441}
442
443 //----------------------------------//
444 // Management of derived quantities //
445 //----------------------------------//
446
448
450
451 if (p_xa_barycenter != 0x0) delete p_xa_barycenter ;
452
453 set_der_0x0() ;
454}
455
456
457
458
460
462
463 p_xa_barycenter = 0x0 ;
464
465}
466
468
470
471 del_deriv() ;
472
473}
474
475
476 //--------------//
477 // Assignment //
478 //--------------//
479
480// Assignment to another Etoile_bin
481// --------------------------------
483
484 // Assignment of quantities common to the derived classes of Etoile
485 Etoile::operator=(et) ;
486
487 // Assignement of proper quantities of class Etoile_bin
489
490 assert(et.ref_triad == ref_triad) ;
491
492 psi0 = et.psi0 ;
493 d_psi = et.d_psi ;
494 wit_w = et.wit_w ;
495 loggam = et.loggam ;
496 logn_comp = et.logn_comp ;
500 beta_comp = et.beta_comp ;
505 w_shift = et.w_shift ;
506 khi_shift = et.khi_shift ;
507 tkij_auto = et.tkij_auto ;
508 tkij_comp = et.tkij_comp ;
511 bsn = et.bsn ;
514 ssjm1_beta = et.ssjm1_beta ;
515 ssjm1_khi = et.ssjm1_khi ;
517 ssjm1_psi = et.ssjm1_psi ;
518 decouple = et.decouple ;
519
520 del_deriv() ; // Deletes all derived quantities
521
522}
523
525
526 del_deriv() ; // sets to 0x0 all the derived quantities
527 return logn_comp ;
528
529}
530
532
533 del_deriv() ; // sets to 0x0 all the derived quantities
534 return pot_centri ;
535
536}
537
539
540 del_deriv() ; // sets to 0x0 all the derived quantities
541 return w_shift ;
542
543}
544
546
547 del_deriv() ; // sets to 0x0 all the derived quantities
548 return khi_shift ;
549
550}
551
552
553 //--------------//
554 // Outputs //
555 //--------------//
556
557// Save in a file
558// --------------
559void Etoile_bin::sauve(FILE* fich) const {
560
561 Etoile::sauve(fich) ;
562
563 fwrite(&irrotational, sizeof(bool), 1, fich) ;
564
565 if (irrotational) {
566 gam_euler.sauve(fich) ; // required to construct d_psi from psi0
567 psi0.sauve(fich) ;
568 }
569
570 w_shift.sauve(fich) ;
571 khi_shift.sauve(fich) ;
572
573 ssjm1_logn.sauve(fich) ;
574 ssjm1_beta.sauve(fich) ;
575 ssjm1_khi.sauve(fich) ;
576 ssjm1_wshift.sauve(fich) ;
577}
578
579// Printing
580// --------
581
582ostream& Etoile_bin::operator>>(ostream& ost) const {
583
584 using namespace Unites ;
585
586 Etoile::operator>>(ost) ;
587
588 ost << endl ;
589 ost << "Star in a binary system" << endl ;
590 ost << "-----------------------" << endl ;
591
592 if (irrotational) {
593 ost << "irrotational configuration" << endl ;
594 }
595 else {
596 ost << "corotating configuration" << endl ;
597 }
598
599 ost << "Absolute abscidia of the stellar center: " <<
600 mp.get_ori_x() / km << " km" << endl ;
601
602 ost << "Absolute abscidia of the barycenter of the baryon density : " <<
603 xa_barycenter() / km << " km" << endl ;
604
605 double r_0 = 0.5 * ( ray_eq() + ray_eq_pi() ) ;
606 double d_ns = fabs( mp.get_ori_x() ) + ray_eq_pi() - r_0 ;
607 double d_tilde = 2 * d_ns / r_0 ;
608
609 ost << "d_tilde : " << d_tilde << endl ;
610
611 ost << "Orientation with respect to the absolute frame : " <<
612 mp.get_rot_phi() << " rad" << endl ;
613
614 ost << "Central value of gam_euler : "
615 << gam_euler()(0, 0, 0, 0) << endl ;
616
617 ost << "Central u_euler (U^X, U^Y, U^Z) [c] : "
618 << u_euler(0)(0, 0, 0, 0) << " "
619 << u_euler(1)(0, 0, 0, 0) << " "
620 << u_euler(2)(0, 0, 0, 0) << endl ;
621
622 if (irrotational) {
623 ost << "Central d_psi (X, Y, Z) [c] : "
624 << d_psi(0)(0, 0, 0, 0) << " "
625 << d_psi(1)(0, 0, 0, 0) << " "
626 << d_psi(2)(0, 0, 0, 0) << endl ;
627
628 ost << "Central vel. / co-orb. (W^X, W^Y, W^Z) [c] : "
629 << wit_w(0)(0, 0, 0, 0) << " "
630 << wit_w(1)(0, 0, 0, 0) << " "
631 << wit_w(2)(0, 0, 0, 0) << endl ;
632
633 ost << "Max vel. / co-orb. (W^X, W^Y, W^Z) [c] : "
634 << max(max(wit_w(0))) << " "
635 << max(max(wit_w(1))) << " "
636 << max(max(wit_w(2))) << endl ;
637
638 ost << "Min vel. / co-orb. (W^X, W^Y, W^Z) [c] : "
639 << min(min(wit_w(0))) << " "
640 << min(min(wit_w(1))) << " "
641 << min(min(wit_w(2))) << endl ;
642
643 double r_surf = mp.val_r(0,1.,M_PI/4,M_PI/4) ;
644
645 ost << "Velocity at (r_surf,pi/4,pi/4) / co-orb. [c] : "
646 << wit_w(0).val_point(r_surf,M_PI/4,M_PI/4) << " "
647 << wit_w(1).val_point(r_surf,M_PI/4,M_PI/4) << " "
648 << wit_w(2).val_point(r_surf,M_PI/4,M_PI/4) << endl ;
649
650 ost << "Central value of loggam : "
651 << loggam()(0, 0, 0, 0) << endl ;
652 }
653
654
655 ost << "Central value of log(N) auto, comp : "
656 << logn_auto()(0, 0, 0, 0) << " "
657 << logn_comp()(0, 0, 0, 0) << endl ;
658
659 ost << "Central value of beta=log(AN) auto, comp : "
660 << beta_auto()(0, 0, 0, 0) << " "
661 << beta_comp()(0, 0, 0, 0) << endl ;
662
663 ost << "Central value of shift (N^X, N^Y, N^Z) [c] : "
664 << shift(0)(0, 0, 0, 0) << " "
665 << shift(1)(0, 0, 0, 0) << " "
666 << shift(2)(0, 0, 0, 0) << endl ;
667
668 ost << " ... shift_auto part of it [c] : "
669 << shift_auto(0)(0, 0, 0, 0) << " "
670 << shift_auto(1)(0, 0, 0, 0) << " "
671 << shift_auto(2)(0, 0, 0, 0) << endl ;
672
673 ost << " ... shift_comp part of it [c] : "
674 << shift_comp(0)(0, 0, 0, 0) << " "
675 << shift_comp(1)(0, 0, 0, 0) << " "
676 << shift_comp(2)(0, 0, 0, 0) << endl ;
677
678 ost << " ... w_shift (NB: components in the star Cartesian frame) [c] : "
679 << endl
680 << w_shift(0)(0, 0, 0, 0) << " "
681 << w_shift(1)(0, 0, 0, 0) << " "
682 << w_shift(2)(0, 0, 0, 0) << endl ;
683
684 ost << "Central value of khi_shift [km c] : "
685 << khi_shift()(0, 0, 0, 0) / km << endl ;
686
687 ost << endl << "Central value of (B^X, B^Y, B^Z)/N [c] : "
688 << bsn(0)(0, 0, 0, 0) << " "
689 << bsn(1)(0, 0, 0, 0) << " "
690 << bsn(2)(0, 0, 0, 0) << endl ;
691
692 ost << endl <<
693 "Central (d/dX,d/dY,d/dZ)(logn_auto) [km^{-1}] : "
694 << d_logn_auto(0)(0, 0, 0, 0) * km << " "
695 << d_logn_auto(1)(0, 0, 0, 0) * km << " "
696 << d_logn_auto(2)(0, 0, 0, 0) * km << endl ;
697
698 ost << "Central (d/dX,d/dY,d/dZ)(logn_comp) [km^{-1}] : "
699 << d_logn_comp(0)(0, 0, 0, 0) * km << " "
700 << d_logn_comp(1)(0, 0, 0, 0) * km << " "
701 << d_logn_comp(2)(0, 0, 0, 0) * km << endl ;
702
703 ost << endl <<
704 "Central (d/dX,d/dY,d/dZ)(beta_auto) [km^{-1}] : "
705 << d_beta_auto(0)(0, 0, 0, 0) * km << " "
706 << d_beta_auto(1)(0, 0, 0, 0) * km << " "
707 << d_beta_auto(2)(0, 0, 0, 0) * km << endl ;
708
709 ost << "Central (d/dX,d/dY,d/dZ)(beta_comp) [km^{-1}] : "
710 << d_beta_comp(0)(0, 0, 0, 0) * km << " "
711 << d_beta_comp(1)(0, 0, 0, 0) * km << " "
712 << d_beta_comp(2)(0, 0, 0, 0) * km << endl ;
713
714
715 ost << endl << "Central A^2 K^{ij} [c/km] : " << endl ;
716 ost << " A^2 K^{xx} auto, comp : "
717 << tkij_auto(0, 0)(0, 0, 0, 0) * km << " "
718 << tkij_comp(0, 0)(0, 0, 0, 0) * km << endl ;
719 ost << " A^2 K^{xy} auto, comp : "
720 << tkij_auto(0, 1)(0, 0, 0, 0) * km << " "
721 << tkij_comp(0, 1)(0, 0, 0, 0) * km << endl ;
722 ost << " A^2 K^{xz} auto, comp : "
723 << tkij_auto(0, 2)(0, 0, 0, 0) * km << " "
724 << tkij_comp(0, 2)(0, 0, 0, 0) * km << endl ;
725 ost << " A^2 K^{yy} auto, comp : "
726 << tkij_auto(1, 1)(0, 0, 0, 0) * km << " "
727 << tkij_comp(1, 1)(0, 0, 0, 0) * km << endl ;
728 ost << " A^2 K^{yz} auto, comp : "
729 << tkij_auto(1, 2)(0, 0, 0, 0) * km << " "
730 << tkij_comp(1, 2)(0, 0, 0, 0) * km << endl ;
731 ost << " A^2 K^{zz} auto, comp : "
732 << tkij_auto(2, 2)(0, 0, 0, 0) * km << " "
733 << tkij_comp(2, 2)(0, 0, 0, 0) * km << endl ;
734
735 ost << endl << "Central A^2 K_{ij} K^{ij} [c^2/km^2] : " << endl ;
736 ost << " A^2 K_{ij} K^{ij} auto, comp : "
737 << akcar_auto()(0, 0, 0, 0) * km*km << " "
738 << akcar_comp()(0, 0, 0, 0) * km*km << endl ;
739
740
741 return ost ;
742}
743
744 //-------------------------//
745 // Computational routines //
746 //-------------------------//
747
748Tenseur Etoile_bin::sprod(const Tenseur& t1, const Tenseur& t2) const {
749
750 int val1 = t1.get_valence() ;
751
752 // Both indices should be contravariant or both covariant :
753 if (t1.get_type_indice(val1-1) == CON) {
754 assert( t2.get_type_indice(0) == CON ) ;
755 return a_car * flat_scalar_prod(t1, t2) ;
756 }
757 else{
758 assert(t1.get_type_indice(val1-1) == COV) ;
759 assert( t2.get_type_indice(0) == COV ) ;
760 return flat_scalar_prod(t1, t2)/a_car ;
761 }
762}
763
765
766 if (!irrotational) {
768 return ;
769 }
770
771 // Specific relativistic enthalpy ---> hhh
772 //----------------------------------
773
774 Tenseur hhh = exp(unsurc2 * ent) ; // = 1 at the Newtonian limit
775
776 // Computation of W^i = - A^2 h Gamma_n B^i/N
777 //----------------------------------------------
778
779 Tenseur www = - a_car * hhh * gam_euler * bsn ;
780
781
782 // Constant value of W^i at the center of the star
783 //-------------------------------------------------
784
785 Tenseur v_orb(mp, 1, COV, ref_triad) ;
786
787 v_orb.set_etat_qcq() ;
788 for (int i=0; i<3; i++) {
789 v_orb.set(i) = www(i)(0, 0, 0, 0) ;
790 }
791
792 v_orb.set_triad( *(www.get_triad()) ) ;
793 int nzm1 = mp.get_mg()->get_nzone() - 1 ;
794 for (int l=nzet; l<=nzm1; l++)
795 for (int i=0; i<=2; i++)
796 v_orb.set(i).annule(l) ;
797
798
799 // Gradient of psi
800 //----------------
801
802 Tenseur d_psi0 = psi0.gradient() ;
803
804 d_psi0.change_triad( ref_triad ) ;
805
806 d_psi = d_psi0 + v_orb ;
807
808
809 // C^1 continuation of d_psi outside the star
810 // (to ensure a smooth enthalpy field accross the stellar surface)
811 // ----------------------------------------------------------------
812
813 if (d_psi0.get_etat() == ETATQCQ ) {
814 d_psi.annule(nzet, nzm1) ;
815 for (int i=0; i<3; i++) {
816 d_psi.set(i).va.set_base( d_psi0(i).va.base ) ;
817 d_psi.set(i) = raccord_c1(d_psi(i), nzet) ;
818 }
819 }
820 else{
821 d_psi.annule(nzm1) ;
822 }
823}
824
825
827
828 Tenseur d_khi = khi_shift.gradient() ;
829
830 if (d_khi.get_etat() == ETATQCQ) {
831 d_khi.dec2_dzpuis() ; // divide by r^2 in the external compactified
832 // domain
833 }
834
835 // x_k dW^k/dx_i
836
837 Tenseur x_d_w = skxk( w_shift.gradient() ) ;
838 x_d_w.dec_dzpuis() ;
839
840 double lambda = double(1) / double(3) ;
841
842 // The final computation is done component by component because
843 // d_khi and x_d_w are covariant comp. whereas w_shift is
844 // contravariant
845
847
848 for (int i=0; i<3; i++) {
849 shift_auto.set(i) = (lambda+2)/2./(lambda+1) * w_shift(i)
850 - (lambda/2./(lambda+1)) * (d_khi(i) + x_d_w(i)) ;
851 }
852
854
855 // The final components of shift_auto should be those with respect
856 // to the absolute frame (X,Y,Z) :
857
859
860}
861
862
863void Etoile_bin::relaxation(const Etoile_bin& star_jm1, double relax_ent,
864 double relax_met, int mer, int fmer_met) {
865
866 double relax_ent_jm1 = 1. - relax_ent ;
867 double relax_met_jm1 = 1. - relax_met ;
868
869 ent = relax_ent * ent + relax_ent_jm1 * star_jm1.ent ;
870
871 if ( (mer != 0) && (mer % fmer_met == 0)) {
872
873 logn_auto = relax_met * logn_auto + relax_met_jm1 * star_jm1.logn_auto ;
874
875 logn_auto_regu = relax_met * logn_auto_regu
876 + relax_met_jm1 * star_jm1.logn_auto_regu ;
877
878 logn_auto_div = relax_met * logn_auto_div
879 + relax_met_jm1 * star_jm1.logn_auto_div ;
880
881 d_logn_auto_div = relax_met * d_logn_auto_div
882 + relax_met_jm1 * star_jm1.d_logn_auto_div ;
883
884 beta_auto = relax_met * beta_auto + relax_met_jm1 * star_jm1.beta_auto ;
885
886 shift_auto = relax_met * shift_auto
887 + relax_met_jm1 * star_jm1.shift_auto ;
888
889 }
890
891 del_deriv() ;
892
894
895}
896}
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition base_vect.h:105
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
void sauve(FILE *) const
Save in a file.
Definition cmp.C:561
Valeur va
The numerical value of the Cmp
Definition cmp.h:464
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:348
Equation of state base class.
Definition eos.h:190
Class for stars in binary system.
Definition etoile.h:814
Tenseur shift_comp
Part of the shift vector generated principaly by the companion star.
Definition etoile.h:895
Tenseur psi0
Scalar potential of the non-translational part of fluid 4-velocity (in the irrotational case)
Definition etoile.h:833
Tenseur d_beta_auto
Gradient of beta_auto (Cartesian components with respect to ref_triad )
Definition etoile.h:879
Tenseur d_logn_auto
Gradient of logn_auto (Cartesian components with respect to ref_triad )
Definition etoile.h:859
double * p_xa_barycenter
Absolute coordinate X of the barycenter of the baryon density.
Definition etoile.h:1006
virtual void sauve(FILE *) const
Save in a file.
Definition etoile_bin.C:559
Tenseur bsn
3-vector shift, divided by N , of the rotating coordinates, .
Definition etoile.h:950
Tenseur akcar_comp
Part of the scalar generated by shift_auto and shift_comp , i.e.
Definition etoile.h:944
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer.
Definition etoile_bin.C:467
Cmp ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Definition etoile.h:973
Tenseur d_psi
Gradient of (in the irrotational case) (Cartesian components with respect to ref_triad )
Definition etoile.h:838
void relaxation(const Etoile_bin &star_prev, double relax_ent, double relax_met, int mer, int fmer_met)
Performs a relaxation on ent , logn_auto , beta_auto and shift_auto .
Definition etoile_bin.C:863
const Base_vect & ref_triad
Reference triad ("absolute frame"), with respect to which the components of all the member Tenseur 's...
Definition etoile.h:828
void fait_shift_auto()
Computes shift_auto from w_shift and khi_shift according to Shibata's prescription [Prog.
Definition etoile_bin.C:826
bool irrotational
true for an irrotational star, false for a corotating one
Definition etoile.h:822
Tenseur w_shift
Vector used in the decomposition of shift_auto , following Shibata's prescription [Prog.
Definition etoile.h:908
Tenseur logn_comp
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principaly by ...
Definition etoile.h:854
Tenseur & set_w_shift()
Read/write of w_shift.
Definition etoile_bin.C:538
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition etoile_bin.C:459
Tenseur d_logn_comp
Gradient of logn_comp (Cartesian components with respect to ref_triad )
Definition etoile.h:869
Tenseur beta_comp
Part of the logarithm of AN generated principaly by the companion star.
Definition etoile.h:874
Tenseur ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for by means...
Definition etoile.h:983
Cmp ssjm1_logn
Effective source at the previous step for the resolution of the Poisson equation for logn_auto by mea...
Definition etoile.h:959
virtual void del_deriv() const
Deletes all the derived quantities.
Definition etoile_bin.C:447
Cmp ssjm1_psi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Definition etoile.h:989
Tenseur pot_centri
Centrifugal potential.
Definition etoile.h:953
Tenseur loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition etoile.h:849
Tenseur shift_auto
Part of the shift vector generated principaly by the star.
Definition etoile.h:889
virtual double xa_barycenter() const
Absolute coordinate X of the barycenter of the baryon density, defined according to the formula.
Tenseur_sym tkij_auto
Part of the extrinsic curvature tensor generated by shift_auto .
Definition etoile.h:925
Tenseur akcar_auto
Part of the scalar generated by shift_auto , i.e.
Definition etoile.h:938
Cmp ssjm1_beta
Effective source at the previous step for the resolution of the Poisson equation for beta_auto by mea...
Definition etoile.h:965
Cmp decouple
Function used to construct the part of generated by the star from the total .
Definition etoile.h:1000
Tenseur wit_w
Spatial projection of the fluid 3-velocity with respect to the co-orbiting observer.
Definition etoile.h:844
Tenseur & set_logn_comp()
Read/write the part of the lapse logarithm (gravitational potential at the Newtonian limit) generated...
Definition etoile_bin.C:524
Tenseur_sym tkij_comp
Part of the extrinsic curvature tensor generated by shift_comp .
Definition etoile.h:932
Tenseur & set_khi_shift()
Read/write of khi_shift.
Definition etoile_bin.C:545
virtual Tenseur sprod(const Tenseur &t1, const Tenseur &t2) const
Performs the scalar product of two tensors by contracting the last index of t1 with the first index o...
Definition etoile_bin.C:748
Tenseur d_beta_comp
Gradient of beta_comp (Cartesian components with respect to ref_triad )
Definition etoile.h:884
Etoile_bin(Map &mp_i, int nzet_i, bool relat, const Eos &eos_i, bool irrot, const Base_vect &ref_triad_i)
Standard constructor.
Definition etoile_bin.C:207
Tenseur d_logn_auto_regu
Gradient of logn_auto_regu (Cartesian components with respect to ref_triad )
Definition etoile.h:864
virtual ~Etoile_bin()
Destructor.
Definition etoile_bin.C:437
Tenseur & set_pot_centri()
Read/write the centrifugal potential.
Definition etoile_bin.C:531
void fait_d_psi()
Computes the gradient of the total velocity potential .
Definition etoile_bin.C:764
void operator=(const Etoile_bin &)
Assignment to another Etoile_bin.
Definition etoile_bin.C:482
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition etoile_bin.C:582
Tenseur khi_shift
Scalar used in the decomposition of shift_auto , following Shibata's prescription [Prog.
Definition etoile.h:918
Base class for stars *** DEPRECATED : use class Star instead ***.
Definition etoile.h:424
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition etoile.C:396
Tenseur logn_auto_div
Divergent part (if k_div!=0 ) of the logarithm of the part of the lapse N generated principaly by t...
Definition etoile.h:497
int nzet
Number of domains of *mp occupied by the star.
Definition etoile.h:432
void operator=(const Etoile &)
Assignment to another Etoile.
Definition etoile.C:430
Tenseur logn_auto_regu
Regular part of the logarithm of the part of the lapse N generated principaly by the star.
Definition etoile.h:491
double ray_eq_pi() const
Coordinate radius at , [r_unit].
double ray_eq() const
Coordinate radius at , [r_unit].
Tenseur logn_auto
Total of the logarithm of the part of the lapse N generated principaly by the star.
Definition etoile.h:484
virtual void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition etoile.C:566
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition etoile.h:474
Tenseur gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition etoile.h:471
Map & mp
Mapping associated with the star.
Definition etoile.h:429
virtual void del_deriv() const
Deletes all the derived quantities.
Definition etoile.C:378
Tenseur d_logn_auto_div
Gradient of logn_auto_div (if k_div!=0 )
Definition etoile.h:501
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition etoile.C:511
virtual void sauve(FILE *) const
Save in a file.
Definition etoile.C:483
Tenseur shift
Total shift vector.
Definition etoile.h:512
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Definition etoile.h:457
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer.
Definition etoile.C:410
Tenseur beta_auto
Logarithm of the part of the product AN generated principaly by by the star.
Definition etoile.h:506
Tenseur a_car
Total conformal factor .
Definition etoile.h:515
double unsurc2
: unsurc2=1 for a relativistic star, 0 for a Newtonian one.
Definition etoile.h:442
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_ori_x() const
Returns the x coordinate of the origin.
Definition map.h:768
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition map.h:775
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
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
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
int get_type_indice(int i) const
Returns the type of the index number i .
Definition tenseur.h:726
void sauve(FILE *) const
Save in a file.
Definition tenseur.C:1325
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tenseur.C:636
void annule(int l)
Sets the Tenseur to zero in a given domain.
Definition tenseur.C:900
void dec2_dzpuis()
dzpuis -= 2 ;
Definition tenseur.C:1130
void dec_dzpuis()
dzpuis -= 1 ;
Definition tenseur.C:1104
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition tenseur.C:1542
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition tenseur.C:645
int get_valence() const
Returns the valence.
Definition tenseur.h:710
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 set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition tenseur.C:674
void set_etat_nondef()
Sets the logical state to ETATNONDEF (undefined state).
Definition tenseur.C:650
int get_etat() const
Returns the logical state.
Definition tenseur.h:707
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition valeur.C:810
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition cmp_math.C:458
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition cmp_math.C:435
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...
Tenseur skxk(const Tenseur &)
Contraction of the last index of (*this) with or , depending on the type of S .
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.