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
28char 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// --------------------
85namespace Lorene {
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
135
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// -----------------------
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
243
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
252
255
258
261
264
267
268 // All other fields are initialized to zero :
269 // ----------------------------------------
270 a_car = 0 ;
271 bbb = 0 ;
272 b_car = 0 ;
273 uuu = 0 ;
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
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
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// --------------------------------
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
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 ;
387 ssjm1_khi = et.ssjm1_khi ;
389
390 del_deriv() ; // Deletes all derived quantities
391
392}
393
394
395 //--------------//
396 // Outputs //
397 //--------------//
398
399// Save in a file
400// --------------
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) ;
414
421
422
423}
424
425// Printing
426// --------
427
429
430 using namespace Unites ;
431
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
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
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
640double Star_rot::get_omega_c() const {
641
642 return omega ;
643
644}
645
646// display_poly
647// ------------
648
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
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) ;
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}
Polytropic equation of state (relativistic case).
Definition eos.h:757
Equation of state base class.
Definition eos.h:190
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Base class for coordinate mappings.
Definition map.h:670
Coord cosp
Definition map.h:724
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
Coord y
y coordinate centered on the grid
Definition map.h:727
Coord sint
Definition map.h:721
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...
Coord sinp
Definition map.h:723
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
Coord z
z coordinate centered on the grid
Definition map.h:728
Coord x
x coordinate centered on the grid
Definition map.h:726
Coord cost
Definition map.h:722
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
virtual void sauve(FILE *) const
Save in a file.
Definition scalar.C:686
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
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
void div_rsint()
Division by everywhere; dzpuis is not changed.
const Tbl & domain(int l) const
Read-only of the value in a given domain.
Definition scalar.h:625
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Class for isolated rotating stars.
Definition star_rot.h:86
double * p_angu_mom
Angular momentum.
Definition star_rot.h:232
virtual void display_poly(ostream &) const
Display in polytropic units.
Definition star_rot.C:649
virtual double mean_radius() const
Mean star radius from the area .
double * p_r_isco
Circumferential radius of the ISCO.
Definition star_rot.h:243
double * p_aplat
Flatening r_pole/r_eq.
Definition star_rot.h:237
Star_rot(Map &mp_i, int nzet_i, bool relat, const Eos &eos_i)
Standard constructor.
Definition star_rot.C:86
double * p_mom_quad
Quadrupole moment.
Definition star_rot.h:242
Scalar ssjm1_tggg
Effective source at the previous step for the resolution of the Poisson equation for tggg .
Definition star_rot.h:208
Sym_tensor tkij
Tensor related to the extrinsic curvature tensor by .
Definition star_rot.h:167
virtual double mom_quad() const
Quadrupole moment.
double * p_z_eqf
Forward redshift factor at equator.
Definition star_rot.h:239
Scalar tggg
Metric potential .
Definition star_rot.h:137
Scalar tnphi
Component of the shift vector.
Definition star_rot.h:118
void fait_shift()
Computes shift from w_shift and khi_shift according to Shibata's prescription [Prog.
Definition star_rot.C:703
virtual ~Star_rot()
Destructor.
Definition star_rot.C:287
virtual double z_pole() const
Redshift factor at North pole.
virtual double r_isco(ostream *ost=0x0) const
Circumferential radius of the innermost stable circular orbit (ISCO).
virtual double z_eqb() const
Backward redshift factor at equator.
double unsurc2
: unsurc2=1 for a relativistic star, 0 for a Newtonian one.
Definition star_rot.h:99
Vector w_shift
Vector used in the decomposition of shift , following Shibata's prescription [Prog.
Definition star_rot.h:150
Scalar b_car
Square of the metric factor B.
Definition star_rot.h:110
void operator=(const Star_rot &)
Assignment to another Star_rot.
Definition star_rot.C:359
double * p_lspec_isco
Specific angular momentum of a particle on the ISCO.
Definition star_rot.h:248
Scalar bbb
Metric factor B.
Definition star_rot.h:107
virtual double mass_b() const
Baryon mass.
virtual double get_omega_c() const
Returns the central value of the rotation angular velocity ([f_unit] )
Definition star_rot.C:640
double omega
Rotation angular velocity ([f_unit] )
Definition star_rot.h:101
Scalar ak_car
Scalar .
Definition star_rot.h:186
virtual double r_circ() const
Circumferential radius.
Scalar uuu
Norm of u_euler.
Definition star_rot.h:121
virtual double z_eqf() const
Forward redshift factor at equator.
Scalar nphi
Metric coefficient .
Definition star_rot.h:113
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition star_rot.C:322
double * p_f_isco
Orbital frequency of the ISCO.
Definition star_rot.h:244
virtual double f_isco() const
Orbital frequency at the innermost stable circular orbit (ISCO).
virtual double mass_g() const
Gravitational mass.
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition star_rot.C:428
virtual void sauve(FILE *) const
Save in a file.
Definition star_rot.C:401
double * p_area
Integrated surface area.
Definition star_rot.h:238
virtual double aplat() const
Flatening r_pole/r_eq.
double * p_grv3
Error on the virial identity GRV3.
Definition star_rot.h:235
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
double * p_espec_isco
Specific energy of a particle on the ISCO.
Definition star_rot.h:246
double * p_tsw
Ratio T/W.
Definition star_rot.h:233
Scalar nuf
Part of the Metric potential = logn generated by the matter terms.
Definition star_rot.h:126
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one.
Definition star_rot.h:94
double * p_f_eq
Orbital frequency at the equator.
Definition star_rot.h:249
Scalar ssjm1_dzeta
Effective source at the previous step for the resolution of the Poisson equation for dzeta .
Definition star_rot.h:203
virtual double tsw() const
Ratio T/W.
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer.
Definition star_rot.C:345
virtual double grv3(ostream *ost=0x0) const
Error on the virial identity GRV3.
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
Scalar dzeta
Metric potential .
Definition star_rot.h:134
virtual double area() const
Integrated surface area in .
Scalar a_car
Square of the metric factor A.
Definition star_rot.h:104
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
Scalar nuq
Part of the Metric potential = logn generated by the quadratic terms.
Definition star_rot.h:131
double * p_grv2
Error on the virial identity GRV2.
Definition star_rot.h:234
virtual void del_deriv() const
Deletes all the derived quantities.
Definition star_rot.C:297
Vector ssjm1_wshift
Effective source at the previous step for the resolution of the vector Poisson equation for .
Definition star_rot.h:225
double * p_r_circ
Circumferential radius.
Definition star_rot.h:236
virtual double grv2() const
Error on the virial identity GRV2.
virtual double angu_mom() const
Angular momentum.
double * p_z_eqb
Backward redshift factor at equator.
Definition star_rot.h:240
virtual void partial_display(ostream &) const
Printing of some informations, excluding all global quantities.
Definition star_rot.C:590
double * p_z_pole
Redshift factor at North pole.
Definition star_rot.h:241
Scalar khi_shift
Scalar used in the decomposition of shift , following Shibata's prescription [Prog.
Definition star_rot.h:160
void fait_nphi()
Computes tnphi and nphi from the Cartesian components of the shift, stored in shift .
Definition star_rot.C:748
Base class for stars.
Definition star.h:175
Scalar ener
Total energy density in the fluid frame.
Definition star.h:193
Scalar logn
Logarithm of the lapse N .
Definition star.h:222
Scalar nn
Lapse function N .
Definition star.h:225
virtual void del_hydro_euler()
Sets to ETATNONDEF (undefined state) the hydrodynamical quantities relative to the Eulerian observer.
Definition star.C:330
const Eos & eos
Equation of state of the stellar matter.
Definition star.h:185
Scalar nbar
Baryon density in the fluid frame.
Definition star.h:192
virtual void del_deriv() const
Deletes all the derived quantities.
Definition star.C:298
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Definition star.C:417
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition star.C:316
double ray_eq() const
Coordinate radius at , [r_unit].
virtual void sauve(FILE *) const
Save in a file.
Definition star.C:397
Scalar press
Fluid pressure.
Definition star.h:194
Scalar ent
Log-enthalpy.
Definition star.h:190
Map & mp
Mapping associated with the star.
Definition star.h:180
void operator=(const Star &)
Assignment to another Star.
Definition star.C:351
int nzet
Number of domains of *mp occupied by the star.
Definition star.h:183
Vector beta
Shift vector.
Definition star.h:228
Basic array class.
Definition tbl.h:161
Tensor handling.
Definition tensor.h:288
Tensor field of valence 1.
Definition vector.h:188
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
virtual void sauve(FILE *) const
Save in a binary file.
Definition tensor.C:906
void set_triad(const Base_vect &new_triad)
Assigns a new vectorial basis (triad) of decomposition.
Definition tensor.C:519
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
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition tensor.C:497
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.