95 #include "star_rot_dirac.h" 97 #include "utilitaires.h" 102 double fact_omega,
int ,
const Tbl& ent_limit,
103 const Itbl& icontrol,
const Tbl& control,
104 double mbar_wanted,
double aexp_mass,
Tbl& diff){
113 char display_bold[]=
"x[1m" ; display_bold[0] = 27 ;
114 char display_normal[] =
"x[0m" ; display_normal[0] = 27 ;
126 assert( ( type_t == SYM) || (type_t == NONSYM) ) ;
128 int i_b = mg->
get_nr(l_b) - 1 ;
129 int j_b = (type_t == SYM ? mg->
get_nt(l_b) - 1 : mg->
get_nt(l_b)/2 ) ;
133 double ent_b = ent_limit(
nzet-1) ;
138 int mer_max = icontrol(0) ;
139 int mer_rot = icontrol(1) ;
140 int mer_change_omega = icontrol(2) ;
141 int mer_fix_omega = icontrol(3) ;
142 int mer_mass = icontrol(4) ;
143 int delta_mer_kep = icontrol(5) ;
144 int mer_hij = icontrol(6) ;
147 if (mer_change_omega < mer_rot) {
148 cout <<
"Star_rot_Dirac::equilibrium: mer_change_omega < mer_rot !" 150 cout <<
" mer_change_omega = " << mer_change_omega << endl ;
151 cout <<
" mer_rot = " << mer_rot << endl ;
154 if (mer_fix_omega < mer_change_omega) {
155 cout <<
"Star_rot_Dirac::equilibrium: mer_fix_omega < mer_change_omega !" 157 cout <<
" mer_fix_omega = " << mer_fix_omega << endl ;
158 cout <<
" mer_change_omega = " << mer_change_omega << endl ;
164 bool change_ent = true ;
167 mer_mass =
abs(mer_mass) ;
171 double precis = control(0) ;
172 double omega_ini = control(1) ;
173 double relax = control(2) ;
174 double relax_prev = double(1) - relax ;
180 double& diff_ent = diff.
set(0) ;
190 double accrois_omega = (omega0 - omega_ini) /
191 double(mer_fix_omega - mer_change_omega) ;
211 ofstream fichconv(
"convergence.d") ;
212 fichconv <<
"# diff_ent GRV2 max_triax vit_triax" << endl ;
214 ofstream fichfreq(
"frequency.d") ;
215 fichfreq <<
"# f [Hz]" << endl ;
217 ofstream fichevol(
"evolution.d") ;
219 "# |dH/dr_eq/dH/dr_pole| r_pole/r_eq ent_c" 223 double err_grv2 = 1 ;
231 for(
int mer=0 ; (diff_ent > precis) && (mer<mer_max) ; mer++) {
233 cout <<
"-----------------------------------------------" << endl ;
234 cout <<
"step: " << mer << endl ;
235 cout <<
"diff_ent = " << display_bold << diff_ent << display_normal
237 cout <<
"err_grv2 = " << err_grv2 << endl ;
244 if (mer >= mer_rot) {
246 if (mer < mer_change_omega) {
250 if (mer <= mer_fix_omega) {
251 omega = omega_ini + accrois_omega *
252 (mer - mer_change_omega) ;
288 if (mer > mer_fix_omega + delta_mer_kep) {
290 omega *= fact_omega ;
293 bool omega_trop_grand = false ;
300 bool superlum = true ;
329 for (
int l=0; l<
nzet; l++) {
330 for (
int i=0; i<mg->
get_nr(l); i++) {
335 cout <<
"U > c for l, i : " << l <<
" " << i
336 <<
" U = " <<
sqrt( u2 ) << endl ;
341 cout <<
"**** VELOCITY OF LIGHT REACHED ****" << endl ;
342 omega /= fact_omega ;
343 cout <<
"New rotation frequency : " 344 <<
omega/(2.*M_PI) * f_unit <<
" Hz" << endl ;
345 omega_trop_grand = true ;
373 double mlngamma_c = 0 ;
377 double alpha_r2 = ( ent_c - ent_b + mlngamma_c - mlngamma_b
378 + ln_q_c - ln_q_b) / ( ln_f_b - ln_f_c ) ;
380 alpha_r =
sqrt(alpha_r2) ;
382 cout <<
"alpha_r = " << alpha_r << endl ;
391 logn = alpha_r2 * ln_f_new + ln_q_new ;
398 ent = (ent_c + logn_c + mlngamma_c) -
logn - mlngamma ;
407 for (
int l=0; l<
nzet; l++) {
408 int imax = mg->
get_nr(l) - 1 ;
409 if (l == l_b) imax-- ;
410 for (
int i=0; i<imax; i++) {
413 cout <<
"ent < 0 for l, i : " << l <<
" " << i
420 cout <<
"**** KEPLERIAN VELOCITY REACHED ****" << endl ;
421 omega /= fact_omega ;
422 cout <<
"New rotation frequency : " 423 <<
omega/(2.*M_PI) * f_unit <<
" Hz" << endl ;
424 omega_trop_grand = true ;
429 if ( omega_trop_grand ) {
431 fact_omega =
sqrt( fact_omega ) ;
432 cout <<
"**** New fact_omega : " << fact_omega << endl ;
464 hij_new.set_etat_zero() ;
484 logn = relax *
logn + relax_prev * logn_prev ;
486 qqq = relax *
qqq + relax_prev * qqq_prev ;
500 fichfreq <<
" " <<
omega / (2*M_PI) * f_unit ;
501 fichevol <<
" " << ent_c ;
508 if (mer > mer_mass) {
511 if (mbar_wanted > 0.) {
512 xx =
mass_b() / mbar_wanted - 1. ;
513 cout <<
"Discrep. baryon mass <-> wanted bar. mass : " << xx
517 xx =
mass_g() / fabs(mbar_wanted) - 1. ;
518 cout <<
"Discrep. grav. mass <-> wanted grav. mass : " << xx
521 double xprog = ( mer > 2*mer_mass) ? 1. :
522 double(mer-mer_mass)/double(mer_mass) ;
524 double ax = .5 * ( 2. + xx ) / (1. + xx ) ;
525 double fact =
pow(ax, aexp_mass) ;
526 cout <<
" xprog, xx, ax, fact : " << xprog <<
" " <<
527 xx <<
" " << ax <<
" " << fact << endl ;
533 if (mer%4 == 0)
omega *= fact ;
544 diff_ent = diff_ent_tbl(0) ;
545 for (
int l=1; l<
nzet; l++) {
546 diff_ent += diff_ent_tbl(l) ;
550 fichconv <<
" " <<
log10( fabs(diff_ent) + 1.e-16 ) ;
551 fichconv <<
" " <<
log10( fabs(err_grv2) + 1.e-16 ) ;
void solve_logn_f(Scalar &ln_f_new) const
Solution of the two scalar Poisson equations for rotating stars in Dirac gauge.
Cmp log(const Cmp &)
Neperian logarithm.
double omega
Rotation angular velocity ([f_unit] )
Map & mp
Mapping associated with the star.
Cmp sqrt(const Cmp &)
Square root.
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
double & set(int i)
Read/write of a particular element (index i) (1D case)
Tensor field of valence 0 (or component of a tensorial field).
Sym_tensor_trans hh
is defined by .
int relat_type
Relativistic flag.
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Basic integer array class.
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Tensor field of valence 1.
void update_metric()
Computes metric quantities from known potentials.
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
int nzet
Number of domains of *mp occupied by the star.
Scalar gam_euler
Lorentz factor between the fluid and Eulerian observers.
virtual double mass_g() const
Gravitational mass.
virtual void homothetie(double lambda)=0
Sets a new radial scale.
virtual double grv2() const
Error on the virial identity GRV2.
int get_nzone() const
Returns the number of domains.
Vector u_euler
Fluid 3-velocity with respect to the Eulerian observer.
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
Cmp pow(const Cmp &, int)
Power .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Scalar logn
Logarithm of the lapse N .
virtual double mass_b() const
Baryonic mass.
Transverse symmetric tensors of rank 2.
void equilibrium(double ent_c, double omega0, double fact_omega, int nzadapt, const Tbl &ent_limit, const Itbl &icontrol, const Tbl &control, double mbar_wanted, double aexp_mass, Tbl &diff)
Computes an equilibrium configuration.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
void solve_hij(Sym_tensor_trans &hij_new) const
Solution of the tensor Poisson equation for rotating stars in Dirac gauge.
Scalar nn
Lapse function N .
Cmp log10(const Cmp &)
Basis 10 logarithm.
Cmp abs(const Cmp &)
Absolute value.
void solve_shift(Vector &shift_new) const
Solution of the shift equation for rotating stars in Dirac gauge.
void solve_logn_q(Scalar &ln_q_new) const
Solution of the two scalar Poisson equations for rotating stars in Dirac gauge.
void solve_qqq(Scalar &q_new) const
Solution of the two scalar Poisson equations for rotating stars in Dirac gauge.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Scalar & set(int)
Read/write access to a component.
void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
const Metric_flat & flat
flat metric (spherical components)
void annule_hard()
Sets the Tbl to zero in a hard way.