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