23char bin_ns_bh_coal_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Bin_ns_bh/bin_ns_bh_coal.C,v 1.12 2014/10/13 08:52:43 j_novak Exp $" ;
82void Bin_ns_bh::coal (
double precis,
double relax,
int itemax_equil,
83 int itemax_mp_et,
double ent_c_init,
double seuil_masses,
84 double dist,
double m1,
double m2,
double spin_cible,
85 double scale_ome_local,
const int sortie,
int bound_nn,
98 double scale_linear = (mass_ns+mass_bh)/2.*distance*
omega ;
100 char name_iteration[40] ;
101 char name_correction[40] ;
102 char name_viriel[40] ;
104 char name_linear[40] ;
106 char name_error_m1[40] ;
107 char name_error_m2[40] ;
112 char name_ome_local[40] ;
114 sprintf(name_iteration,
"ite.dat") ;
115 sprintf(name_correction,
"cor.dat") ;
116 sprintf(name_viriel,
"vir.dat") ;
117 sprintf(name_ome,
"ome.dat") ;
118 sprintf(name_linear,
"linear.dat") ;
119 sprintf(name_axe,
"axe.dat") ;
120 sprintf(name_error_m1,
"error_m_bh.dat") ;
121 sprintf(name_error_m2,
"error_m_ns.dat") ;
122 sprintf(name_hor,
"hor.dat") ;
123 sprintf(name_entc,
"entc.dat") ;
124 sprintf(name_dist,
"distance.dat") ;
125 sprintf(name_spin,
"spin.dat") ;
126 sprintf(name_ome_local,
"ome_local.dat") ;
128 ofstream fiche_iteration(name_iteration) ;
129 fiche_iteration.precision(8) ;
131 ofstream fiche_correction(name_correction) ;
132 fiche_correction.precision(8) ;
134 ofstream fiche_viriel(name_viriel) ;
135 fiche_viriel.precision(8) ;
137 ofstream fiche_ome(name_ome) ;
138 fiche_ome.precision(8) ;
140 ofstream fiche_linear(name_linear) ;
141 fiche_linear.precision(8) ;
143 ofstream fiche_axe(name_axe) ;
144 fiche_axe.precision(8) ;
146 ofstream fiche_error_m1 (name_error_m1) ;
147 fiche_error_m1.precision(8) ;
149 ofstream fiche_error_m2 (name_error_m2) ;
150 fiche_error_m2.precision(8) ;
152 ofstream fiche_hor (name_hor) ;
153 fiche_hor.precision(8) ;
155 ofstream fiche_entc (name_entc) ;
156 fiche_entc.precision(8) ;
158 ofstream fiche_dist (name_dist) ;
159 fiche_dist.precision(8) ;
161 ofstream fiche_spin (name_spin) ;
162 fiche_spin.precision(8) ;
164 ofstream fiche_ome_local (name_ome_local) ;
165 fiche_ome_local.precision(8) ;
168 bool search_masses = false ;
169 double ent_c = ent_c_init ;
171 Cmp shift_bh_old (
hole.
mp) ;
172 Cmp shift_ns_old (
star.
mp) ;
187 spin =
hole.local_momentum() ;
190 fiche_spin << conte <<
" " << spin/m1/m1 << endl ;
193 double conv_spin = fabs(spin-spin_old) ;
194 double error_spin = spin - spin_cible ;
195 double rel_diff_spin = (spin_cible==0) ? fabs(error_spin) :
196 fabs(error_spin)/spin_cible ;
197 if ((conv_spin*2<rel_diff_spin) && (search_masses) &&
199 double func = scale_ome_local*
log((2+error_spin)/(2+2*error_spin));
203 old_mass_ns = mass_ns ;
220 diff.set_etat_qcq() ;
225 relax, itemax_mp_et, relax, diff) ;
229 hole.equilibrium (
star, precis, relax, bound_nn, lim_nn) ;
230 cout <<
"Apres equilibrium" << endl ;
233 cout <<
"Apres star::update_metric" << endl ;
236 cout <<
"Apres star::update_metric_der_comp" << endl ;
238 cout <<
"Apres Bin_ns_bh::fait_tkij" << endl ;
242 for (
int i=1 ; i<nz_bh ; i++)
243 if (diff_bh(i) > erreur)
244 erreur = diff_bh(i) ;
246 for (
int i=0 ; i<nz_ns ; i++)
247 if (diff_ns(i) > erreur)
248 erreur = diff_ns(i) ;
250 if (erreur<seuil_masses)
251 search_masses = true ;
255 cout <<
"Avant viriel" << endl ;
256 double error_viriel = viriel() ;
257 cout <<
"Apres viriel" << endl ;
258 double error_linear = linear_momentum_systeme_inf()(1)/scale_linear ;
259 cout <<
"Apres linear" << endl ;
261 cout <<
"Apres Mbh" << endl ;
262 double error_m2 = 1 - mass_ns/m2 ;
263 cout <<
"Apres Mns" << endl ;
266 fiche_iteration << conte <<
" " << erreur << endl ;
267 fiche_correction << conte <<
" " <<
hole.
regul << endl ;
268 fiche_viriel << conte <<
" " << error_viriel << endl ;
269 fiche_linear << conte <<
" " << error_linear << endl ;
270 fiche_error_m1 << conte <<
" " << error_m1 << endl ;
271 fiche_error_m2 << conte <<
" " << error_m2 << endl ;
272 fiche_hor << conte <<
" " <<
hole.
mp.
val_r(0, 1, 0,0) << endl ;
273 fiche_entc << conte <<
" " << ent_c << endl ;
274 fiche_dist << conte <<
" " << distance << endl ;
275 fiche_ome << conte <<
" " <<
omega << endl ;
276 fiche_axe << conte <<
" " << axe_pos << endl ;
280 double scaling_axe =
pow((2+error_linear)/(2+2*error_linear), 0.1) ;
281 axe_pos *= scaling_axe ;
286 double new_ome =
star.compute_angul() ;
294 double error_dist = (distance-dist)/dist ;
295 double scale_d =
pow((2+error_dist)/(2+2*error_dist), 0.2) ;
296 distance *= scale_d ;
302 double scaling_r =
pow((2-error_m1)/(2-2*error_m1), 0.25) ;
308 double convergence = fabs(mass_ns - old_mass_ns)/mass_ns ;
309 double rel_diff = fabs(error_m2) ;
310 if ((search_masses) && (convergence*2 < rel_diff)) {
311 double scaling_ent =
pow((2-error_m2)/(2-2*error_m2), 1) ;
312 ent_c *= scaling_ent ;
318 cout <<
"PAS TOTAL : " << conte <<
" DIFFERENCE : " << erreur << endl ;
328 fiche_iteration.close() ;
329 fiche_correction.close() ;
330 fiche_viriel.close() ;
332 fiche_linear.close() ;
334 fiche_error_m1.close() ;
335 fiche_error_m2.close() ;
340 fiche_ome_local.close() ;
double omega_local
local angular velocity
double regul
Intensity of the correction on the shift vector.
double masse_adm_seul() const
Calculates the ADM mass of the black hole using : .
Map_af & set_mp()
Read/write of the mapping.
double area() const
Computes the area of the throat.
double get_rot_state() const
Returns the state of rotation.
Map_af & mp
Affine mapping.
const Tenseur & get_shift_auto() const
Returns the part of generated by the hole.
double get_rayon() const
Returns the radius of the horizon.
void set_omega_local(double ome)
Sets the local angular velocity to ome .
void set_rayon(double ray)
Sets the radius of the horizon to ray .
void fait_tkij(int bound_nn=-1, double lim_nn=0)
Computation of the extrinsic curvature tensor for both {\tt star} and {\tt bhole}.
void set_omega(double)
Sets the orbital angular velocity [{\tt f_unit}].
double x_axe
Absolute X coordinate of the rotation axis.
Et_bin_nsbh star
The neutron star.
Bhole hole
The black hole.
double omega
Angular velocity with respect to an asymptotically inertial observer.
void update_metric(const Bhole &comp)
Computes metric coefficients from known potentials, when the companion is a black hole.
virtual void equilibrium_nsbh(double ent_c, int mermax, int mermax_poisson, double relax_poisson, int mermax_potvit, double relax_potvit, double thres_adapt, const Tbl &fact, Tbl &diff)
Computes an equilibrium configuration in a NS-BH binary system.
virtual void kinematics(double omega, double x_axe)
Computes the quantities bsn and pot_centri .
void update_metric_der_comp(const Bhole &comp)
Computes the derivative of metric functions related to the companion black hole.
virtual double mass_b() const
Baryon mass.
const Tenseur & get_shift_auto() const
Returns the part of the shift vector generated principaly by the star.
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
virtual double mass_g() const
Gravitational mass.
void fait_d_psi()
Computes the gradient of the total velocity potential .
Map & mp
Mapping associated with the star.
Map & set_mp()
Read/write of the mapping.
void homothetie_interne(double lambda)
Sets a new radial scale at the bondary between the nucleus and the first shell.
virtual double val_r(int l, double xi, double theta, double pphi) const
Returns the value of the radial coordinate r for a given in a given domain.
void set_ori(double xa0, double ya0, double za0)
Sets a new origin.
double get_ori_x() const
Returns the x coordinate of the origin.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
int get_nzone() const
Returns the number of domains.
int get_etat() const
Returns the logical state.
Cmp sqrt(const Cmp &)
Square root.
Cmp pow(const Cmp &, int)
Power .
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Cmp log(const Cmp &)
Neperian logarithm.
Standard units of space, time and mass.