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