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