24 char binhor_coal_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Bin_hor/binhor_coal.C,v 1.15 2014/10/13 08:52:42 j_novak Exp $" ;
92 #include "graphique.h" 97 double lim_nn,
int bound_psi) {
110 cout <<
"Static black holes : " << endl ;
121 for (
int i=1 ; i<nz ; i++)
122 if (diff(i) > erreur)
125 cout <<
"Step : " << conte <<
" Difference : " << erreur << endl ;
134 int nb_it,
int bound_nn,
double lim_nn,
135 int bound_psi,
int bound_beta,
double omega_eff,
137 ostream& fich_iteration, ostream& fich_correction,
138 ostream& fich_viriel, ostream& fich_kss,
139 int step,
int search_mass,
double mass_irr,
144 double precis = 1e-7 ;
147 cout <<
"OMEGA INCREASED STEP BY STEP." << endl ;
149 double inc_homme = (angu_vel - homme)/nb_ome ;
150 for (
int pas = 0 ; pas <nb_ome ; pas ++) {
153 if (omega_eff == alpha*homme ) verif = true ;
158 omega_eff = alpha*homme ;
161 solve_shift (precis, relax, bound_beta, omega_eff) ;
168 if (search_mass == 1 && step >= 30) {
171 double error_m = (mass_irr - mass_area) / mass_irr ;
172 double scaling_r =
pow((2-error_m)/(2-2*error_m), 1.) ;
188 cout <<
"Angular momentum computed at the horizon : " <<
ang_mom_hor()
193 for (
int i=1 ; i<nz ; i++)
194 if (diff(i) > erreur)
207 for (
int k=0 ; k<nnp ; k++)
208 for (
int j=0 ; j<nnt ; j++){
216 fich_iteration << step <<
" " <<
log10(erreur) <<
" " << homme << endl ;
217 fich_correction << step <<
" " <<
log10(
hole1.
regul) <<
" " << homme << endl ;
219 fich_viriel << step <<
" " <<
viriel() <<
" " << homme <<
" " <<
hole1.
omega_hor() - alpha*homme <<
" " << omega_eff << endl ;
220 fich_kss << step <<
" " << max_kss <<
" " << min_kss << endl ;
223 cout <<
"STEP : " << step <<
" DIFFERENCE : " << erreur << endl ;
230 cout <<
"OMEGA FIXED" << endl ;
233 for (
int pas = 0 ; pas <nb_it ; pas ++) {
237 solve_shift (precis, relax, bound_beta, omega_eff) ;
244 if (search_mass == 1 && step >= 30) {
247 double error_m = (mass_irr - mass_area) / mass_irr ;
248 double scaling_r =
pow((2-error_m)/(2-2*error_m), 1.) ;
258 for (
int i=1 ; i<nz ; i++)
259 if (diff(i) > erreur)
272 for (
int k=0 ; k<nnp ; k++)
273 for (
int j=0 ; j<nnt ; j++){
282 fich_iteration << step <<
" " <<
log10(erreur) <<
" " << homme << endl ;
283 fich_correction << step <<
" " <<
log10(
hole1.
regul) <<
" " << homme << endl ;
285 fich_viriel << step <<
" " <<
viriel() <<
" " << homme <<
" " <<
hole1.
omega_hor() - alpha*homme <<
" " << omega_eff << endl ;
286 fich_kss << step <<
" " << max_kss <<
" " << min_kss << endl ;
289 cout <<
"STEP : " << step <<
" DIFFERENCE : " << erreur << endl ;
294 fich_iteration <<
"#----------------------------" << endl ;
295 fich_correction <<
"#-----------------------------" << endl ;
296 fich_viriel <<
"#------------------------------" << endl ;
void solve_shift(double precis, double relax, int bound_beta, double omega_eff)
Solves the equation for the shift, using the Oohara-Nakarmure scheme : The fields are the total value...
const Metric & get_gam() const
metric
double omega_hor() const
Orbital velocity.
double ang_mom_hor() const
Calculates the angular momentum of the black hole using the formula at the horizon.
Metric tgam
3 metric tilde
void extrinsic_curvature()
Calculation of the extrinsic curvature tensor.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Cmp sqrt(const Cmp &)
Square root.
double area_hor() const
Area of the horizon.
const Sym_tensor & get_k_dd() const
k_dd
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Tensor field of valence 0 (or component of a tensorial field).
Single_hor hole1
Black hole one.
virtual const Vector & radial_vect() const
Returns the radial vector normal to a spherical slicing and pointing toward spatial infinity...
void set_statiques(double precis, double relax, int bound_nn, double lim_nn, int bound_psi)
Initialize the systeme to Misner Lindquist solution, that is solving for N and in the case ...
double viriel() const
Computes the viriel error, that is the difference between the ADM and the Komar masses, calculated by the asymptotic behaviours of respectively and N .
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
void set_omega(double ome)
Sets the orbital velocity to ome.
double coal(double ang_vel, double relax, int nb_om, int nb_it, int bound_nn, double lim_nn, int bound_psi, int bound_beta, double omega_eff, double alpha, ostream &fich_iteration, ostream &fich_correction, ostream &fich_viriel, ostream &fich_kss, int step, int search_mass, double mass_irr, const int sortie=0)
Solves the equation for a particular angular velocity, the systeme being initialized to Misner-Lindqu...
void init_bin_hor()
Initialisation of the system.
double regul
Intensity of the correction on the shift vector.
void solve_psi(double precis, double relax, int bound_psi)
Solves the equation for the conformal factor : The fields are the total values excpet those with subs...
Vector beta_auto
Shift function .
int get_nzone() const
Returns the number of domains.
Cmp pow(const Cmp &, int)
Power .
Scalar n_auto
Lapse function .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Single_hor hole2
Black hole two.
void solve_lapse(double precis, double relax, int bound_nn, double lim_nn)
Solves the equation for the lapse : The fields are the total values except those with subscript ...
Map_af & mp
Affine mapping.
Cmp log10(const Cmp &)
Basis 10 logarithm.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
void homothetie_interne(double lambda)
Sets a new radial scale at the bondary between the nucleus and the first shell.
void init_met_trK()
Sets the 3-metric tilde to the flat metric and gamt_point, trK and trK_point to zero.
double radius
Radius of the horizon in LORENE's units.
Metric_flat ff
3 metric flat
double get_omega() const
Returns the angular velocity.
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).