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