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
29char 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// --------------------
136namespace Lorene {
137Etoile_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 ;
170 nphi = 0 ;
171 tnphi = 0 ;
172 nuf = 0 ;
173 nuq = 0 ;
174 tggg = 0 ;
175
177 for (int i=0; i<3; i++) {
178 w_shift.set(i) = 0 ;
179 }
180
182 khi_shift.set() = 0 ;
183
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// -----------------------
239Etoile_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
336
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 ;
428 ssjm1_khi = et.ssjm1_khi ;
430
431 del_deriv() ; // Deletes all derived quantities
432
433}
434
435 //--------------//
436 // Outputs //
437 //--------------//
438
439// Save in a file
440// --------------
441void 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
466ostream& 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
630void 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
681
682 return omega ;
683
684}
685
686
687// display_poly
688// ------------
689
690void 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
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
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 // --------------------
794
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).
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
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
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
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
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tenseur.h:704
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
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.