85 ifstream file(file_name) ;
87 cerr <<
"Problem in opening the file " << file_name << endl ;
93 int nrfile, nthetafile;
94 file >> nrfile >> nthetafile ;
98 if (
rHor<0. || nrfile<0. || nthetafile<0.){
99 cerr <<
"In scalarBH::scalarBH(Map,char*): " 100 <<
"Bad parameters rHor, nrfile, nthetafile" << endl;
107 cout <<
"nr, ntheta from file = " << nrfile <<
" " << nthetafile << endl;
108 cout <<
"rHor from file = " <<
rHor << endl;
110 cout <<
"Scalar field values provided" << endl;
112 cout <<
"Scalar field values not provided, put to zero" << endl;
114 cerr <<
"In scalarBH::scalarBH(Map,char*): " 115 <<
"Bad parameter isphi" << endl;
119 double* Xfile =
new double[nrfile * nthetafile] ;
120 double* thetafile =
new double[nrfile * nthetafile] ;
121 double* f0file =
new double[nrfile * nthetafile] ;
122 double* f1file =
new double[nrfile * nthetafile] ;
123 double* f2file =
new double[nrfile * nthetafile] ;
124 double* wwfile =
new double[nrfile * nthetafile] ;
125 double* sfieldfile =
new double[nrfile * nthetafile] ;
127 cout <<
"Reading metric data... ";
128 for (
int ii=0;ii<nrfile*nthetafile;ii++){
131 file >> thetafile[ii] ;
136 file >> sfieldfile[ii] ;
143 cout <<
"done" << endl;
146 double Xbefmax = Xfile[nrfile-2];
156 double delta_theta = 0.1*fabs(thetafile[nrfile] - thetafile[0]);
158 cout <<
"Starting interpolating Lorene grid... " ;
167 for (
int l=l_min; l<nz; l++) {
172 for (
int k=0; k<np; k++){
173 for (
int j=0; j<nt; j++){
174 double th0 = theta(l, k, j, 0);
175 for (
int i=0; i<nr; i++){
176 double r0 = rr(l, k, j, i);
178 if (r0 < __infinity){
179 x0 =
sqrt(r0*r0 - rH2);
190 double thc = thetafile[ith];
191 while (fabs(th0 - thc) > delta_theta){
194 thc = thetafile[ith];
197 double xc = Xfile[ir];
199 cerr <<
"In scalarBH::ScalarBH(): " 200 "r should be zero here" << endl;
210 }
else if(xx0 > Xbefmax && xx0< 1.){
223 double f0interp, f1interp, f2interp, winterp, sfieldinterp;
225 double xcinf = Xfile[ir-1], xcsup = Xfile[ir+1];
227 if (ir-3>0) {irext1=ir-2; irext2=ir-3;}
228 else if (ir+3<nrfile) {irext1=ir+2; irext2=ir+3;}
230 cerr <<
"scalarBH::scalarBH(): bad radial indice" << endl;
233 double xcext1 = Xfile[irext1], xcext2 = Xfile[irext2];
241 if (fabs(thc-th0)>delta_theta){
242 cerr <<
"scalarBH::ScalarBH(): theta problem in grid" << endl;
243 cerr <<
"Theta info: " << thc <<
" " << th0 << endl;
246 if (xx0 <= Xbefmax &&
247 (xx0 < xcinf || xx0 > xc || xx0 > xcsup)){
248 cerr <<
"scalarBH::ScalarBH(): rad problem in grid" << endl;
249 cerr <<
"Radial info: " << xcinf <<
" " << xx0 <<
" " 250 << xc <<
" " << xcsup << endl;
252 }
else if (xx0 > Xbefmax &&
253 (xx0 < xcinf || xx0 < xc || xx0 > xcsup)){
254 cerr <<
"scalarBH::ScalarBH(): special rad " 255 "problem in grid" << endl;
256 cerr <<
"Radial info: " << xcinf <<
" " << xx0 <<
" " 257 << xc <<
" " << xcsup << endl;
261 double polyrinf = (xx0-xc)*(xx0-xcsup)*(xx0-xcext1)*(xx0-xcext2)/((xcinf-xc)*(xcinf-xcsup)*(xcinf-xcext1)*(xcinf-xcext2));
262 double polyrmid = (xx0-xcinf)*(xx0-xcsup)*(xx0-xcext1)*(xx0-xcext2)/((xc-xcinf)*(xc-xcsup)*(xc-xcext1)*(xc-xcext2));
263 double polyrsup = (xx0-xcinf)*(xx0-xc)*(xx0-xcext1)*(xx0-xcext2)/((xcsup-xcinf)*(xcsup-xc)*(xcsup-xcext1)*(xcsup-xcext2));
264 double polyrext1 = (xx0-xcinf)*(xx0-xc)*(xx0-xcsup)*(xx0-xcext2)/((xcext1-xcinf)*(xcext1-xc)*(xcext1-xcsup)*(xcext1-xcext2));
265 double polyrext2 = (xx0-xcinf)*(xx0-xc)*(xx0-xcsup)*(xx0-xcext1)/((xcext2-xcinf)*(xcext2-xc)*(xcext2-xcsup)*(xcext2-xcext1));
268 double f0ext1 = f0file[irext1], f0ext2 = f0file[irext2],
269 f0inf = f0file[ir-1],
271 f0sup = f0file[ir+1], f1ext1=f1file[irext1],
272 f1ext2=f1file[irext2],
273 f1inf = f1file[ir-1], f1mid = f1file[ir],
274 f1sup = f1file[ir+1], f2ext1=f2file[irext1],
275 f2ext2=f2file[irext2],
276 f2inf = f2file[ir-1], f2mid = f2file[ir],
277 f2sup = f2file[ir+1], wext1=wwfile[irext1],
278 wext2=wwfile[irext2],
279 winf = wwfile[ir-1], wmid = wwfile[ir],
280 wsup = wwfile[ir+1], sfext1=sfieldfile[irext1],
281 sfext2=sfieldfile[irext2],
282 sfinf = sfieldfile[ir-1],
283 sfmid = sfieldfile[ir], sfsup = sfieldfile[ir+1];
290 f0interp = f0ext1*polyrext1 + f0ext2*polyrext2 + f0inf*polyrinf
291 + f0mid*polyrmid + f0sup*polyrsup;
292 f1interp = f1ext1*polyrext1 + f1ext2*polyrext2 + f1inf*polyrinf
293 + f1mid*polyrmid + f1sup*polyrsup;
294 f2interp = f2ext1*polyrext1 + f2ext2*polyrext2 + f2inf*polyrinf
295 + f2mid*polyrmid + f2sup*polyrsup;
296 winterp = wext1*polyrext1 + wext2*polyrext2 + winf*polyrinf
297 + wmid*polyrmid + wsup*polyrsup;
298 sfieldinterp = sfext1*polyrext1 + sfext2*polyrext2
300 + sfmid*polyrmid + sfsup*polyrsup;
307 f0interp = f0file[ir] ;
308 f1interp = f1file[ir] ;
309 f2interp = f2file[ir] ;
310 winterp = wwfile[ir] ;
311 sfieldinterp = sfieldfile[ir] ;
340 cout <<
"done." << endl;
345 cout <<
"Starting updating metric... " ;
347 cout <<
"done." << endl;
450 ost << endl <<
"Black hole with scalar hair (class ScalarBH) " << endl ;
504 void ScalarBH::update_metric() {
509 NN.set_domain(0) = 1;
517 gam.set(1,1) =
exp(2*
ff1)/NN ;
518 gam.
set(1,1).std_spectral_base() ;
521 gam.set(2,2) =
exp(2*
ff1);
522 gam.set(2,2).std_spectral_base() ;
524 gam.set(3,3) =
exp(2*
ff2) ;
525 gam.set(3,3).std_spectral_base() ;
533 Scalar nphi_ortho(
ww) ;
534 nphi_ortho.mult_rsint() ;
virtual ~ScalarBH()
Destructor.
Scalar ff2
Metric field F_2 of Herdeiro & Radu (2015)
Cmp exp(const Cmp &)
Exponential.
virtual void sauve(FILE *) const
Save in a file.
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.
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).
Base class for coordinate mappings.
Scalar ww
Metric field W of Herdeiro & Radu (2015)
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
virtual void extrinsic_curvature()
Computation of the extrinsic curvature.
virtual void del_deriv() const
Deletes all the derived quantities.
virtual void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Coord tet
coordinate centered on the grid
void operator=(const ScalarBH &)
Assignment to another ScalarBH.
Scalar ff1
Metric field F_1 of Herdeiro & Radu (2015)
Base class for stationary compact objects (under development).
void operator=(const Compobj &)
Assignment to another Compobj.
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Vector beta
Shift vector .
ScalarBH(Map &mp_i, const char *file_name)
Standard constructor.
int get_nzone() const
Returns the number of domains.
Scalar ff0
Metric field F_0 of Herdeiro & Radu (2015)
double rHor
Event horizon coordinate radius.
virtual void del_deriv() const
Deletes all the derived quantities.
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Tbl & set(int l)
Read/write of the value in a given domain.
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Black hole with scalar hair spacetime (under development).
Scalar nn
Lapse function N .
Scalar sfield
Scalar field (modulus of Phi)
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.
Map & mp
Mapping describing the coordinate system (r,theta,phi)
Coord r
r coordinate centered on the grid