LORENE
bin_bhns_omega_tp.C
1/*
2 * Methods of class Bin_bhns to compute an orbital angular velocity
3 * from two points at the stellar surface
4 *
5 * (see file bin_bhns.h for documentation).
6 *
7 */
8
9/*
10 * Copyright (c) 2006-2007 Keisuke Taniguchi
11 *
12 * This file is part of LORENE.
13 *
14 * LORENE is free software; you can redistribute it and/or modify
15 * it under the terms of the GNU General Public License version 2
16 * as published by the Free Software Foundation.
17 *
18 * LORENE is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU General Public License for more details.
22 *
23 * You should have received a copy of the GNU General Public License
24 * along with LORENE; if not, write to the Free Software
25 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26 *
27 */
28
29char star_bhns_omega_tp_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bin_bhns/bin_bhns_omega_tp.C,v 1.4 2014/10/13 08:52:41 j_novak Exp $" ;
30
31/*
32 * $Id: bin_bhns_omega_tp.C,v 1.4 2014/10/13 08:52:41 j_novak Exp $
33 * $Log: bin_bhns_omega_tp.C,v $
34 * Revision 1.4 2014/10/13 08:52:41 j_novak
35 * Lorene classes and functions now belong to the namespace Lorene.
36 *
37 * Revision 1.3 2014/10/06 15:13:00 j_novak
38 * Modified #include directives to use c++ syntax.
39 *
40 * Revision 1.2 2008/05/15 19:00:27 k_taniguchi
41 * Change of some parameters.
42 *
43 * Revision 1.1 2007/06/22 01:10:00 k_taniguchi
44 * *** empty log message ***
45 *
46 *
47 * $Header: /cvsroot/Lorene/C++/Source/Bin_bhns/bin_bhns_omega_tp.C,v 1.4 2014/10/13 08:52:41 j_novak Exp $
48 *
49 */
50
51// C++ headers
52//#include <>
53
54// C headers
55#include <cmath>
56
57// Lorene headers
58#include "bin_bhns.h"
59#include "unites.h"
60#include "utilitaires.h"
61
62 //---------------------------------------------------//
63 // Orbaital angular velocity from two points //
64 //---------------------------------------------------//
65
66namespace Lorene {
68
69 // Fundamental constants and units
70 // -------------------------------
71 using namespace Unites ;
72
73 if (p_omega_two_points == 0x0) { // a new computation is required
74
75 double omega_two ;
76
77 const Scalar& lapconf = star.get_lapconf_tot() ;
78 const Scalar& confo = star.get_confo_tot() ;
79 const Scalar& psi4 = star.get_psi4() ;
80 const Vector& shift = star.get_shift_tot() ;
81 const Scalar& gam = star.get_gam() ;
82
83 int ii = (star.get_mp()).get_mg()->get_nr(0) - 1 ;
84 int jj = (star.get_mp()).get_mg()->get_nt(0) - 1 ;
85 int ka = 0 ;
86 int kb = (star.get_mp()).get_mg()->get_np(0) / 2 ;
87
88 double psi4_a = psi4.val_grid_point(0,ka,jj,ii) ;
89 double psi4_b = psi4.val_grid_point(0,kb,jj,ii) ;
90 double con2_a = confo.val_grid_point(0,ka,jj,ii)
91 * confo.val_grid_point(0,ka,jj,ii) ;
92 double con2_b = confo.val_grid_point(0,kb,jj,ii)
93 * confo.val_grid_point(0,kb,jj,ii) ;
94 double gam2_a = gam.val_grid_point(0,ka,jj,ii)
95 * gam.val_grid_point(0,ka,jj,ii) ;
96 double gam2_b = gam.val_grid_point(0,kb,jj,ii)
97 * gam.val_grid_point(0,kb,jj,ii) ;
98 double lap2_a = lapconf.val_grid_point(0,ka,jj,ii)
99 * lapconf.val_grid_point(0,ka,jj,ii) ;
100 double lap2_b = lapconf.val_grid_point(0,kb,jj,ii)
101 * lapconf.val_grid_point(0,kb,jj,ii) ;
102 double shiftx_a = shift(1).val_grid_point(0,ka,jj,ii) ;
103 double shiftx_b = shift(1).val_grid_point(0,kb,jj,ii) ;
104 double shifty_a = shift(2).val_grid_point(0,ka,jj,ii) ;
105 double shifty_b = shift(2).val_grid_point(0,kb,jj,ii) ;
106 double shiftz_a = shift(3).val_grid_point(0,ka,jj,ii) ;
107 double shiftz_b = shift(3).val_grid_point(0,kb,jj,ii) ;
108
109 double xns_rot = (star.get_mp()).get_ori_x() - x_rot ;
110 double yns_rot = (star.get_mp()).get_ori_y() - y_rot ;
111
112 double ra = star.ray_eq() ;
113 double rb = star.ray_eq_pi() ;
114
115 if (hole.is_kerrschild()) {
116
117 cout << "!!!!! WARNING: Not yet available !!!!!" << endl ;
118 abort() ;
119 /*
120 double y_separ = (star.get_mp()).get_ori_y() ;
121 double xbh_rot = (hole.get_mp()).get_ori_x() - x_rot ;
122 double mass = ggrav * hole.get_mass_bh() ;
123 double rbh_a = sqrt( (ra+separ)*(ra+separ) + y_separ*y_separ ) ;
124 double rbh_b = sqrt( (-rb+separ)*(-rb+separ) + y_separ*y_separ ) ;
125
126 double msr_a = 2.*mass / pow(rbh_a, 3.) ;
127 double msr_b = 2.*mass / pow(rbh_b, 3.) ;
128
129 double sa = shiftx_a*shiftx_a+shifty_a*shifty_a+shiftz_a*shiftz_a
130 + msr_a * ((ra+separ)*shiftx_a + y_separ*shifty_a)
131 * ((ra+separ)*shiftx_a + y_separ*shifty_a) ;
132
133 double sb = shiftx_b*shiftx_b+shifty_b*shifty_b+shiftz_b*shiftz_b
134 + msr_b * ((-rb+separ)*shiftx_b + y_separ*shifty_b)
135 * ((-rb+separ)*shiftx_b + y_separ*shifty_b) ;
136
137 double ta = -shiftx_a*yns_rot + shifty_a*(ra+xns_rot)
138 + msr_a * ((ra+separ)*shiftx_a + y_separ*shifty_a)
139 * y_separ * xbh_rot ;
140
141 double tb = -shiftx_b*yns_rot + shifty_b*(-rb+xns_rot)
142 + msr_b * ((-rb+separ)*shiftx_b + y_separ*shifty_b)
143 * y_separ * xbh_rot ;
144
145 double ua = yns_rot*yns_rot + (ra+xns_rot)*(ra+xns_rot)
146 + msr_a * y_separ * y_separ * xbh_rot * xbh_rot ;
147
148 double ub = yns_rot*yns_rot + (-rb+xns_rot)*(-rb+xns_rot)
149 + msr_b * y_separ * y_separ * xbh_rot * xbh_rot ;
150
151 // Coefficients : Omega^2 * aaa + 2*Omega * bbb + ccc = 0
152 // ------------------------------------------------------
153
154 double aaa = psi4_a * gam2_a * ua - psi4_b * gam2_b * ub ;
155 double bbb = psi4_a * gam2_a * ta - psi4_b * gam2_b * tb ;
156 double ccc = psi4_a * gam2_a * sa - psi4_b * gam2_b * sb
157 - lap2_a * gam2_a + lap2_b * gam2_b ;
158
159 // Term inside the square root : ddd = bbb*bbb - aaa*ccc
160 // -----------------------------------------------------
161
162 double ddd = bbb*bbb - aaa*ccc ;
163
164 if ( ddd < 0 ) {
165 cout <<
166 "!!! WARNING : Omega (from two points) does not exist !!!"
167 << endl ;
168
169 omega_two = 0. ;
170 }
171 else {
172
173 double omega_1 = (-bbb + sqrt(ddd)) / aaa ;
174 double omega_2 = (-bbb - sqrt(ddd)) / aaa ;
175
176 cout << "Bin_bhns::omega_two_points:" << endl ;
177 cout << " omega_1 : " << omega_1 * f_unit << " [rad/s]"
178 << endl ;
179 cout << " omega_2 : " << omega_2 * f_unit << " [rad/s]"
180 << endl ;
181
182 omega_two = omega_1 ;
183
184 }
185 */
186 }
187 else { // Isotropic coordinates with the maximal slicing
188
189 double sa = shiftx_a*shiftx_a+shifty_a*shifty_a+shiftz_a*shiftz_a ;
190 double sb = shiftx_b*shiftx_b+shifty_b*shifty_b+shiftz_b*shiftz_b ;
191
192 double ta = -shiftx_a*yns_rot + shifty_a*(ra+xns_rot) ;
193 double tb = -shiftx_b*yns_rot + shifty_b*(-rb+xns_rot) ;
194
195 double ua = yns_rot*yns_rot + (ra+xns_rot)*(ra+xns_rot) ;
196 double ub = yns_rot*yns_rot + (-rb+xns_rot)*(-rb+xns_rot) ;
197
198 // Coefficients : Omega^2 * aaa + 2*Omega * bbb + ccc = 0
199 // ------------------------------------------------------
200
201 double aaa = psi4_a * gam2_a * ua - psi4_b * gam2_b * ub ;
202 double bbb = psi4_a * gam2_a * ta - psi4_b * gam2_b * tb ;
203 double ccc = psi4_a * gam2_a * sa - psi4_b * gam2_b * sb
204 - lap2_a * gam2_a / con2_a + lap2_b * gam2_b / con2_b ;
205
206 // Term inside the square root : ddd = bbb*bbb - aaa*ccc
207 // -----------------------------------------------------
208
209 double ddd = bbb*bbb - aaa*ccc ;
210
211 if ( ddd < 0 ) {
212 cout <<
213 "!!! WARNING : Omega (from two points) does not exist !!!"
214 << endl ;
215
216 omega_two = 0. ;
217 }
218 else {
219
220 double omega_1 = (-bbb + sqrt(ddd)) / aaa ;
221 double omega_2 = (-bbb - sqrt(ddd)) / aaa ;
222
223 cout << "Bin_bhns::omega_two_points:" << endl ;
224 cout << " omega_1 : " << omega_1 * f_unit << " [rad/s]"
225 << endl ;
226 cout << " omega_2 : " << omega_2 * f_unit << " [rad/s]"
227 << endl ;
228
229 omega_two = omega_1 ;
230
231 }
232
233 }
234
235 p_omega_two_points = new double( omega_two ) ;
236
237 }
238
239 return *p_omega_two_points ;
240
241}
242}
double y_rot
Absolute Y coordinate of the rotation axis.
Definition bin_bhns.h:89
Hole_bhns hole
Black hole.
Definition bin_bhns.h:72
double * p_omega_two_points
Orbital angular velocity derived from another method.
Definition bin_bhns.h:137
double omega_two_points() const
Orbital angular velocity derived from another method.
double x_rot
Absolute X coordinate of the rotation axis.
Definition bin_bhns.h:86
Star_bhns star
Neutron star.
Definition bin_bhns.h:75
bool is_kerrschild() const
Returns true for a Kerr-Schild background, false for a Conformally flat one.
Definition blackhole.h:218
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
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
const Scalar & get_lapconf_tot() const
Returns the total lapconf function.
Definition star_bhns.h:347
const Scalar & get_gam() const
Returns the Lorentz factor gam.
Definition star_bhns.h:330
const Scalar & get_confo_tot() const
Returns the total conformal factor.
Definition star_bhns.h:391
const Vector & get_shift_tot() const
Returns the part total shift vector.
Definition star_bhns.h:372
const Scalar & get_psi4() const
Returns the fourth power of the total conformal factor.
Definition star_bhns.h:404
const Map & get_mp() const
Returns the mapping.
Definition star.h:355
double ray_eq() const
Coordinate radius at , [r_unit].
double ray_eq_pi() const
Coordinate radius at , [r_unit].
Tensor field of valence 1.
Definition vector.h:188
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.