LORENE
etoile_rot.C
1 /*
2  * Methods for the class Etoile_rot
3  *
4  * (see file etoile.h for documentation)
5  */
6 
7 /*
8  * Copyright (c) 2000-2001 Eric Gourgoulhon
9  *
10  * This file is part of LORENE.
11  *
12  * LORENE is free software; you can redistribute it and/or modify
13  * it under the terms of the GNU General Public License as published by
14  * the Free Software Foundation; either version 2 of the License, or
15  * (at your option) any later version.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 
29 char etoile_rot_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/etoile_rot.C,v 1.7 2015/12/03 14:17:24 j_novak Exp $" ;
30 
31 /*
32  * $Id: etoile_rot.C,v 1.7 2015/12/03 14:17:24 j_novak Exp $
33  * $Log: etoile_rot.C,v $
34  * Revision 1.7 2015/12/03 14:17:24 j_novak
35  * Check added for the computation of area (thanks S. Koeppel).
36  *
37  * Revision 1.6 2015/06/10 14:37:44 a_sourie
38  * Corrected the formula for the quadrupole.
39  *
40  * Revision 1.5 2014/10/13 08:52:59 j_novak
41  * Lorene classes and functions now belong to the namespace Lorene.
42  *
43  * Revision 1.4 2004/03/25 10:29:07 j_novak
44  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
45  *
46  * Revision 1.3 2001/12/06 15:11:43 jl_zdunik
47  * Introduction of the new function f_eq() in the class Etoile_rot
48  *
49  * Revision 1.2 2001/12/04 21:27:53 e_gourgoulhon
50  *
51  * All writing/reading to a binary file are now performed according to
52  * the big endian convention, whatever the system is big endian or
53  * small endian, thanks to the functions fwrite_be and fread_be
54  *
55  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
56  * LORENE
57  *
58  * Revision 2.17 2001/10/24 15:36:20 eric
59  * Ajout de la fonction display_poly.
60  *
61  * Revision 2.16 2001/10/16 14:49:02 eric
62  * Appel de get_omega_c() pour avoir la valeur centrale de Omega.
63  * Affichage different si rotation differentielle.
64  *
65  * Revision 2.15 2001/09/13 08:32:01 eric
66  * Ajout du facteur de compacite M/R dans l'affichage.
67  *
68  * Revision 2.14 2001/06/20 14:20:56 novak
69  * Appel a Etoile_rot::set_der0x0 dans del_deriv (au lieu de set_der0x0
70  * tout court).
71  *
72  * Revision 2.13 2001/03/26 09:30:58 jlz
73  * New members p_espec_isco and p_lspec_isco.
74  *
75  * Revision 2.12 2000/11/20 21:42:02 eric
76  * Appel de fait_nphi() dans le constructeur par lecture de fichier.
77  *
78  * Revision 2.11 2000/11/18 23:18:30 eric
79  * Modifs affichage.
80  *
81  * Revision 2.10 2000/11/18 21:09:57 eric
82  * Ajout des membres p_r_isco et p_f_isco.
83  *
84  * Revision 2.9 2000/11/07 16:33:08 eric
85  * Modif affichage.
86  *
87  * Revision 2.8 2000/10/12 15:37:01 eric
88  * Ajout de la fonction fait_nphi().
89  *
90  * Revision 2.7 2000/09/18 16:15:12 eric
91  * Ajout du membre tkij.
92  *
93  * Revision 2.6 2000/08/31 15:38:00 eric
94  * Bases spectrales standards pour bbb et b_car dans le constructeur
95  * standard (initialisation a la metrique plate).
96  *
97  * Revision 2.5 2000/08/31 11:25:45 eric
98  * Ajout des membres tnphi et ak_car.
99  *
100  * Revision 2.4 2000/08/25 12:28:29 eric
101  * Modif affichage.
102  *
103  * Revision 2.3 2000/08/18 14:01:59 eric
104  * Ajout de partial_display
105  *
106  * Revision 2.2 2000/08/17 12:40:04 eric
107  * *** empty log message ***
108  *
109  * Revision 2.1 2000/07/21 16:31:26 eric
110  * *** empty log message ***
111  *
112  * Revision 1.1 2000/07/20 15:32:37 eric
113  * Initial revision
114  *
115  *
116  * $Header: /cvsroot/Lorene/C++/Source/Etoile/etoile_rot.C,v 1.7 2015/12/03 14:17:24 j_novak Exp $
117  *
118  */
119 
120 // Headers C
121 #include "math.h"
122 
123 // Headers Lorene
124 #include "etoile.h"
125 #include "eos.h"
126 #include "nbr_spx.h"
127 #include "utilitaires.h"
128 #include "unites.h"
129 
130  //--------------//
131  // Constructors //
132  //--------------//
133 
134 // Standard constructor
135 // --------------------
136 namespace Lorene {
137 Etoile_rot::Etoile_rot(Map& mpi, int nzet_i, bool relat, const Eos& eos_i)
138  : Etoile(mpi, nzet_i, relat, eos_i),
139  bbb(mpi),
140  b_car(mpi),
141  nphi(mpi),
142  tnphi(mpi),
143  uuu(mpi),
144  logn(logn_auto),
145  nuf(mpi),
146  nuq(mpi),
147  dzeta(beta_auto),
148  tggg(mpi),
149  w_shift(mpi, 1, CON, mp.get_bvect_cart()),
150  khi_shift(mpi),
151  tkij(mpi, 2, COV, mp.get_bvect_cart()),
152  ak_car(mpi),
153  ssjm1_nuf(mpi),
154  ssjm1_nuq(mpi),
155  ssjm1_dzeta(mpi),
156  ssjm1_tggg(mpi),
157  ssjm1_khi(mpi),
158  ssjm1_wshift(mpi, 1, CON, mp.get_bvect_cart())
159 {
160 
161  // Initialization to a static state :
162  omega = 0 ;
163  uuu = 0 ;
164 
165  // Initialization to a flat metric :
166  bbb = 1 ;
167  bbb.set_std_base() ;
168  b_car = 1 ;
169  b_car.set_std_base() ;
170  nphi = 0 ;
171  tnphi = 0 ;
172  nuf = 0 ;
173  nuq = 0 ;
174  tggg = 0 ;
175 
176  w_shift.set_etat_qcq() ;
177  for (int i=0; i<3; i++) {
178  w_shift.set(i) = 0 ;
179  }
180 
182  khi_shift.set() = 0 ;
183 
184  tkij.set_etat_zero() ;
185 
186  ak_car = 0 ;
187 
188  ssjm1_nuf = 0 ;
189  ssjm1_nuq = 0 ;
190  ssjm1_dzeta = 0 ;
191  ssjm1_tggg = 0 ;
192  ssjm1_khi = 0 ;
193 
195  for (int i=0; i<3; i++) {
196  ssjm1_wshift.set(i) = 0 ;
197  }
198 
199  // Pointers of derived quantities initialized to zero :
200  set_der_0x0() ;
201 
202 }
203 
204 // Copy constructor
205 // ----------------
206 
208  : Etoile(et),
209  bbb(et.bbb),
210  b_car(et.b_car),
211  nphi(et.nphi),
212  tnphi(et.tnphi),
213  uuu(et.uuu),
214  logn(logn_auto),
215  nuf(et.nuf),
216  nuq(et.nuq),
217  dzeta(beta_auto),
218  tggg(et.tggg),
219  w_shift(et.w_shift),
220  khi_shift(et.khi_shift),
221  tkij(et.tkij),
222  ak_car(et.ak_car),
223  ssjm1_nuf(et.ssjm1_nuf),
224  ssjm1_nuq(et.ssjm1_nuq),
225  ssjm1_dzeta(et.ssjm1_dzeta),
226  ssjm1_tggg(et.ssjm1_tggg),
227  ssjm1_khi(et.ssjm1_khi),
228  ssjm1_wshift(et.ssjm1_wshift)
229 {
230  omega = et.omega ;
231 
232  // Pointers of derived quantities initialized to zero :
233  set_der_0x0() ;
234 }
235 
236 
237 // Constructor from a file
238 // -----------------------
239 Etoile_rot::Etoile_rot(Map& mpi, const Eos& eos_i, FILE* fich)
240  : Etoile(mpi, eos_i, fich),
241  bbb(mpi),
242  b_car(mpi),
243  nphi(mpi),
244  tnphi(mpi),
245  uuu(mpi),
246  logn(logn_auto),
247  nuf(mpi),
248  nuq(mpi),
249  dzeta(beta_auto),
250  tggg(mpi),
251  w_shift(mpi, 1, CON, mp.get_bvect_cart()),
252  khi_shift(mpi),
253  tkij(mpi, 2, COV, mp.get_bvect_cart()),
254  ak_car(mpi),
255  ssjm1_nuf(mpi),
256  ssjm1_nuq(mpi),
257  ssjm1_dzeta(mpi),
258  ssjm1_tggg(mpi),
259  ssjm1_khi(mpi),
260  ssjm1_wshift(mpi, 1, CON, mp.get_bvect_cart())
261 {
262 
263  // Etoile parameters
264  // -----------------
265 
266  // omega is read in the file:
267  fread_be(&omega, sizeof(double), 1, fich) ;
268 
269 
270  // Read of the saved fields:
271  // ------------------------
272 
273  Tenseur nuf_file(mp, fich) ;
274  nuf = nuf_file ;
275 
276  Tenseur nuq_file(mp, fich) ;
277  nuq = nuq_file ;
278 
279  Tenseur tggg_file(mp, fich) ;
280  tggg = tggg_file ;
281 
282  Tenseur w_shift_file(mp, mp.get_bvect_cart(), fich) ;
283  w_shift = w_shift_file ;
284 
285  Tenseur khi_shift_file(mp, fich) ;
286  khi_shift = khi_shift_file ;
287 
288  fait_shift() ; // constructs shift from w_shift and khi_shift
289  fait_nphi() ; // constructs N^phi from (N^x,N^y,N^z)
290 
291  Cmp ssjm1_nuf_file(mp, *(mp.get_mg()), fich) ;
292  ssjm1_nuf = ssjm1_nuf_file ;
293 
294  Cmp ssjm1_nuq_file(mp, *(mp.get_mg()), fich) ;
295  ssjm1_nuq = ssjm1_nuq_file ;
296 
297  Cmp ssjm1_dzeta_file(mp, *(mp.get_mg()), fich) ;
298  ssjm1_dzeta = ssjm1_dzeta_file ;
299 
300  Cmp ssjm1_tggg_file(mp, *(mp.get_mg()), fich) ;
301  ssjm1_tggg = ssjm1_tggg_file ;
302 
303  Cmp ssjm1_khi_file(mp, *(mp.get_mg()), fich) ;
304  ssjm1_khi = ssjm1_khi_file ;
305 
306  Tenseur ssjm1_wshift_file(mp, mp.get_bvect_cart(), fich) ;
307  ssjm1_wshift = ssjm1_wshift_file ;
308 
309  // All other fields are initialized to zero :
310  // ----------------------------------------
311  bbb = 0 ;
312  b_car = 0 ;
313  uuu = 0 ;
314 
315  // Pointers of derived quantities initialized to zero
316  // --------------------------------------------------
317  set_der_0x0() ;
318 
319 }
320 
321  //------------//
322  // Destructor //
323  //------------//
324 
326 
327  del_deriv() ;
328 
329 }
330 
331  //----------------------------------//
332  // Management of derived quantities //
333  //----------------------------------//
334 
335 void Etoile_rot::del_deriv() const {
336 
337  Etoile::del_deriv() ;
338 
339  if (p_angu_mom != 0x0) delete p_angu_mom ;
340  if (p_tsw != 0x0) delete p_tsw ;
341  if (p_grv2 != 0x0) delete p_grv2 ;
342  if (p_grv3 != 0x0) delete p_grv3 ;
343  if (p_r_circ != 0x0) delete p_r_circ ;
344  if (p_area != 0x0) delete p_area ;
345  if (p_aplat != 0x0) delete p_aplat ;
346  if (p_z_eqf != 0x0) delete p_z_eqf ;
347  if (p_z_eqb != 0x0) delete p_z_eqb ;
348  if (p_z_pole != 0x0) delete p_z_pole ;
349  if (p_mom_quad != 0x0) delete p_mom_quad ;
350  if (p_mom_quad_old != 0x0) delete p_mom_quad_old ;
351  if (p_mom_quad_Bo != 0x0) delete p_mom_quad_Bo ;
352  if (p_r_isco != 0x0) delete p_r_isco ;
353  if (p_f_isco != 0x0) delete p_f_isco ;
354  if (p_lspec_isco != 0x0) delete p_lspec_isco ;
355  if (p_espec_isco != 0x0) delete p_espec_isco ;
356  if (p_f_eq != 0x0) delete p_f_eq ;
357 
359 }
360 
361 
362 
363 
365 
367 
368  p_angu_mom = 0x0 ;
369  p_tsw = 0x0 ;
370  p_grv2 = 0x0 ;
371  p_grv3 = 0x0 ;
372  p_r_circ = 0x0 ;
373  p_area = 0x0 ;
374  p_aplat = 0x0 ;
375  p_z_eqf = 0x0 ;
376  p_z_eqb = 0x0 ;
377  p_z_pole = 0x0 ;
378  p_mom_quad = 0x0 ;
379  p_mom_quad_old = 0x0 ;
380  p_mom_quad_Bo = 0x0 ;
381  p_r_isco = 0x0 ;
382  p_f_isco = 0x0 ;
383  p_lspec_isco = 0x0 ;
384  p_espec_isco = 0x0 ;
385  p_f_eq = 0x0 ;
386 
387 }
388 
390 
392 
393  del_deriv() ;
394 
395 }
396 
397 
398  //--------------//
399  // Assignment //
400  //--------------//
401 
402 // Assignment to another Etoile_rot
403 // --------------------------------
405 
406  // Assignment of quantities common to all the derived classes of Etoile
407  Etoile::operator=(et) ;
408 
409  // Assignement of proper quantities of class Etoile_rot
410  omega = et.omega ;
411 
412  bbb = et.bbb ;
413  b_car = et.b_car ;
414  nphi = et.nphi ;
415  tnphi = et.tnphi ;
416  uuu = et.uuu ;
417  nuf = et.nuf ;
418  nuq = et.nuq ;
419  tggg = et.tggg ;
420  w_shift = et.w_shift ;
421  khi_shift = et.khi_shift ;
422  tkij = et.tkij ;
423  ak_car = et.ak_car ;
424  ssjm1_nuf = et.ssjm1_nuf ;
425  ssjm1_nuq = et.ssjm1_nuq ;
426  ssjm1_dzeta = et.ssjm1_dzeta ;
427  ssjm1_tggg = et.ssjm1_tggg ;
428  ssjm1_khi = et.ssjm1_khi ;
429  ssjm1_wshift = et.ssjm1_wshift ;
430 
431  del_deriv() ; // Deletes all derived quantities
432 
433 }
434 
435  //--------------//
436  // Outputs //
437  //--------------//
438 
439 // Save in a file
440 // --------------
441 void Etoile_rot::sauve(FILE* fich) const {
442 
443  Etoile::sauve(fich) ;
444 
445  fwrite_be(&omega, sizeof(double), 1, fich) ;
446 
447  nuf.sauve(fich) ;
448  nuq.sauve(fich) ;
449  tggg.sauve(fich) ;
450  w_shift.sauve(fich) ;
451  khi_shift.sauve(fich) ;
452 
453  ssjm1_nuf.sauve(fich) ;
454  ssjm1_nuq.sauve(fich) ;
455  ssjm1_dzeta.sauve(fich) ;
456  ssjm1_tggg.sauve(fich) ;
457  ssjm1_khi.sauve(fich) ;
458  ssjm1_wshift.sauve(fich) ;
459 
460 
461 }
462 
463 // Printing
464 // --------
465 
466 ostream& Etoile_rot::operator>>(ostream& ost) const {
467 
468  using namespace Unites ;
469 
470  Etoile::operator>>(ost) ;
471 
472  double omega_c = get_omega_c() ;
473 
474  ost << endl ;
475  if (omega != __infinity) {
476  ost << "Uniformly rotating star" << endl ;
477  ost << "-----------------------" << endl ;
478 
479  double freq = omega / (2.*M_PI) ;
480  ost << "Omega : " << omega * f_unit
481  << " rad/s f : " << freq * f_unit << " Hz" << endl ;
482  ost << "Rotation period : " << 1000. / (freq * f_unit) << " ms"
483  << endl ;
484 
485  }
486  else {
487  ost << "Differentially rotating star" << endl ;
488  ost << "----------------------------" << endl ;
489 
490  double freq = omega_c / (2.*M_PI) ;
491  ost << "Central value of Omega : " << omega_c * f_unit
492  << " rad/s f : " << freq * f_unit << " Hz" << endl ;
493  ost << "Central rotation period : " << 1000. / (freq * f_unit) << " ms"
494  << endl ;
495 
496  }
497 
498 
499  double nphi_c = nphi()(0, 0, 0, 0) ;
500  if ( (omega_c==0) && (nphi_c==0) ) {
501  ost << "Central N^phi : " << nphi_c << endl ;
502  }
503  else{
504  ost << "Central N^phi/Omega : " << nphi_c / omega_c << endl ;
505  }
506 
507  ost << "Error on the virial identity GRV2 : " << endl ;
508  ost << "GRV2 = " << grv2() << endl ;
509  ost << "Error on the virial identity GRV3 : " << endl ;
510  double xgrv3 = grv3(&ost) ;
511  ost << "GRV3 = " << xgrv3 << endl ;
512 
513  double mom_quad_38si = mom_quad() * rho_unit * (pow(r_unit, double(5.))
514  / double(1.e38) ) ;
515  ost << "Quadrupole moment Q : " << mom_quad_38si << " 10^38 kg m^2"
516  << endl ;
517  ost << "Q / (M R_circ^2) : "
518  << mom_quad() / ( mass_g() * pow( r_circ(), 2. ) ) << endl ;
519  ost << "c^4 Q / (G^2 M^3) : "
520  << mom_quad() / ( pow(qpig/(4*M_PI), 2.) * pow(mass_g(), 3.) )
521  << endl ;
522 
523  ost << "Angular momentum J : "
524  << angu_mom()/( qpig / (4* M_PI) * msol*msol) << " G M_sol^2 / c"
525  << endl ;
526  ost << "c J / (G M^2) : "
527  << angu_mom()/( qpig / (4* M_PI) * pow(mass_g(), 2.) ) << endl ;
528 
529  if (omega != __infinity) {
530  double mom_iner = angu_mom() / omega ;
531  double mom_iner_38si = mom_iner * rho_unit * (pow(r_unit, double(5.))
532  / double(1.e38) ) ;
533  ost << "Moment of inertia: " << mom_iner_38si << " 10^38 kg m^2"
534  << endl ;
535  }
536 
537  ost << "Ratio T/W : " << tsw() << endl ;
538  ost << "Circumferential equatorial radius R_circ : "
539  << r_circ()/km << " km" << endl ;
540  if (mp.get_mg()->get_np(0) == 1) {
541  ost << "Surface area : " << area()/(km*km) << " km^2" << endl ;
542  ost << "Mean radius R_mean : "
543  << mean_radius()/km << " km" << endl ;
544  } else {
545  ost <<
546  "Skipping surface statements due to number of points in phi direction np == 1"
547  << endl;
548  }
549  ost << "Coordinate equatorial radius r_eq : " << ray_eq()/km << " km"
550  << endl ;
551  ost << "Flattening r_pole/r_eq : " << aplat() << endl ;
552 
553  double compact = qpig/(4.*M_PI) * mass_g() / r_circ() ;
554  ost << "Compaction parameter M_g / R_circ : " << compact << endl ;
555 
556  int lsurf = nzet - 1;
557  int nt = mp.get_mg()->get_nt(lsurf) ;
558  int nr = mp.get_mg()->get_nr(lsurf) ;
559  ost << "Equatorial value of the velocity U: "
560  << uuu()(nzet-1, 0, nt-1, nr-1) << " c" << endl ;
561  ost << "Redshift at the equator (forward) : " << z_eqf() << endl ;
562  ost << "Redshift at the equator (backward): " << z_eqb() << endl ;
563  ost << "Redshift at the pole : " << z_pole() << endl ;
564 
565 
566  ost << "Central value of log(N) : "
567  << logn()(0, 0, 0, 0) << endl ;
568 
569  ost << "Central value of dzeta=log(AN) : "
570  << dzeta()(0, 0, 0, 0) << endl ;
571 
572  if ( (omega_c==0) && (nphi_c==0) ) {
573  ost << "Central N^phi : " << nphi_c << endl ;
574  }
575  else{
576  ost << "Central N^phi/Omega : " << nphi_c / omega_c << endl ;
577  }
578 
579  ost << " ... w_shift (NB: components in the star Cartesian frame) [c] : "
580  << endl
581  << w_shift(0)(0, 0, 0, 0) << " "
582  << w_shift(1)(0, 0, 0, 0) << " "
583  << w_shift(2)(0, 0, 0, 0) << endl ;
584 
585  ost << "Central value of khi_shift [km c] : "
586  << khi_shift()(0, 0, 0, 0) / km << endl ;
587 
588  ost << "Central value of B^2 : " << b_car()(0,0,0,0) << endl ;
589 
590  Tbl diff_a_b = diffrel( a_car(), b_car() ) ;
591  ost <<
592  "Relative discrepancy in each domain between the metric coef. A^2 and B^2 : "
593  << endl ;
594  for (int l=0; l<diff_a_b.get_taille(); l++) {
595  ost << diff_a_b(l) << " " ;
596  }
597  ost << endl ;
598 
599  // Approximate formula for R_isco = 6 R_g (1-(2/3)^1.5 j )
600  // up to the first order in j
601  double jdimless = angu_mom() / ( ggrav * pow(mass_g(), 2.) ) ;
602  double r_grav = ggrav * mass_g() ;
603  double r_isco_appr = 6. * r_grav * ( 1. - pow(2./3.,1.5) * jdimless ) ;
604 
605  // Approximate formula for the ISCO frequency
606  // freq_ms = 6^{-1.5}/2pi/R_g (1+11*6^(-1.5) j )
607  // up to the first order in j
608  double f_isco_appr = ( 1. + 11. /6. /sqrt(6.) * jdimless ) / r_grav /
609  (12. * M_PI ) / sqrt(6.) ;
610 
611  ost << endl << "Innermost stable circular orbit (ISCO) : " << endl ;
612  double xr_isco = r_isco(&ost) ;
613  ost <<" circumferential radius r_isco = "
614  << xr_isco / km << " km" << endl ;
615  ost <<" (approx. 6M + 1st order in j : "
616  << r_isco_appr / km << " km)" << endl ;
617  ost <<" (approx. 6M : "
618  << 6. * r_grav / km << " km)" << endl ;
619  ost <<" orbital frequency f_isco = "
620  << f_isco() * f_unit << " Hz" << endl ;
621  ost <<" (approx. 1st order in j : "
622  << f_isco_appr * f_unit << " Hz)" << endl ;
623 
624 
625  return ost ;
626 
627 }
628 
629 
630 void Etoile_rot::partial_display(ostream& ost) const {
631 
632  using namespace Unites ;
633 
634  double omega_c = get_omega_c() ;
635  double freq = omega_c / (2.*M_PI) ;
636  ost << "Central Omega : " << omega_c * f_unit
637  << " rad/s f : " << freq * f_unit << " Hz" << endl ;
638  ost << "Rotation period : " << 1000. / (freq * f_unit) << " ms"
639  << endl ;
640  ost << endl << "Central enthalpy : " << ent()(0,0,0,0) << " c^2" << endl ;
641  ost << "Central proper baryon density : " << nbar()(0,0,0,0)
642  << " x 0.1 fm^-3" << endl ;
643  ost << "Central proper energy density : " << ener()(0,0,0,0)
644  << " rho_nuc c^2" << endl ;
645  ost << "Central pressure : " << press()(0,0,0,0)
646  << " rho_nuc c^2" << endl ;
647 
648  ost << "Central value of log(N) : "
649  << logn()(0, 0, 0, 0) << endl ;
650  ost << "Central lapse N : " << nnn()(0,0,0,0) << endl ;
651  ost << "Central value of dzeta=log(AN) : "
652  << dzeta()(0, 0, 0, 0) << endl ;
653  ost << "Central value of A^2 : " << a_car()(0,0,0,0) << endl ;
654  ost << "Central value of B^2 : " << b_car()(0,0,0,0) << endl ;
655 
656  double nphi_c = nphi()(0, 0, 0, 0) ;
657  if ( (omega_c==0) && (nphi_c==0) ) {
658  ost << "Central N^phi : " << nphi_c << endl ;
659  }
660  else{
661  ost << "Central N^phi/Omega : " << nphi_c / omega_c
662  << endl ;
663  }
664 
665 
666  int lsurf = nzet - 1;
667  int nt = mp.get_mg()->get_nt(lsurf) ;
668  int nr = mp.get_mg()->get_nr(lsurf) ;
669  ost << "Equatorial value of the velocity U: "
670  << uuu()(nzet-1, 0, nt-1, nr-1) << " c" << endl ;
671 
672  ost << endl
673  << "Coordinate equatorial radius r_eq = "
674  << ray_eq()/km << " km" << endl ;
675  ost << "Flattening r_pole/r_eq : " << aplat() << endl ;
676 
677 }
678 
679 
680 double Etoile_rot::get_omega_c() const {
681 
682  return omega ;
683 
684 }
685 
686 
687 // display_poly
688 // ------------
689 
690 void Etoile_rot::display_poly(ostream& ost) const {
691 
692  using namespace Unites ;
693 
694  const Eos_poly* p_eos_poly = dynamic_cast<const Eos_poly*>( &eos ) ;
695 
696  if (p_eos_poly != 0x0) {
697 
698  double kappa = p_eos_poly->get_kap() ;
699  double gamma = p_eos_poly->get_gam() ; ;
700 
701  // kappa^{n/2}
702  double kap_ns2 = pow( kappa, 0.5 /(gamma-1) ) ;
703 
704  // Polytropic unit of length in terms of r_unit :
705  double r_poly = kap_ns2 / sqrt(ggrav) ;
706 
707  // Polytropic unit of time in terms of t_unit :
708  double t_poly = r_poly ;
709 
710  // Polytropic unit of mass in terms of m_unit :
711  double m_poly = r_poly / ggrav ;
712 
713  // Polytropic unit of angular momentum in terms of j_unit :
714  double j_poly = r_poly * r_poly / ggrav ;
715 
716  // Polytropic unit of density in terms of rho_unit :
717  double rho_poly = 1. / (ggrav * r_poly * r_poly) ;
718 
719  ost.precision(10) ;
720  ost << endl << "Quantities in polytropic units : " << endl ;
721  ost << "==============================" << endl ;
722  ost << " ( r_poly = " << r_poly / km << " km )" << endl ;
723  ost << " n_c : " << nbar()(0, 0, 0, 0) / rho_poly << endl ;
724  ost << " e_c : " << ener()(0, 0, 0, 0) / rho_poly << endl ;
725  ost << " Omega_c : " << get_omega_c() * t_poly << endl ;
726  ost << " P_c : " << 2.*M_PI / get_omega_c() / t_poly << endl ;
727  ost << " M_bar : " << mass_b() / m_poly << endl ;
728  ost << " M : " << mass_g() / m_poly << endl ;
729  ost << " J : " << angu_mom() / j_poly << endl ;
730  ost << " r_eq : " << ray_eq() / r_poly << endl ;
731  ost << " R_circ : " << r_circ() / r_poly << endl ;
732 
733 
734  }
735 
736 
737 }
738 
739 
740 
741 
742 
743 
744  //-------------------------//
745  // Computational routines //
746  //-------------------------//
747 
749 
750  Tenseur d_khi = khi_shift.gradient() ;
751 
752  if (d_khi.get_etat() == ETATQCQ) {
753  d_khi.dec2_dzpuis() ; // divide by r^2 in the external compactified
754  // domain
755  }
756 
757  // x_k dW^k/dx_i
758 
759  Tenseur x_d_w = skxk( w_shift.gradient() ) ;
760  x_d_w.dec_dzpuis() ;
761 
762  double lambda = double(1) / double(3) ;
763 
764  // The final computation is done component by component because
765  // d_khi and x_d_w are covariant comp. whereas w_shift is
766  // contravariant
767 
768  shift.set_etat_qcq() ;
769 
770  for (int i=0; i<3; i++) {
771  shift.set(i) = (lambda+2)/2./(lambda+1) * w_shift(i)
772  - (lambda/2./(lambda+1)) * (d_khi(i) + x_d_w(i)) ;
773  }
774 
775  shift.set_triad( *(w_shift.get_triad()) ) ;
776 
777 }
778 
779 
780 
782 
783  if ( shift.get_etat() == ETATZERO ) {
784  tnphi = 0 ;
785  nphi = 0 ;
786  return ;
787  }
788 
789  assert( shift.get_etat() == ETATQCQ ) ;
790 
791  // Computation of tnphi
792  // --------------------
793  tnphi.set_etat_qcq() ;
794 
795  mp.comp_p_from_cartesian( shift(0), shift(1), tnphi.set() ) ;
796 
797  // Computation of nphi
798  // -------------------
799 
800  nphi = tnphi ;
801  (nphi.set()).div_rsint() ;
802 
803 }
804 }
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
Polytropic equation of state (relativistic case).
Definition: eos.h:757
double get_gam() const
Returns the adiabatic index (cf. Eq. (3))
Definition: eos_poly.C:256
double get_kap() const
Returns the pressure coefficient (cf.
Definition: eos_poly.C:260
Equation of state base class.
Definition: eos.h:190
Class for isolated rotating stars *** DEPRECATED : use class Star_rot instead ***.
Definition: etoile.h:1496
Tenseur ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for .
Definition: etoile.h:1625
Tenseur uuu
Norm of u_euler.
Definition: etoile.h:1518
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: etoile_rot.C:335
double omega
Rotation angular velocity ([f_unit] )
Definition: etoile.h:1501
double * p_z_pole
Redshift factor at North pole.
Definition: etoile.h:1640
virtual double mom_quad() const
Quadrupole moment.
double * p_z_eqf
Forward redshift factor at equator.
Definition: etoile.h:1638
void operator=(const Etoile_rot &)
Assignment to another Etoile_rot.
Definition: etoile_rot.C:404
virtual double r_circ() const
Circumferential radius.
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition: etoile_rot.C:364
Tenseur & logn
Metric potential = logn_auto.
Definition: etoile.h:1521
double * p_mom_quad_old
Part of the quadrupole moment.
Definition: etoile.h:1642
Cmp ssjm1_dzeta
Effective source at the previous step for the resolution of the Poisson equation for dzeta .
Definition: etoile.h:1603
Tenseur nuq
Part of the Metric potential = logn generated by the quadratic terms.
Definition: etoile.h:1531
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: etoile_rot.C:466
double * p_z_eqb
Backward redshift factor at equator.
Definition: etoile.h:1639
virtual double mass_g() const
Gravitational mass.
virtual double f_isco() const
Orbital frequency at the innermost stable circular orbit (ISCO).
Definition: et_rot_isco.C:267
double * p_r_circ
Circumferential radius.
Definition: etoile.h:1635
Tenseur khi_shift
Scalar used in the decomposition of shift , following Shibata's prescription [Prog.
Definition: etoile.h:1560
virtual double aplat() const
Flatening r_pole/r_eq.
virtual double grv3(ostream *ost=0x0) const
Error on the virial identity GRV3.
Etoile_rot(Map &mp_i, int nzet_i, bool relat, const Eos &eos_i)
Standard constructor.
Definition: etoile_rot.C:137
virtual double tsw() const
Ratio T/W.
Tenseur tggg
Metric potential .
Definition: etoile.h:1537
Tenseur nuf
Part of the Metric potential = logn generated by the matter terms.
Definition: etoile.h:1526
Cmp ssjm1_tggg
Effective source at the previous step for the resolution of the Poisson equation for tggg .
Definition: etoile.h:1608
Tenseur nphi
Metric coefficient .
Definition: etoile.h:1510
virtual double mass_b() const
Baryon mass.
virtual double r_isco(ostream *ost=0x0) const
Circumferential radius of the innermost stable circular orbit (ISCO).
Definition: et_rot_isco.C:84
Tenseur bbb
Metric factor B.
Definition: etoile.h:1504
double * p_area
Surface area.
Definition: etoile.h:1636
Tenseur ak_car
Scalar .
Definition: etoile.h:1586
virtual ~Etoile_rot()
Destructor.
Definition: etoile_rot.C:325
Tenseur & dzeta
Metric potential = beta_auto.
Definition: etoile.h:1534
virtual double mean_radius() const
Mean radius.
double * p_grv3
Error on the virial identity GRV3.
Definition: etoile.h:1634
double * p_grv2
Error on the virial identity GRV2.
Definition: etoile.h:1633
double * p_aplat
Flatening r_pole/r_eq.
Definition: etoile.h:1637
Cmp ssjm1_nuf
Effective source at the previous step for the resolution of the Poisson equation for nuf by means of ...
Definition: etoile.h:1592
double * p_mom_quad_Bo
Part of the quadrupole moment.
Definition: etoile.h:1643
virtual double grv2() const
Error on the virial identity GRV2.
void fait_nphi()
Computes tnphi and nphi from the Cartesian components of the shift, stored in shift .
Definition: etoile_rot.C:781
double * p_mom_quad
Quadrupole moment.
Definition: etoile.h:1641
double * p_f_eq
Orbital frequency at the equator.
Definition: etoile.h:1650
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer.
Definition: etoile_rot.C:389
virtual double area() const
Surface area.
double * p_angu_mom
Angular momentum.
Definition: etoile.h:1631
Cmp ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Definition: etoile.h:1616
virtual double angu_mom() const
Angular momentum.
Tenseur_sym tkij
Tensor related to the extrinsic curvature tensor by .
Definition: etoile.h:1567
virtual double z_eqf() const
Forward redshift factor at equator.
virtual void display_poly(ostream &) const
Display in polytropic units.
Definition: etoile_rot.C:690
double * p_f_isco
Orbital frequency of the ISCO.
Definition: etoile.h:1645
double * p_lspec_isco
Specific angular momentum of a particle on the ISCO.
Definition: etoile.h:1649
double * p_tsw
Ratio T/W.
Definition: etoile.h:1632
Tenseur b_car
Square of the metric factor B.
Definition: etoile.h:1507
virtual double z_eqb() const
Backward redshift factor at equator.
Tenseur tnphi
Component of the shift vector.
Definition: etoile.h:1515
double * p_espec_isco
Specific energy of a particle on the ISCO.
Definition: etoile.h:1647
double * p_r_isco
Circumferential radius of the ISCO.
Definition: etoile.h:1644
virtual void sauve(FILE *) const
Save in a file.
Definition: etoile_rot.C:441
virtual double get_omega_c() const
Returns the central value of the rotation angular velocity ([f_unit] )
Definition: etoile_rot.C:680
Cmp ssjm1_nuq
Effective source at the previous step for the resolution of the Poisson equation for nuq by means of ...
Definition: etoile.h:1598
virtual double z_pole() const
Redshift factor at North pole.
Tenseur w_shift
Vector used in the decomposition of shift , following Shibata's prescription [Prog.
Definition: etoile.h:1550
virtual void partial_display(ostream &) const
Printing of some informations, excluding all global quantities.
Definition: etoile_rot.C:630
void fait_shift()
Computes shift from w_shift and khi_shift according to Shibata's prescription [Prog.
Definition: etoile_rot.C:748
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
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
double ray_eq() const
Coordinate radius at , [r_unit].
Tenseur nnn
Total lapse function.
Definition: etoile.h:509
Tenseur nbar
Baryon density in the fluid frame.
Definition: etoile.h:459
const Eos & eos
Equation of state of the stellar matter.
Definition: etoile.h:451
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 ener
Total energy density in the fluid frame.
Definition: etoile.h:460
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: etoile.C:511
Tenseur press
Fluid pressure.
Definition: etoile.h:461
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 a_car
Total conformal factor .
Definition: etoile.h:515
Base class for coordinate mappings.
Definition: map.h:670
virtual void comp_p_from_cartesian(const Scalar &v_x, const Scalar &v_y, Scalar &v_p) const =0
Computes the Spherical component (with respect to bvect_spher ) of a vector given by its cartesian c...
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
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
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Basic array class.
Definition: tbl.h:161
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:301
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
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 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_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1170
void set_etat_zero()
Sets the logical state to ETATZERO (zero state).
Definition: tenseur.C:645
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition: tenseur.C:674
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tenseur.h:704
int get_etat() const
Returns the logical state.
Definition: tenseur.h:707
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition: fread_be.C:69
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition: fwrite_be.C:70
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.