LORENE
bin_bhns_shift_ana.C
1 /*
2  * Method of class Bin_bhns to compute the analytic shift vector
3  *
4  * (see file bin_bhns.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2005 Keisuke Taniguchi
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 
28 
29 
30 /*
31  * $Id: bin_bhns_shift_ana.C,v 1.4 2016/12/05 16:17:45 j_novak Exp $
32  * $Log: bin_bhns_shift_ana.C,v $
33  * Revision 1.4 2016/12/05 16:17:45 j_novak
34  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
35  *
36  * Revision 1.3 2014/10/13 08:52:41 j_novak
37  * Lorene classes and functions now belong to the namespace Lorene.
38  *
39  * Revision 1.2 2014/10/06 15:13:00 j_novak
40  * Modified #include directives to use c++ syntax.
41  *
42  * Revision 1.1 2007/06/22 01:11:08 k_taniguchi
43  * *** empty log message ***
44  *
45  *
46  * $Header: /cvsroot/Lorene/C++/Source/Bin_bhns/bin_bhns_shift_ana.C,v 1.4 2016/12/05 16:17:45 j_novak Exp $
47  *
48  */
49 
50 // C++ headers
51 //#include <>
52 
53 // C headers
54 #include <cmath>
55 
56 // Lorene headers
57 #include "bin_bhns.h"
58 #include "utilitaires.h"
59 #include "unites.h"
60 
61 namespace Lorene {
62 void Bin_bhns::shift_analytic(double reduce_shift_bh, double reduce_shift_ns)
63 {
64 
65  using namespace Unites ;
66 
67  double massbh = hole.get_mass_bh() ;
68  double massns = star.mass_g_bhns() ;
69  double mass_bh = ggrav * massbh ;
70  double mass_ns = ggrav * massns ;
71 
72  double mass_tot = mass_bh + mass_ns ;
73 
74  double comb = mass_bh * mass_ns * omega * separ / mass_tot ;
75 
76  //-----------------------------------------//
77  // Shift vector for the black hole //
78  //-----------------------------------------//
79 
80  const Vector& shift_bh_old = hole.get_shift_auto_rs() ;
81 
82  const Map& mp_bh = hole.get_mp() ;
83  Scalar xx_bh(mp_bh) ;
84  xx_bh = mp_bh.x ;
85  xx_bh.std_spectral_base() ;
86  Scalar yy_bh(mp_bh) ;
87  yy_bh = mp_bh.y ;
88  yy_bh.std_spectral_base() ;
89  Scalar zz_bh(mp_bh) ;
90  zz_bh = mp_bh.z ;
91  zz_bh.std_spectral_base() ;
92  Scalar rr_bh(mp_bh) ;
93  rr_bh = mp_bh.r ;
94  rr_bh.std_spectral_base() ;
95  Scalar st_bh(mp_bh) ;
96  st_bh = mp_bh.sint ;
97  st_bh.std_spectral_base() ;
98  Scalar ct_bh(mp_bh) ;
99  ct_bh = mp_bh.cost ;
100  ct_bh.std_spectral_base() ;
101  Scalar sp_bh(mp_bh) ;
102  sp_bh = mp_bh.sinp ;
103  sp_bh.std_spectral_base() ;
104  Scalar cp_bh(mp_bh) ;
105  cp_bh = mp_bh.cosp ;
106  cp_bh.std_spectral_base() ;
107 
108  double rad_bh = rr_bh.val_grid_point(1, 0, 0, 0) ;
109 
110  Scalar x_bh_ex(mp_bh) ;
111  Scalar y_bh_ex(mp_bh) ;
112  Scalar z_bh_ex(mp_bh) ;
113 
114  if (hole.is_irrotational()) {
115 
116  // x component
117  // -----------
118  x_bh_ex = 0.2 * comb * rad_bh * rad_bh
119  * st_bh * st_bh * cp_bh * sp_bh / pow(rr_bh, 3.) ;
120  x_bh_ex.annule_domain(0) ;
121  x_bh_ex.std_spectral_base() ;
122 
123  (hole.set_shift_auto_rs()).set(1) = shift_bh_old(1)
124  + reduce_shift_bh * x_bh_ex ;
125 
126  // y component
127  // -----------
128  y_bh_ex = 0.5 * comb * (7. + 0.2*rad_bh*rad_bh/rr_bh/rr_bh) / rr_bh
129  + 0.5 * comb * pow(st_bh*sp_bh,2.)
130  * (1. - 0.6*rad_bh*rad_bh/rr_bh/rr_bh) / rr_bh ;
131  y_bh_ex.annule_domain(0) ;
132  y_bh_ex.std_spectral_base() ;
133 
134  (hole.set_shift_auto_rs()).set(2) = shift_bh_old(2)
135  + reduce_shift_bh * y_bh_ex ;
136 
137  // z component
138  // -----------
139  z_bh_ex = 0.5 * comb * st_bh * sp_bh * ct_bh
140  * (1.-0.6*rad_bh*rad_bh/rr_bh/rr_bh) / rr_bh ;
141  z_bh_ex.annule_domain(0) ;
142  z_bh_ex.std_spectral_base() ;
143 
144  (hole.set_shift_auto_rs()).set(3) = shift_bh_old(3)
145  + reduce_shift_bh * z_bh_ex ;
146 
147  (hole.set_shift_auto_rs()).std_spectral_base() ;
148 
149  }
150  else { // Corotational
151 
152  // x component
153  // -----------
154  x_bh_ex = - 0.6 * mass_ns * omega * rad_bh * rad_bh
155  * st_bh * sp_bh / pow(rr_bh, 2.)
156  + 0.5 * comb * st_bh * st_bh * cp_bh * sp_bh
157  * (1. - 0.6*rad_bh*rad_bh/rr_bh/rr_bh) / rr_bh
158  - 0.6*mass_bh*omega*rad_bh*rad_bh*pow(st_bh,3.)*pow(cp_bh,2.)*sp_bh
159  /pow(rr_bh, 2.) ;
160  x_bh_ex.annule_domain(0) ;
161  x_bh_ex.std_spectral_base() ;
162 
163  (hole.set_shift_auto_rs()).set(1) = shift_bh_old(1)
164  + reduce_shift_bh * x_bh_ex ;
165 
166  // y component
167  // -----------
168  y_bh_ex = 0.5 * comb * (7. + 0.2*rad_bh*rad_bh/rr_bh/rr_bh) / rr_bh
169  - 0.6 * mass_bh * omega * rad_bh * rad_bh * st_bh * cp_bh
170  / pow(rr_bh, 2.)
171  + 0.5 * comb * pow(st_bh*sp_bh,2.)
172  * (1. - 0.6*rad_bh*rad_bh/rr_bh/rr_bh) / rr_bh
173  - 0.6*mass_bh*omega*rad_bh*rad_bh*pow(st_bh,3.)*cp_bh*pow(sp_bh,2.)
174  /pow(rr_bh, 2.) ;
175  y_bh_ex.annule_domain(0) ;
176  y_bh_ex.std_spectral_base() ;
177 
178  (hole.set_shift_auto_rs()).set(2) = shift_bh_old(2)
179  + reduce_shift_bh * y_bh_ex ;
180 
181  // z component
182  // -----------
183  z_bh_ex = 0.5 * comb * st_bh * cp_bh * ct_bh
184  * (1. - 0.6*rad_bh*rad_bh/rr_bh/rr_bh) / rr_bh
185  - 0.6*mass_bh*omega*rad_bh*rad_bh*st_bh*st_bh*cp_bh*sp_bh*ct_bh
186  / pow(rr_bh, 2.) ;
187  z_bh_ex.annule_domain(0) ;
188  z_bh_ex.std_spectral_base() ;
189 
190  (hole.set_shift_auto_rs()).set(3) = shift_bh_old(3)
191  + reduce_shift_bh * z_bh_ex ;
192 
193  (hole.set_shift_auto_rs()).std_spectral_base() ;
194 
195  }
196 
197 
198  //-------------------------------------------//
199  // Shift vector for the neutron star //
200  //-------------------------------------------//
201  int nzet = star.get_nzet() ;
202  int nz_ns = (star.get_mp()).get_mg()->get_nzone() ;
203 
204  const Map& mp_ns = star.get_mp() ;
205  Scalar xx_ns(mp_ns) ;
206  xx_ns = mp_ns.x ;
207  xx_ns.std_spectral_base() ;
208  Scalar yy_ns(mp_ns) ;
209  yy_ns = mp_ns.y ;
210  yy_ns.std_spectral_base() ;
211  Scalar zz_ns(mp_ns) ;
212  zz_ns = mp_ns.z ;
213  zz_ns.std_spectral_base() ;
214  Scalar rr_ns(mp_ns) ;
215  rr_ns = mp_ns.r ;
216  rr_ns.std_spectral_base() ;
217  Scalar st_ns(mp_ns) ;
218  st_ns = mp_ns.sint ;
219  st_ns.std_spectral_base() ;
220  Scalar ct_ns(mp_ns) ;
221  ct_ns = mp_ns.cost ;
222  ct_ns.std_spectral_base() ;
223  Scalar sp_ns(mp_ns) ;
224  sp_ns = mp_ns.sinp ;
225  sp_ns.std_spectral_base() ;
226  Scalar cp_ns(mp_ns) ;
227  cp_ns = mp_ns.cosp ;
228  cp_ns.std_spectral_base() ;
229 
230  double rad_ns = rr_ns.val_grid_point(1, 0, 0, 0) ;
231 
232  Scalar x_ns_in(mp_ns) ;
233  Scalar x_ns_ex(mp_ns) ;
234  Scalar y_ns_in(mp_ns) ;
235  Scalar y_ns_ex(mp_ns) ;
236  Scalar z_ns_in(mp_ns) ;
237  Scalar z_ns_ex(mp_ns) ;
238 
239  if (star.is_irrotational()) {
240 
241  // x component
242  // -----------
243  x_ns_in = - 0.2 * comb * xx_ns * yy_ns / pow(rad_ns, 3.) ;
244  x_ns_in.annule(nzet, nz_ns-1) ;
245  x_ns_in.std_spectral_base() ;
246 
247  x_ns_ex = - 0.2 * comb * rad_ns * rad_ns
248  * st_ns * st_ns * cp_ns * sp_ns / pow(rr_ns, 3.) ;
249  x_ns_ex.annule(0, nzet-1) ;
250  x_ns_ex.std_spectral_base() ;
251 
252  (star.set_shift_auto()).set(1) = reduce_shift_ns
253  * (x_ns_in + x_ns_ex) ;
254 
255  // y component
256  // -----------
257  y_ns_in = - 0.5 * comb * (11. - 3.8*rr_ns*rr_ns/rad_ns/rad_ns) / rad_ns
258  - 0.2 * comb * yy_ns * yy_ns / pow(rad_ns, 3.) ;
259  y_ns_in.annule(nzet, nz_ns-1) ;
260  y_ns_in.std_spectral_base() ;
261 
262  y_ns_ex = - 0.5 * comb * (7. + 0.2*rad_ns*rad_ns/rr_ns/rr_ns) / rr_ns
263  - 0.5 * comb * pow(st_ns*sp_ns,2.)
264  * (1. - 0.6*rad_ns*rad_ns/rr_ns/rr_ns) / rr_ns ;
265  y_ns_ex.annule(0, nzet-1) ;
266  y_ns_ex.std_spectral_base() ;
267 
268  (star.set_shift_auto()).set(2) = reduce_shift_ns
269  * (y_ns_in + y_ns_ex) ;
270 
271  // z component
272  // -----------
273  z_ns_in = - 0.2 * comb * yy_ns * zz_ns / pow(rad_ns, 3.) ;
274  z_ns_in.annule(nzet, nz_ns-1) ;
275  z_ns_in.std_spectral_base() ;
276 
277  z_ns_ex = - 0.5 * comb * st_ns * sp_ns * ct_ns
278  * (1.-0.6*rad_ns*rad_ns/rr_ns/rr_ns) / rr_ns ;
279  z_ns_ex.annule(0, nzet-1) ;
280  z_ns_ex.std_spectral_base() ;
281 
282  (star.set_shift_auto()).set(3) = reduce_shift_ns
283  * (z_ns_in + z_ns_ex) ;
284 
285  }
286  else { // Corotational
287 
288  // x component
289  // -----------
290  x_ns_in = 1.5 * mass_ns * omega * yy_ns
291  * (1. - 0.6*rr_ns*rr_ns/rad_ns/rad_ns) / rad_ns
292  - 0.2 * comb * xx_ns * yy_ns / pow(rad_ns, 3.)
293  + 0.6 * mass_ns * omega * xx_ns * xx_ns * yy_ns / pow(rad_ns, 3.) ;
294  x_ns_in.annule(nzet, nz_ns-1) ;
295  x_ns_in.std_spectral_base() ;
296 
297  x_ns_ex = 0.6 * mass_ns * omega * rad_ns * rad_ns
298  * st_ns * sp_ns / pow(rr_ns, 2.)
299  - 0.5 * comb * st_ns * st_ns * cp_ns * sp_ns
300  * (1. - 0.6*rad_ns*rad_ns/rr_ns/rr_ns) / rr_ns
301  + 0.6*mass_ns*omega*rad_ns*rad_ns*pow(st_ns,3.)*pow(cp_ns,2.)*sp_ns
302  /pow(rr_ns, 2.) ;
303  x_ns_ex.annule(0, nzet-1) ;
304  x_ns_ex.std_spectral_base() ;
305 
306  (star.set_shift_auto()).set(1) = reduce_shift_ns
307  * (x_ns_in + x_ns_ex) ;
308 
309  // y component
310  // -----------
311  y_ns_in = - 0.5 * comb * (11. - 3.8*rr_ns*rr_ns/rad_ns/rad_ns) / rad_ns
312  + 1.5 * mass_ns * omega * xx_ns
313  * (1. - 0.6*rr_ns*rr_ns/rad_ns/rad_ns) / rad_ns
314  - 0.2 * comb * yy_ns * yy_ns / pow(rad_ns, 3.)
315  + 0.6 * mass_ns * omega * xx_ns * yy_ns * yy_ns / pow(rad_ns, 3.) ;
316  y_ns_in.annule(nzet, nz_ns-1) ;
317  y_ns_in.std_spectral_base() ;
318 
319  y_ns_ex = - 0.5 * comb * (7. + 0.2*rad_ns*rad_ns/rr_ns/rr_ns) / rr_ns
320  + 0.6 * mass_ns * omega * rad_ns * rad_ns * st_ns * cp_ns
321  / pow(rr_ns, 2.)
322  - 0.5 * comb * pow(st_ns*sp_ns,2.)
323  * (1. - 0.6*rad_ns*rad_ns/rr_ns/rr_ns) / rr_ns
324  + 0.6*mass_ns*omega*rad_ns*rad_ns*pow(st_ns,3.)*cp_ns*pow(sp_ns,2.)
325  /pow(rr_ns, 2.) ;
326  y_ns_ex.annule(0, nzet-1) ;
327  y_ns_ex.std_spectral_base() ;
328 
329  (star.set_shift_auto()).set(2) = reduce_shift_ns
330  * (y_ns_in + y_ns_ex) ;
331 
332  // z component
333  // -----------
334  z_ns_in = - 0.2 * comb * yy_ns * zz_ns / pow(rad_ns, 3.)
335  + 0.6 * mass_ns * omega * xx_ns * yy_ns * zz_ns / pow(rad_ns, 3.) ;
336  z_ns_in.annule(nzet, nz_ns-1) ;
337  z_ns_in.std_spectral_base() ;
338 
339  z_ns_ex = - 0.5 * comb * st_ns * cp_ns * ct_ns
340  * (1. - 0.6*rad_ns*rad_ns/rr_ns/rr_ns) / rr_ns
341  + 0.6*mass_ns*omega*rad_ns*rad_ns*st_ns*st_ns*cp_ns*sp_ns*ct_ns
342  / pow(rr_ns, 2.) ;
343  z_ns_ex.annule(0, nzet-1) ;
344  z_ns_ex.std_spectral_base() ;
345 
346  (star.set_shift_auto()).set(3) = reduce_shift_ns
347  * (z_ns_in + z_ns_ex) ;
348 
349  }
350 
351  (star.set_shift_auto()).std_spectral_base() ;
352 
353 }
354 }
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:675
double get_mass_bh() const
Returns the gravitational mass of BH [{ m_unit}].
Definition: blackhole.h:221
Hole_bhns hole
Black hole.
Definition: bin_bhns.h:72
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:397
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Base class for coordinate mappings.
Definition: map.h:688
Vector & set_shift_auto_rs()
Read/write of the shift vector generated by the black hole.
Definition: hole_bhns.C:529
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
const Map & get_mp() const
Returns the mapping.
Definition: star.h:355
void shift_analytic(double reduce_shift_bh, double reduce_shift_ns)
Sets some analytical template for the initial shift vector.
Tensor field of valence 1.
Definition: vector.h:188
double separ
Absolute orbital separation between two centers of BH and NS.
Definition: bin_bhns.h:83
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:643
Coord sint
Definition: map.h:739
double omega
Angular velocity with respect to an asymptotically inertial observer.
Definition: bin_bhns.h:80
const Vector & get_shift_auto_rs() const
Returns the part of the shift vector from numerical computation.
Definition: hole_bhns.h:405
Star_bhns star
Neutron star.
Definition: bin_bhns.h:75
bool is_irrotational() const
Returns true for an irrotational black hole, false for a corotating one.
Definition: hole_bhns.h:359
const Map & get_mp() const
Returns the mapping.
Definition: blackhole.h:213
Coord sinp
Definition: map.h:741
Vector & set_shift_auto()
Read/write of the shift vector generated by the neutron star.
Definition: star_bhns.C:438
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
int get_nzet() const
Returns the number of domains occupied by the star.
Definition: star.h:358
Coord y
y coordinate centered on the grid
Definition: map.h:745
Coord cosp
Definition: map.h:742
Coord x
x coordinate centered on the grid
Definition: map.h:744
bool is_irrotational() const
Returns true for an irrotational motion, false for a corotating one.
Definition: star_bhns.h:302
Coord z
z coordinate centered on the grid
Definition: map.h:746
Coord r
r coordinate centered on the grid
Definition: map.h:736
Coord cost
Definition: map.h:740