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 char 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 
66 namespace 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 }
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:220
Lorene prototypes.
Definition: app_hor.h:64
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:387
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:637
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:186
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:108
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