LORENE
star_rot.C
1 /*
2  * Methods of class Star_rot
3  *
4  * (see file star_rot.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2010 Eric Gourgoulhon
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 version 2
15  * as published by the Free Software Foundation.
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 char star_rot_C[] = "$Header: /cvsroot/Lorene/C++/Source/Star/star_rot.C,v 1.9 2015/05/19 09:30:56 j_novak Exp $" ;
29 
30 /*
31  * $Id: star_rot.C,v 1.9 2015/05/19 09:30:56 j_novak Exp $
32  * $Log: star_rot.C,v $
33  * Revision 1.9 2015/05/19 09:30:56 j_novak
34  * New methods for computing the area of the star and its mean radius.
35  *
36  * Revision 1.8 2014/10/13 08:53:39 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.7 2014/10/06 15:13:17 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.6 2010/02/08 11:45:58 j_novak
43  * Better computation of fait_shift()
44  *
45  * Revision 1.5 2010/02/08 10:56:30 j_novak
46  * Added a few things missing for the reading from resulting file.
47  *
48  * Revision 1.4 2010/02/02 12:45:16 e_gourgoulhon
49  * Improved the display (operator>>)
50  *
51  * Revision 1.3 2010/01/25 22:33:35 e_gourgoulhon
52  * Debugging...
53  *
54  * Revision 1.2 2010/01/25 18:15:32 e_gourgoulhon
55  * Added member unsurc2
56  *
57  * Revision 1.1 2010/01/24 16:09:39 e_gourgoulhon
58  * New class Star_rot.
59  *
60  *
61  * $Header: /cvsroot/Lorene/C++/Source/Star/star_rot.C,v 1.9 2015/05/19 09:30:56 j_novak Exp $
62  *
63  */
64 
65 
66 // C headers
67 #include <cmath>
68 #include <cassert>
69 
70 // Lorene headers
71 #include "star_rot.h"
72 #include "eos.h"
73 #include "unites.h"
74 #include "utilitaires.h"
75 #include "nbr_spx.h"
76 
77 
78 
79  //--------------//
80  // Constructors //
81  //--------------//
82 
83 // Standard constructor
84 // --------------------
85 namespace Lorene {
86 Star_rot::Star_rot(Map& mpi, int nzet_i, bool relat, const Eos& eos_i)
87  : Star(mpi, nzet_i, eos_i),
88  relativistic(relat),
89  a_car(mpi),
90  bbb(mpi),
91  b_car(mpi),
92  nphi(mpi),
93  tnphi(mpi),
94  uuu(mpi),
95  nuf(mpi),
96  nuq(mpi),
97  dzeta(mpi),
98  tggg(mpi),
99  w_shift(mpi, CON, mp.get_bvect_cart()),
100  khi_shift(mpi),
101  tkij(mpi, COV, mp.get_bvect_cart()),
102  ak_car(mpi),
103  ssjm1_nuf(mpi),
104  ssjm1_nuq(mpi),
105  ssjm1_dzeta(mpi),
106  ssjm1_tggg(mpi),
107  ssjm1_khi(mpi),
108  ssjm1_wshift(mpi, CON, mp.get_bvect_cart())
109 {
110 
111  // Parameter 1/c^2 is deduced from relativistic:
112  unsurc2 = relativistic ? double(1) : double(0) ;
113 
114  // Initialization to a static state :
115  omega = 0 ;
116  uuu = 0 ;
117 
118  // Initialization to a flat metric :
119  a_car = 1 ;
120  bbb = 1 ;
122  b_car = 1 ;
123  nphi = 0 ;
124  tnphi = 0 ;
125  nuf = 0 ;
126  nuq = 0 ;
127  dzeta = 0 ;
128  tggg = 0 ;
129 
131  khi_shift = 0 ;
132 
133  beta.set_etat_zero() ;
135 
136  tkij.set_etat_zero() ;
137 
138  ak_car = 0 ;
139 
140  ssjm1_nuf = 0 ;
141  ssjm1_nuq = 0 ;
142  ssjm1_dzeta = 0 ;
143  ssjm1_tggg = 0 ;
144  ssjm1_khi = 0 ;
145 
147 
148  // Pointers of derived quantities initialized to zero :
149  set_der_0x0() ;
150 
151 }
152 
153 // Copy constructor
154 // ----------------
155 
157  : Star(et),
158  relativistic(et.relativistic),
159  unsurc2(et.unsurc2),
160  omega(et.omega),
161  a_car(et.a_car),
162  bbb(et.bbb),
163  b_car(et.b_car),
164  nphi(et.nphi),
165  tnphi(et.tnphi),
166  uuu(et.uuu),
167  nuf(et.nuf),
168  nuq(et.nuq),
169  dzeta(et.dzeta),
170  tggg(et.tggg),
171  w_shift(et.w_shift),
172  khi_shift(et.khi_shift),
173  tkij(et.tkij),
174  ak_car(et.ak_car),
175  ssjm1_nuf(et.ssjm1_nuf),
176  ssjm1_nuq(et.ssjm1_nuq),
177  ssjm1_dzeta(et.ssjm1_dzeta),
178  ssjm1_tggg(et.ssjm1_tggg),
179  ssjm1_khi(et.ssjm1_khi),
180  ssjm1_wshift(et.ssjm1_wshift)
181 {
182  // Pointers of derived quantities initialized to zero :
183  set_der_0x0() ;
184 }
185 
186 
187 // Constructor from a file
188 // -----------------------
189 Star_rot::Star_rot(Map& mpi, const Eos& eos_i, FILE* fich)
190  : Star(mpi, eos_i, fich),
191  a_car(mpi),
192  bbb(mpi),
193  b_car(mpi),
194  nphi(mpi),
195  tnphi(mpi),
196  uuu(mpi),
197  nuf(mpi),
198  nuq(mpi),
199  dzeta(mpi),
200  tggg(mpi),
201  w_shift(mpi, CON, mp.get_bvect_cart()),
202  khi_shift(mpi),
203  tkij(mpi, COV, mp.get_bvect_cart()),
204  ak_car(mpi),
205  ssjm1_nuf(mpi),
206  ssjm1_nuq(mpi),
207  ssjm1_dzeta(mpi),
208  ssjm1_tggg(mpi),
209  ssjm1_khi(mpi),
210  ssjm1_wshift(mpi, CON, mp.get_bvect_cart())
211 {
212 
213  // Star parameters
214  // -----------------
215 
216  // relativistic is read in the file:
217  fread(&relativistic, sizeof(bool), 1, fich) ; //## to be checked !
218 
219  // Parameter 1/c^2 is deduced from relativistic:
220  unsurc2 = relativistic ? double(1) : double(0) ;
221 
222  // omega is read in the file:
223  fread_be(&omega, sizeof(double), 1, fich) ;
224 
225 
226  // Read of the saved fields:
227  // ------------------------
228 
229  Scalar nuf_file(mp, *(mp.get_mg()), fich) ;
230  nuf = nuf_file ;
231 
232  Scalar nuq_file(mp, *(mp.get_mg()), fich) ;
233  nuq = nuq_file ;
234 
235  Scalar dzeta_file(mp, *(mp.get_mg()), fich) ;
236  dzeta = dzeta_file ;
237 
238  Scalar tggg_file(mp, *(mp.get_mg()), fich) ;
239  tggg = tggg_file ;
240 
241  Vector w_shift_file(mp, mp.get_bvect_cart(), fich) ;
242  w_shift = w_shift_file ;
243 
244  Scalar khi_shift_file(mp, *(mp.get_mg()), fich) ;
245  khi_shift = khi_shift_file ;
246 
247  fait_shift() ; // constructs shift from w_shift and khi_shift
248  fait_nphi() ; // constructs N^phi from (N^x,N^y,N^z)
249 
250  Scalar ssjm1_nuf_file(mp, *(mp.get_mg()), fich) ;
251  ssjm1_nuf = ssjm1_nuf_file ;
252 
253  Scalar ssjm1_nuq_file(mp, *(mp.get_mg()), fich) ;
254  ssjm1_nuq = ssjm1_nuq_file ;
255 
256  Scalar ssjm1_dzeta_file(mp, *(mp.get_mg()), fich) ;
257  ssjm1_dzeta = ssjm1_dzeta_file ;
258 
259  Scalar ssjm1_tggg_file(mp, *(mp.get_mg()), fich) ;
260  ssjm1_tggg = ssjm1_tggg_file ;
261 
262  Scalar ssjm1_khi_file(mp, *(mp.get_mg()), fich) ;
263  ssjm1_khi = ssjm1_khi_file ;
264 
265  Vector ssjm1_wshift_file(mp, mp.get_bvect_cart(), fich) ;
266  ssjm1_wshift = ssjm1_wshift_file ;
267 
268  // All other fields are initialized to zero :
269  // ----------------------------------------
270  a_car = 0 ;
271  bbb = 0 ;
272  b_car = 0 ;
273  uuu = 0 ;
274  tkij.set_etat_zero() ;
275  ak_car = 0 ;
276 
277  // Pointers of derived quantities initialized to zero
278  // --------------------------------------------------
279  set_der_0x0() ;
280 
281 }
282 
283  //------------//
284  // Destructor //
285  //------------//
286 
288 
289  del_deriv() ;
290 
291 }
292 
293  //----------------------------------//
294  // Management of derived quantities //
295  //----------------------------------//
296 
297 void Star_rot::del_deriv() const {
298 
299  Star::del_deriv() ;
300 
301  if (p_angu_mom != 0x0) delete p_angu_mom ;
302  if (p_tsw != 0x0) delete p_tsw ;
303  if (p_grv2 != 0x0) delete p_grv2 ;
304  if (p_grv3 != 0x0) delete p_grv3 ;
305  if (p_r_circ != 0x0) delete p_r_circ ;
306  if (p_area != 0x0) delete p_area ;
307  if (p_aplat != 0x0) delete p_aplat ;
308  if (p_z_eqf != 0x0) delete p_z_eqf ;
309  if (p_z_eqb != 0x0) delete p_z_eqb ;
310  if (p_z_pole != 0x0) delete p_z_pole ;
311  if (p_mom_quad != 0x0) delete p_mom_quad ;
312  if (p_r_isco != 0x0) delete p_r_isco ;
313  if (p_f_isco != 0x0) delete p_f_isco ;
314  if (p_lspec_isco != 0x0) delete p_lspec_isco ;
315  if (p_espec_isco != 0x0) delete p_espec_isco ;
316  if (p_f_eq != 0x0) delete p_f_eq ;
317 
319 }
320 
321 
322 void Star_rot::set_der_0x0() const {
323 
325 
326  p_angu_mom = 0x0 ;
327  p_tsw = 0x0 ;
328  p_grv2 = 0x0 ;
329  p_grv3 = 0x0 ;
330  p_r_circ = 0x0 ;
331  p_area = 0x0 ;
332  p_aplat = 0x0 ;
333  p_z_eqf = 0x0 ;
334  p_z_eqb = 0x0 ;
335  p_z_pole = 0x0 ;
336  p_mom_quad = 0x0 ;
337  p_r_isco = 0x0 ;
338  p_f_isco = 0x0 ;
339  p_lspec_isco = 0x0 ;
340  p_espec_isco = 0x0 ;
341  p_f_eq = 0x0 ;
342 
343 }
344 
346 
348 
349  del_deriv() ;
350 
351 }
352 
353  //--------------//
354  // Assignment //
355  //--------------//
356 
357 // Assignment to another Star_rot
358 // --------------------------------
359 void Star_rot::operator=(const Star_rot& et) {
360 
361  // Assignment of quantities common to all the derived classes of Star
362  Star::operator=(et) ;
363 
364  // Assignement of proper quantities of class Star_rot
365  relativistic = et.relativistic ;
366  unsurc2 = et.unsurc2 ;
367  omega = et.omega ;
368 
369  a_car = et.a_car ;
370  bbb = et.bbb ;
371  b_car = et.b_car ;
372  nphi = et.nphi ;
373  tnphi = et.tnphi ;
374  uuu = et.uuu ;
375  nuf = et.nuf ;
376  nuq = et.nuq ;
377  dzeta = et.dzeta ;
378  tggg = et.tggg ;
379  w_shift = et.w_shift ;
380  khi_shift = et.khi_shift ;
381  tkij = et.tkij ;
382  ak_car = et.ak_car ;
383  ssjm1_nuf = et.ssjm1_nuf ;
384  ssjm1_nuq = et.ssjm1_nuq ;
385  ssjm1_dzeta = et.ssjm1_dzeta ;
386  ssjm1_tggg = et.ssjm1_tggg ;
387  ssjm1_khi = et.ssjm1_khi ;
388  ssjm1_wshift = et.ssjm1_wshift ;
389 
390  del_deriv() ; // Deletes all derived quantities
391 
392 }
393 
394 
395  //--------------//
396  // Outputs //
397  //--------------//
398 
399 // Save in a file
400 // --------------
401 void Star_rot::sauve(FILE* fich) const {
402 
403  Star::sauve(fich) ;
404 
405  fwrite(&relativistic, sizeof(bool), 1, fich) ;
406  fwrite_be(&omega, sizeof(double), 1, fich) ;
407 
408  nuf.sauve(fich) ;
409  nuq.sauve(fich) ;
410  dzeta.sauve(fich) ;
411  tggg.sauve(fich) ;
412  w_shift.sauve(fich) ;
413  khi_shift.sauve(fich) ;
414 
415  ssjm1_nuf.sauve(fich) ;
416  ssjm1_nuq.sauve(fich) ;
417  ssjm1_dzeta.sauve(fich) ;
418  ssjm1_tggg.sauve(fich) ;
419  ssjm1_khi.sauve(fich) ;
420  ssjm1_wshift.sauve(fich) ;
421 
422 
423 }
424 
425 // Printing
426 // --------
427 
428 ostream& Star_rot::operator>>(ostream& ost) const {
429 
430  using namespace Unites ;
431 
432  Star::operator>>(ost) ;
433 
434  double omega_c = get_omega_c() ;
435 
436  ost << endl ;
437 
438  if (omega != __infinity) {
439  ost << "Uniformly rotating star" << endl ;
440  ost << "-----------------------" << endl ;
441 
442  double freq = omega / (2.*M_PI) ;
443  ost << "Omega : " << omega * f_unit
444  << " rad/s f : " << freq * f_unit << " Hz" << endl ;
445  ost << "Rotation period : " << 1000. / (freq * f_unit) << " ms"
446  << endl ;
447 
448  }
449  else {
450  ost << "Differentially rotating star" << endl ;
451  ost << "----------------------------" << endl ;
452 
453  double freq = omega_c / (2.*M_PI) ;
454  ost << "Central value of Omega : " << omega_c * f_unit
455  << " rad/s f : " << freq * f_unit << " Hz" << endl ;
456  ost << "Central rotation period : " << 1000. / (freq * f_unit) << " ms"
457  << endl ;
458  }
459  if (relativistic) {
460  ost << "Relativistic star" << endl ;
461  }
462  else {
463  ost << "Newtonian star" << endl ;
464  }
465  double compact = qpig/(4.*M_PI) * mass_g() / r_circ() ;
466  ost << "Compactness G M_g /(c^2 R_circ) : " << compact << endl ;
467 
468  double nphi_c = nphi.val_grid_point(0, 0, 0, 0) ;
469  if ( (omega_c==0) && (nphi_c==0) ) {
470  ost << "Central N^phi : " << nphi_c << endl ;
471  }
472  else{
473  ost << "Central N^phi/Omega : " << nphi_c / omega_c << endl ;
474  }
475 
476  ost << "Error on the virial identity GRV2 : " << grv2() << endl ;
477  double xgrv3 = grv3(&ost) ;
478  ost << "Error on the virial identity GRV3 : " << xgrv3 << endl ;
479 
480  double mom_quad_38si = mom_quad() * rho_unit * (pow(r_unit, double(5.))
481  / double(1.e38) ) ;
482  ost << "Quadrupole moment Q : " << mom_quad_38si << " 10^38 kg m^2"
483  << endl ;
484  ost << "Q / (M R_circ^2) : "
485  << mom_quad() / ( mass_g() * pow( r_circ(), 2. ) ) << endl ;
486  ost << "c^4 Q / (G^2 M^3) : "
487  << mom_quad() / ( pow(qpig/(4*M_PI), 2.) * pow(mass_g(), 3.) )
488  << endl ;
489 
490  ost << "Angular momentum J : "
491  << angu_mom()/( qpig / (4* M_PI) * msol*msol) << " G M_sol^2 / c"
492  << endl ;
493  ost << "c J / (G M^2) : "
494  << angu_mom()/( qpig / (4* M_PI) * pow(mass_g(), 2.) ) << endl ;
495 
496  if (omega != __infinity) {
497  double mom_iner = angu_mom() / omega ;
498  double mom_iner_38si = mom_iner * rho_unit * (pow(r_unit, double(5.))
499  / double(1.e38) ) ;
500  ost << "Moment of inertia: " << mom_iner_38si << " 10^38 kg m^2"
501  << endl ;
502  }
503 
504  ost << "Ratio T/W : " << tsw() << endl ;
505  ost << "Circumferential equatorial radius R_circ : "
506  << r_circ()/km << " km" << endl ;
507  if (mp.get_mg()->get_np(0) == 1) {
508  ost << "Surface area : " << area()/(km*km) << " km^2" << endl ;
509  ost << "Mean radius : " << mean_radius()/km << " km" << endl ;
510  }
511  ost << "Coordinate equatorial radius r_eq : " << ray_eq()/km << " km"
512  << endl ;
513  ost << "Flattening r_pole/r_eq : " << aplat() << endl ;
514 
515 
516  int lsurf = nzet - 1;
517  int nt = mp.get_mg()->get_nt(lsurf) ;
518  int nr = mp.get_mg()->get_nr(lsurf) ;
519  ost << "Equatorial value of the velocity U: "
520  << uuu.val_grid_point(nzet-1, 0, nt-1, nr-1) << " c" << endl ;
521  ost << "Redshift at the equator (forward) : " << z_eqf() << endl ;
522  ost << "Redshift at the equator (backward): " << z_eqb() << endl ;
523  ost << "Redshift at the pole : " << z_pole() << endl ;
524 
525 
526  ost << "Central value of log(N) : "
527  << logn.val_grid_point(0, 0, 0, 0) << endl ;
528 
529  ost << "Central value of dzeta=log(AN) : "
530  << dzeta.val_grid_point(0, 0, 0, 0) << endl ;
531 
532  if ( (omega_c==0) && (nphi_c==0) ) {
533  ost << "Central N^phi : " << nphi_c << endl ;
534  }
535  else{
536  ost << "Central N^phi/Omega : " << nphi_c / omega_c << endl ;
537  }
538 
539  ost << " ... w_shift (NB: components in the star Cartesian frame) [c] : "
540  << endl
541  << w_shift(1).val_grid_point(0, 0, 0, 0) << " "
542  << w_shift(2).val_grid_point(0, 0, 0, 0) << " "
543  << w_shift(3).val_grid_point(0, 0, 0, 0) << endl ;
544 
545  ost << "Central value of khi_shift [km c] : "
546  << khi_shift.val_grid_point(0, 0, 0, 0) / km << endl ;
547 
548  ost << "Central value of B^2 : " << b_car.val_grid_point(0,0,0,0) << endl ;
549 
550  Tbl diff_a_b = diffrel( a_car, b_car ) ;
551  ost <<
552  "Relative discrepancy in each domain between the metric coef. A^2 and B^2 : "
553  << endl ;
554  for (int l=0; l<diff_a_b.get_taille(); l++) {
555  ost << diff_a_b(l) << " " ;
556  }
557  ost << endl ;
558 
559  // Approximate formula for R_isco = 6 R_g (1-(2/3)^1.5 j )
560  // up to the first order in j
561  double jdimless = angu_mom() / ( ggrav * pow(mass_g(), 2.) ) ;
562  double r_grav = ggrav * mass_g() ;
563  double r_isco_appr = 6. * r_grav * ( 1. - pow(2./3.,1.5) * jdimless ) ;
564 
565  // Approximate formula for the ISCO frequency
566  // freq_ms = 6^{-1.5}/2pi/R_g (1+11*6^(-1.5) j )
567  // up to the first order in j
568  double f_isco_appr = ( 1. + 11. /6. /sqrt(6.) * jdimless ) / r_grav /
569  (12. * M_PI ) / sqrt(6.) ;
570 
571  ost << endl << "Innermost stable circular orbit (ISCO) : " << endl ;
572  double xr_isco = r_isco(&ost) ;
573  ost <<" circumferential radius r_isco = "
574  << xr_isco / km << " km" << endl ;
575  ost <<" (approx. 6M + 1st order in j : "
576  << r_isco_appr / km << " km)" << endl ;
577  ost <<" (approx. 6M : "
578  << 6. * r_grav / km << " km)" << endl ;
579  ost <<" orbital frequency f_isco = "
580  << f_isco() * f_unit << " Hz" << endl ;
581  ost <<" (approx. 1st order in j : "
582  << f_isco_appr * f_unit << " Hz)" << endl ;
583 
584 
585  return ost ;
586 
587 }
588 
589 
590 void Star_rot::partial_display(ostream& ost) const {
591 
592  using namespace Unites ;
593 
594  double omega_c = get_omega_c() ;
595  double freq = omega_c / (2.*M_PI) ;
596  ost << "Central Omega : " << omega_c * f_unit
597  << " rad/s f : " << freq * f_unit << " Hz" << endl ;
598  ost << "Rotation period : " << 1000. / (freq * f_unit) << " ms"
599  << endl ;
600  ost << endl << "Central enthalpy : " << ent.val_grid_point(0,0,0,0) << " c^2" << endl ;
601  ost << "Central proper baryon density : " << nbar.val_grid_point(0,0,0,0)
602  << " x 0.1 fm^-3" << endl ;
603  ost << "Central proper energy density : " << ener.val_grid_point(0,0,0,0)
604  << " rho_nuc c^2" << endl ;
605  ost << "Central pressure : " << press.val_grid_point(0,0,0,0)
606  << " rho_nuc c^2" << endl ;
607 
608  ost << "Central value of log(N) : "
609  << logn.val_grid_point(0, 0, 0, 0) << endl ;
610  ost << "Central lapse N : " << nn.val_grid_point(0,0,0,0) << endl ;
611  ost << "Central value of dzeta=log(AN) : "
612  << dzeta.val_grid_point(0, 0, 0, 0) << endl ;
613  ost << "Central value of A^2 : " << a_car.val_grid_point(0,0,0,0) << endl ;
614  ost << "Central value of B^2 : " << b_car.val_grid_point(0,0,0,0) << endl ;
615 
616  double nphi_c = nphi.val_grid_point(0, 0, 0, 0) ;
617  if ( (omega_c==0) && (nphi_c==0) ) {
618  ost << "Central N^phi : " << nphi_c << endl ;
619  }
620  else{
621  ost << "Central N^phi/Omega : " << nphi_c / omega_c
622  << endl ;
623  }
624 
625 
626  int lsurf = nzet - 1;
627  int nt = mp.get_mg()->get_nt(lsurf) ;
628  int nr = mp.get_mg()->get_nr(lsurf) ;
629  ost << "Equatorial value of the velocity U: "
630  << uuu.val_grid_point(nzet-1, 0, nt-1, nr-1) << " c" << endl ;
631 
632  ost << endl
633  << "Coordinate equatorial radius r_eq = "
634  << ray_eq()/km << " km" << endl ;
635  ost << "Flattening r_pole/r_eq : " << aplat() << endl ;
636 
637 }
638 
639 
640 double Star_rot::get_omega_c() const {
641 
642  return omega ;
643 
644 }
645 
646 // display_poly
647 // ------------
648 
649 void Star_rot::display_poly(ostream& ost) const {
650 
651  using namespace Unites ;
652 
653  const Eos_poly* p_eos_poly = dynamic_cast<const Eos_poly*>( &eos ) ;
654 
655  if (p_eos_poly != 0x0) {
656 
657  double kappa = p_eos_poly->get_kap() ;
658 
659  // kappa^{n/2}
660  double kap_ns2 = pow( kappa, 0.5 /(p_eos_poly->get_gam()-1) ) ;
661 
662  // Polytropic unit of length in terms of r_unit :
663  double r_poly = kap_ns2 / sqrt(ggrav) ;
664 
665  // Polytropic unit of time in terms of t_unit :
666  double t_poly = r_poly ;
667 
668  // Polytropic unit of mass in terms of m_unit :
669  double m_poly = r_poly / ggrav ;
670 
671  // Polytropic unit of angular momentum in terms of j_unit :
672  double j_poly = r_poly * r_poly / ggrav ;
673 
674  // Polytropic unit of density in terms of rho_unit :
675  double rho_poly = 1. / (ggrav * r_poly * r_poly) ;
676 
677  ost.precision(10) ;
678  ost << endl << "Quantities in polytropic units : " << endl ;
679  ost << "==============================" << endl ;
680  ost << " ( r_poly = " << r_poly / km << " km )" << endl ;
681  ost << " n_c : " << nbar.val_grid_point(0, 0, 0, 0) / rho_poly << endl ;
682  ost << " e_c : " << ener.val_grid_point(0, 0, 0, 0) / rho_poly << endl ;
683  ost << " Omega_c : " << get_omega_c() * t_poly << endl ;
684  ost << " P_c : " << 2.*M_PI / get_omega_c() / t_poly << endl ;
685  ost << " M_bar : " << mass_b() / m_poly << endl ;
686  ost << " M : " << mass_g() / m_poly << endl ;
687  ost << " J : " << angu_mom() / j_poly << endl ;
688  ost << " r_eq : " << ray_eq() / r_poly << endl ;
689  ost << " R_circ : " << r_circ() / r_poly << endl ;
690  ost << " R_mean : " << mean_radius() / r_poly << endl ;
691 
692  }
693 
694 
695 }
696 
697 
698 
699  //-------------------------//
700  // Computational methods //
701  //-------------------------//
702 
704 
705  Vector d_khi = khi_shift.derive_con( mp.flat_met_cart() ) ;
706 
707  d_khi.dec_dzpuis(2) ; // divide by r^2 in the external compactified domain
708 
709  // x_k dW^k/dx_i
710  Scalar xx(mp) ;
711  Scalar yy(mp) ;
712  Scalar zz(mp) ;
713  Scalar sintcosp(mp) ;
714  Scalar sintsinp(mp) ;
715  Scalar cost(mp) ;
716  xx = mp.x ;
717  yy = mp.y ;
718  zz = mp.z ;
719  sintcosp = mp.sint * mp.cosp ;
720  sintsinp = mp.sint * mp.sinp ;
721  cost = mp.cost ;
722 
723  int nz = mp.get_mg()->get_nzone() ;
724  Vector xk(mp, COV, mp.get_bvect_cart()) ;
725  xk.set(1) = xx ;
726  xk.set(1).set_domain(nz-1) = sintcosp.domain(nz-1) ;
727  xk.set(1).set_dzpuis(-1) ;
728  xk.set(2) = yy ;
729  xk.set(2).set_domain(nz-1) = sintsinp.domain(nz-1) ;
730  xk.set(2).set_dzpuis(-1) ;
731  xk.set(3) = zz ;
732  xk.set(3).set_domain(nz-1) = cost.domain(nz-1) ;
733  xk.set(3).set_dzpuis(-1) ;
734  xk.std_spectral_base() ;
735 
737 
738  Vector x_d_w = contract(xk, 0, d_w, 0) ;
739  x_d_w.dec_dzpuis() ;
740 
741  double lambda = double(1) / double(3) ;
742 
743  beta = - (lambda+2)/2./(lambda+1) * w_shift
744  + (lambda/2./(lambda+1)) * (d_khi + x_d_w) ;
745 
746 }
747 
749 
750  if ( (beta(1).get_etat() == ETATZERO) && (beta(2).get_etat() == ETATZERO) ) {
751  tnphi = 0 ;
752  nphi = 0 ;
753  return ;
754  }
755 
756  // Computation of tnphi
757  // --------------------
758 
759  mp.comp_p_from_cartesian( -beta(1), -beta(2), tnphi ) ;
760 
761  // Computation of nphi
762  // -------------------
763 
764  nphi = tnphi ;
765  nphi.div_rsint() ;
766 
767 }
768 
769 }
Lorene::Star_rot::khi_shift
Scalar khi_shift
Scalar used in the decomposition of shift , following Shibata's prescription [Prog.
Definition: star_rot.h:160
Lorene::Map::z
Coord z
z coordinate centered on the grid
Definition: map.h:728
Lorene::Star::nbar
Scalar nbar
Baryon density in the fluid frame.
Definition: star.h:192
Lorene::Star::logn
Scalar logn
Logarithm of the lapse N .
Definition: star.h:222
Lorene::Star_rot::p_lspec_isco
double * p_lspec_isco
Specific angular momentum of a particle on the ISCO.
Definition: star_rot.h:248
Lorene::Star_rot::partial_display
virtual void partial_display(ostream &) const
Printing of some informations, excluding all global quantities.
Definition: star_rot.C:590
Lorene::Star_rot::tnphi
Scalar tnphi
Component of the shift vector.
Definition: star_rot.h:118
Lorene::Star_rot::tsw
virtual double tsw() const
Ratio T/W.
Definition: star_rot_global.C:176
Lorene::Star_rot::operator=
void operator=(const Star_rot &)
Assignment to another Star_rot.
Definition: star_rot.C:359
Lorene::Scalar::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition: scalar.C:686
Lorene::Star_rot::z_eqf
virtual double z_eqf() const
Forward redshift factor at equator.
Definition: star_rot_global.C:443
Lorene::Star_rot
Class for isolated rotating stars.
Definition: star_rot.h:86
Lorene::Mg3d::get_np
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Lorene::Star_rot::ssjm1_khi
Scalar ssjm1_khi
Effective source at the previous step for the resolution of the Poisson equation for the scalar by m...
Definition: star_rot.h:216
Lorene::Star_rot::r_isco
virtual double r_isco(ostream *ost=0x0) const
Circumferential radius of the innermost stable circular orbit (ISCO).
Definition: star_rot_isco.C:71
Lorene::Star_rot::unsurc2
double unsurc2
: unsurc2=1 for a relativistic star, 0 for a Newtonian one.
Definition: star_rot.h:99
Lorene::Star_rot::p_tsw
double * p_tsw
Ratio T/W.
Definition: star_rot.h:233
Lorene::Star::operator>>
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: star.C:417
Lorene::Star_rot::mass_g
virtual double mass_g() const
Gravitational mass.
Definition: star_rot_global.C:124
Lorene::Map::cost
Coord cost
Definition: map.h:722
Lorene::Star::ray_eq
double ray_eq() const
Coordinate radius at , [r_unit].
Definition: star_global.C:108
Lorene::Star::del_deriv
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: star.C:298
Lorene::Star::ener
Scalar ener
Total energy density in the fluid frame.
Definition: star.h:193
Lorene::Star_rot::get_omega_c
virtual double get_omega_c() const
Returns the central value of the rotation angular velocity ([f_unit] )
Definition: star_rot.C:640
Lorene::Star_rot::del_deriv
virtual void del_deriv() const
Deletes all the derived quantities.
Definition: star_rot.C:297
Lorene::Tensor
Tensor handling.
Definition: tensor.h:288
Lorene::Star_rot::p_espec_isco
double * p_espec_isco
Specific energy of a particle on the ISCO.
Definition: star_rot.h:246
Lorene::Star_rot::omega
double omega
Rotation angular velocity ([f_unit] )
Definition: star_rot.h:101
Lorene::Star_rot::p_f_eq
double * p_f_eq
Orbital frequency at the equator.
Definition: star_rot.h:249
Lorene::diffrel
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:504
Lorene::Mg3d::get_nt
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::Star::press
Scalar press
Fluid pressure.
Definition: star.h:194
Lorene::Star::nzet
int nzet
Number of domains of *mp occupied by the star.
Definition: star.h:183
Lorene::Star_rot::p_r_isco
double * p_r_isco
Circumferential radius of the ISCO.
Definition: star_rot.h:243
Lorene::Tensor::set_etat_zero
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition: tensor.C:497
Lorene::Star_rot::nuf
Scalar nuf
Part of the Metric potential = logn generated by the matter terms.
Definition: star_rot.h:126
Lorene::Star_rot::bbb
Scalar bbb
Metric factor B.
Definition: star_rot.h:107
Lorene::Map::get_mg
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
Lorene::Star_rot::Star_rot
Star_rot(Map &mp_i, int nzet_i, bool relat, const Eos &eos_i)
Standard constructor.
Definition: star_rot.C:86
Lorene::Star_rot::relativistic
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one.
Definition: star_rot.h:94
Lorene::Tensor::dec_dzpuis
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values in the compactified ext...
Definition: tensor.C:808
Lorene::Star_rot::fait_shift
void fait_shift()
Computes shift from w_shift and khi_shift according to Shibata's prescription [Prog.
Definition: star_rot.C:703
Lorene::Star_rot::z_pole
virtual double z_pole() const
Redshift factor at North pole.
Definition: star_rot_global.C:502
Lorene::Star_rot::~Star_rot
virtual ~Star_rot()
Destructor.
Definition: star_rot.C:287
Lorene::Map::sint
Coord sint
Definition: map.h:721
Lorene::Star_rot::nphi
Scalar nphi
Metric coefficient .
Definition: star_rot.h:113
Lorene::Star_rot::grv2
virtual double grv2() const
Error on the virial identity GRV2.
Definition: star_rot_global.C:204
Lorene::Eos
Equation of state base class.
Definition: eos.h:190
Lorene::Scalar::div_rsint
void div_rsint()
Division by everywhere; dzpuis is not changed.
Definition: scalar_r_manip.C:348
Lorene::Scalar
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:387
Lorene::Star::nn
Scalar nn
Lapse function N .
Definition: star.h:225
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
Lorene::Star_rot::nuq
Scalar nuq
Part of the Metric potential = logn generated by the quadratic terms.
Definition: star_rot.h:131
Lorene::Star_rot::mass_b
virtual double mass_b() const
Baryon mass.
Definition: star_rot_global.C:98
Lorene::Star_rot::ssjm1_dzeta
Scalar ssjm1_dzeta
Effective source at the previous step for the resolution of the Poisson equation for dzeta .
Definition: star_rot.h:203
Lorene::Star_rot::mom_quad
virtual double mom_quad() const
Quadrupole moment.
Definition: star_rot_global.C:521
Lorene::Star_rot::angu_mom
virtual double angu_mom() const
Angular momentum.
Definition: star_rot_global.C:150
Lorene::Star_rot::ssjm1_tggg
Scalar ssjm1_tggg
Effective source at the previous step for the resolution of the Poisson equation for tggg .
Definition: star_rot.h:208
Unites
Standard units of space, time and mass.
Lorene::Map::sinp
Coord sinp
Definition: map.h:723
Lorene::Star_rot::p_mom_quad
double * p_mom_quad
Quadrupole moment.
Definition: star_rot.h:242
Lorene::Map::comp_p_from_cartesian
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...
Lorene::Star_rot::p_f_isco
double * p_f_isco
Orbital frequency of the ISCO.
Definition: star_rot.h:244
Lorene::Scalar::domain
const Tbl & domain(int l) const
Read-only of the value in a given domain.
Definition: scalar.h:625
Lorene::Star::del_hydro_euler
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer.
Definition: star.C:330
Lorene::Star_rot::tkij
Sym_tensor tkij
Tensor related to the extrinsic curvature tensor by .
Definition: star_rot.h:167
Lorene::Star::mp
Map & mp
Mapping associated with the star.
Definition: star.h:180
Lorene::Vector::set
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:296
Lorene::Star_rot::p_z_eqf
double * p_z_eqf
Forward redshift factor at equator.
Definition: star_rot.h:239
Lorene::Star_rot::p_angu_mom
double * p_angu_mom
Angular momentum.
Definition: star_rot.h:232
Lorene::Star_rot::ssjm1_nuf
Scalar ssjm1_nuf
Effective source at the previous step for the resolution of the Poisson equation for nuf by means of ...
Definition: star_rot.h:192
Lorene::Star_rot::p_area
double * p_area
Integrated surface area.
Definition: star_rot.h:238
Lorene::Star::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition: star.C:397
Lorene::Star::ent
Scalar ent
Log-enthalpy.
Definition: star.h:190
Lorene::Eos_poly::get_kap
double get_kap() const
Returns the pressure coefficient (cf.
Definition: eos_poly.C:260
Lorene::Star_rot::del_hydro_euler
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer.
Definition: star_rot.C:345
Lorene::Star_rot::ssjm1_nuq
Scalar ssjm1_nuq
Effective source at the previous step for the resolution of the Poisson equation for nuq by means of ...
Definition: star_rot.h:198
Lorene::Star_rot::p_aplat
double * p_aplat
Flatening r_pole/r_eq.
Definition: star_rot.h:237
Lorene::Tbl
Basic array class.
Definition: tbl.h:161
Lorene::Star_rot::p_z_pole
double * p_z_pole
Redshift factor at North pole.
Definition: star_rot.h:241
Lorene::Eos_poly
Polytropic equation of state (relativistic case).
Definition: eos.h:757
Lorene::Map::flat_met_cart
const Metric_flat & flat_met_cart() const
Returns the flat metric associated with the Cartesian coordinates and with components expressed in th...
Definition: map.C:331
Lorene::Mg3d::get_nzone
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Lorene::Star
Base class for stars.
Definition: star.h:175
Lorene::Star_rot::area
virtual double area() const
Integrated surface area in .
Definition: star_rot_global.C:366
Lorene::Eos_poly::get_gam
double get_gam() const
Returns the adiabatic index (cf. Eq. (3))
Definition: eos_poly.C:256
Lorene::Mg3d::get_nr
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
Lorene::Star_rot::display_poly
virtual void display_poly(ostream &) const
Display in polytropic units.
Definition: star_rot.C:649
Lorene::fread_be
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
Lorene::Star_rot::p_grv3
double * p_grv3
Error on the virial identity GRV3.
Definition: star_rot.h:235
Lorene::Star_rot::set_der_0x0
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition: star_rot.C:322
Lorene::Map::get_bvect_cart
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
Lorene::Star_rot::operator>>
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition: star_rot.C:428
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Lorene::Star::set_der_0x0
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition: star.C:316
Lorene::Scalar::set_dzpuis
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition: scalar.C:808
Lorene::Vector
Tensor field of valence 1.
Definition: vector.h:188
Lorene::Tensor::derive_con
const Tensor & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of this with respect to some metric , by raising the last inde...
Definition: tensor.C:1014
Lorene::Map::x
Coord x
x coordinate centered on the grid
Definition: map.h:726
Lorene::Star_rot::p_grv2
double * p_grv2
Error on the virial identity GRV2.
Definition: star_rot.h:234
Lorene::Star_rot::r_circ
virtual double r_circ() const
Circumferential radius.
Definition: star_rot_global.C:342
Lorene::Scalar::derive_con
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition: scalar_deriv.C:402
Lorene::Star::operator=
void operator=(const Star &)
Assignment to another Star.
Definition: star.C:351
Lorene::Star_rot::f_isco
virtual double f_isco() const
Orbital frequency at the innermost stable circular orbit (ISCO).
Definition: star_rot_isco.C:255
Lorene::Scalar::set_domain
Tbl & set_domain(int l)
Read/write of the value in a given domain.
Definition: scalar.h:615
Lorene::Star_rot::mean_radius
virtual double mean_radius() const
Mean star radius from the area .
Definition: star_rot_global.C:416
Lorene::Scalar::std_spectral_base
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
Definition: scalar.C:784
Lorene::Star_rot::aplat
virtual double aplat() const
Flatening r_pole/r_eq.
Definition: star_rot_global.C:426
Lorene::Star_rot::p_r_circ
double * p_r_circ
Circumferential radius.
Definition: star_rot.h:236
Lorene::Star_rot::ak_car
Scalar ak_car
Scalar .
Definition: star_rot.h:186
Lorene::fwrite_be
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
Lorene::Star_rot::z_eqb
virtual double z_eqb() const
Backward redshift factor at equator.
Definition: star_rot_global.C:471
Lorene::Star_rot::grv3
virtual double grv3(ostream *ost=0x0) const
Error on the virial identity GRV3.
Definition: star_rot_global.C:235
Lorene::Star_rot::w_shift
Vector w_shift
Vector used in the decomposition of shift , following Shibata's prescription [Prog.
Definition: star_rot.h:150
Lorene::Star::eos
const Eos & eos
Equation of state of the stellar matter.
Definition: star.h:185
Lorene::Tbl::get_taille
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:397
Lorene::Map::y
Coord y
y coordinate centered on the grid
Definition: map.h:727
Lorene::Star_rot::p_z_eqb
double * p_z_eqb
Backward redshift factor at equator.
Definition: star_rot.h:240
Lorene::Star_rot::uuu
Scalar uuu
Norm of u_euler.
Definition: star_rot.h:121
Lorene::Scalar::val_grid_point
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition: scalar.h:637
Lorene::Map::cosp
Coord cosp
Definition: map.h:724
Lorene::Tensor::set_triad
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition: tensor.C:519
Lorene::Map
Base class for coordinate mappings.
Definition: map.h:670
Lorene::Star_rot::b_car
Scalar b_car
Square of the metric factor B.
Definition: star_rot.h:110
Lorene::Star_rot::sauve
virtual void sauve(FILE *) const
Save in a file.
Definition: star_rot.C:401
Lorene::Tensor::sauve
virtual void sauve(FILE *) const
Save in a binary file.
Definition: tensor.C:906
Lorene::Star_rot::a_car
Scalar a_car
Square of the metric factor A.
Definition: star_rot.h:104
Lorene::Star_rot::dzeta
Scalar dzeta
Metric potential .
Definition: star_rot.h:134
Lorene::Star_rot::tggg
Scalar tggg
Metric potential .
Definition: star_rot.h:137
Lorene::Star_rot::fait_nphi
void fait_nphi()
Computes tnphi and nphi from the Cartesian components of the shift, stored in shift .
Definition: star_rot.C:748
Lorene::Star_rot::ssjm1_wshift
Vector ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for .
Definition: star_rot.h:225
Lorene::contract
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Definition: tenseur_operateur.C:279
Lorene::Star::beta
Vector beta
Shift vector.
Definition: star.h:228
Lorene::Vector::std_spectral_base
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: vector.C:316