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 
30 char 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
198 namespace Lorene {
199 Cmp raccord_c1(const Cmp& uu, int l1) ;
200 
201  //--------------//
202  // Constructors //
203  //--------------//
204 
205 // Standard constructor
206 // --------------------
207 Etoile_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 
263  w_shift.set_etat_qcq() ;
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),
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),
303  beta_comp(et.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),
320  ssjm1_psi(et.ssjm1_khi),
321  decouple(et.decouple)
322 {
323  set_der_0x0() ;
324 
325 }
326 
327 // Constructor from a file
328 // -----------------------
329 Etoile_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 
447 void Etoile_bin::del_deriv() const {
448 
449  Etoile::del_deriv() ;
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
488  irrotational = et.irrotational ;
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 ;
497  d_logn_auto = et.d_logn_auto ;
499  d_logn_comp = et.d_logn_comp ;
500  beta_comp = et.beta_comp ;
501  d_beta_auto = et.d_beta_auto ;
502  d_beta_comp = et.d_beta_comp ;
503  shift_auto = et.shift_auto ;
504  shift_comp = et.shift_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 ;
509  akcar_auto = et.akcar_auto ;
510  akcar_comp = et.akcar_comp ;
511  bsn = et.bsn ;
512  pot_centri = et.pot_centri ;
513  ssjm1_logn = et.ssjm1_logn ;
514  ssjm1_beta = et.ssjm1_beta ;
515  ssjm1_khi = et.ssjm1_khi ;
516  ssjm1_wshift = et.ssjm1_wshift ;
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 // --------------
559 void 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 
582 ostream& 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 
748 Tenseur 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 
863 void 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 
893  equation_of_state() ;
894 
895 }
896 }
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition: etoile.C:396
Base class for stars *** DEPRECATED : use class Star instead ***.
Definition: etoile.h:424
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: etoile.C:511
int get_type_indice(int i) const
Returns the type of the index number i .
Definition: tenseur.h:723
void annule(int l)
Sets the Tenseur to zero in a given domain.
Definition: tenseur.C:900
Tenseur shift_comp
Part of the shift vector generated principaly by the companion star.
Definition: etoile.h:895
const Base_vect & ref_triad
Reference triad ("absolute frame"), with respect to which the components of all the member Tenseur &#39;s...
Definition: etoile.h:828
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:443
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
virtual void sauve(FILE *) const
Save in a file.
Definition: etoile_bin.C:559
void dec2_dzpuis()
dzpuis -= 2 ;
Definition: tenseur.C:1130
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition: tenseur.C:674
void fait_shift_auto()
Computes shift_auto from w_shift and khi_shift according to Shibata&#39;s prescription [Prog...
Definition: etoile_bin.C:826
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:348
Tenseur pot_centri
Centrifugal potential.
Definition: etoile.h:953
virtual void sauve(FILE *) const
Save in a file.
Definition: etoile.C:483
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
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
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.
Equation of state base class.
Definition: eos.h:190
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
double ray_eq() const
Coordinate radius at , [r_unit].
Base class for coordinate mappings.
Definition: map.h:670
double get_ori_x() const
Returns the x coordinate of the origin.
Definition: map.h:768
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 wit_w
Spatial projection of the fluid 3-velocity with respect to the co-orbiting observer.
Definition: etoile.h:844
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer...
Definition: etoile_bin.C:467
double unsurc2
: unsurc2=1 for a relativistic star, 0 for a Newtonian one.
Definition: etoile.h:442
Class for stars in binary system.
Definition: etoile.h:814
virtual double xa_barycenter() const
Absolute coordinate X of the barycenter of the baryon density, defined according to the formula wher...
double get_rot_phi() const
Returns the angle between the x –axis and X –axis.
Definition: map.h:775
Tenseur d_beta_auto
Gradient of beta_auto (Cartesian components with respect to ref_triad )
Definition: etoile.h:879
Tbl min(const Cmp &)
Minimum values of a Cmp in each domain.
Definition: cmp_math.C:458
Vectorial bases (triads) with respect to which the tensorial components are defined.
Definition: base_vect.h:105
void set_base(const Base_val &)
Sets the bases for spectral expansions (member base )
Definition: valeur.C:810
Tenseur shift
Total shift vector.
Definition: etoile.h:512
Tenseur shift_auto
Part of the shift vector generated principaly by the star.
Definition: etoile.h:889
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
int get_valence() const
Returns the valence.
Definition: tenseur.h:707
Tenseur logn_auto_div
Divergent part (if k_div!=0 ) of the logarithm of the part of the lapse N generated principaly by the...
Definition: etoile.h:497
bool irrotational
true for an irrotational star, false for a corotating one
Definition: etoile.h:822
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: etoile.C:378
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.
Tenseur psi0
Scalar potential of the non-translational part of fluid 4-velocity (in the irrotational case) ...
Definition: etoile.h:833
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
Tenseur w_shift
Vector used in the decomposition of shift_auto , following Shibata&#39;s prescription [Prog...
Definition: etoile.h:908
Tenseur u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: etoile.h:474
Tenseur d_beta_comp
Gradient of beta_comp (Cartesian components with respect to ref_triad )
Definition: etoile.h:884
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
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: etoile_bin.C:447
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition: etoile_bin.C:459
void sauve(FILE *) const
Save in a file.
Definition: tenseur.C:1325
Tenseur logn_comp
Part of the lapse logarithm (gravitational potential at the Newtonian limit) generated principaly by ...
Definition: etoile.h:854
Tenseur skxk(const Tenseur &)
Contraction of the last index of (*this) with or , depending on the type of S .
virtual ~Etoile_bin()
Destructor.
Definition: etoile_bin.C:437
Tenseur d_logn_auto
Gradient of logn_auto (Cartesian components with respect to ref_triad )
Definition: etoile.h:859
Tenseur bsn
3-vector shift, divided by N , of the rotating coordinates, .
Definition: etoile.h:950
void fait_d_psi()
Computes the gradient of the total velocity potential .
Definition: etoile_bin.C:764
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tenseur.h:701
Tenseur gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition: etoile.h:471
Cmp decouple
Function used to construct the part of generated by the star from the total .
Definition: etoile.h:1000
Tenseur akcar_comp
Part of the scalar generated by shift_auto and shift_comp , i.e.
Definition: etoile.h:944
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 beta_comp
Part of the logarithm of AN generated principaly by the companion star.
Definition: etoile.h:874
Map & mp
Mapping associated with the star.
Definition: etoile.h:429
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
virtual void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition: etoile.C:566
Tbl max(const Cmp &)
Maximum values of a Cmp in each domain.
Definition: cmp_math.C:435
int get_etat() const
Returns the logical state.
Definition: tenseur.h:704
Tenseur khi_shift
Scalar used in the decomposition of shift_auto , following Shibata&#39;s prescription [Prog...
Definition: etoile.h:918
Tenseur & set_w_shift()
Read/write of w_shift.
Definition: etoile_bin.C:538
Tenseur_sym tkij_auto
Part of the extrinsic curvature tensor generated by shift_auto .
Definition: etoile.h:925
Tenseur & set_logn_comp()
Read/write the part of the lapse logarithm (gravitational potential at the Newtonian limit) generated...
Definition: etoile_bin.C:524
void operator=(const Etoile_bin &)
Assignment to another Etoile_bin.
Definition: etoile_bin.C:482
int nzet
Number of domains of *mp occupied by the star.
Definition: etoile.h:432
Tenseur a_car
Total conformal factor .
Definition: etoile.h:515
double * p_xa_barycenter
Absolute coordinate X of the barycenter of the baryon density.
Definition: etoile.h:1006
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
const Base_vect_cart & get_bvect_cart() const
Returns the Cartesian basis associated with the coordinates (x,y,z) of the mapping, i.e.
Definition: map.h:791
void dec_dzpuis()
dzpuis -= 1 ;
Definition: tenseur.C:1104
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: etoile_bin.C:582
double ray_eq_pi() 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
Tenseur d_logn_auto_regu
Gradient of logn_auto_regu (Cartesian components with respect to ref_triad )
Definition: etoile.h:864
Tenseur loggam
Logarithm of the Lorentz factor between the fluid and the co-orbiting observer.
Definition: etoile.h:849
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Definition: etoile.h:457
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
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
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tenseur.C:636
void sauve(FILE *) const
Save in a file.
Definition: cmp.C:561
Tenseur & set_pot_centri()
Read/write the centrifugal potential.
Definition: etoile_bin.C:531
Tenseur akcar_auto
Part of the scalar generated by shift_auto , i.e.
Definition: etoile.h:938
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition: tenseur.C:645
Tenseur d_logn_comp
Gradient of logn_comp (Cartesian components with respect to ref_triad )
Definition: etoile.h:869
Tenseur & set_khi_shift()
Read/write of khi_shift.
Definition: etoile_bin.C:545
Tenseur d_logn_auto_div
Gradient of logn_auto_div (if k_div!=0 )
Definition: etoile.h:501
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:461
Tenseur_sym tkij_comp
Part of the extrinsic curvature tensor generated by shift_comp .
Definition: etoile.h:932
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
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:298
void set_etat_nondef()
Sets the logical state to ETATNONDEF (undefined state).
Definition: tenseur.C:650
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition: tenseur.C:1542
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
Tenseur ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for by means...
Definition: etoile.h:983