LORENE
blackhole_rah_iso.C
1 /*
2  * Methods of class Black_hole to compute the radius of the apparent
3  * horizon in isotropic coordinates
4  *
5  * (see file blackhole.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 blackhole_rah_iso_C[] = "$Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_rah_iso.C,v 1.4 2014/10/13 08:52:46 j_novak Exp $" ;
30 
31 /*
32  * $Id: blackhole_rah_iso.C,v 1.4 2014/10/13 08:52:46 j_novak Exp $
33  * $Log: blackhole_rah_iso.C,v $
34  * Revision 1.4 2014/10/13 08:52:46 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.3 2014/10/06 15:13:02 j_novak
38  * Modified #include directives to use c++ syntax.
39  *
40  * Revision 1.2 2008/05/15 19:30:35 k_taniguchi
41  * Change of some parameters.
42  *
43  * Revision 1.1 2007/06/22 01:20:13 k_taniguchi
44  * *** empty log message ***
45  *
46  *
47  * $Header: /cvsroot/Lorene/C++/Source/Black_hole/blackhole_rah_iso.C,v 1.4 2014/10/13 08:52:46 j_novak Exp $
48  *
49  */
50 
51 // C++ headers
52 //#include <>
53 
54 // C headers
55 #include <cmath>
56 
57 // Lorene headers
58 #include "blackhole.h"
59 #include "utilitaires.h"
60 
61 // Local function
62 namespace Lorene {
63 double ff(double, const double) ;
64 
65  //----------------------------------------------------------//
66  // Radius of the apparent horizon (excised surface) //
67  //----------------------------------------------------------//
68 
69 double Black_hole::rah_iso(bool neumann, bool first) const {
70 
71  // Sets C/M^2 for each case of the lapse boundary condition
72  // --------------------------------------------------------
73  double cc ;
74 
75  if (neumann) { // Neumann boundary condition
76  if (first) { // First condition
77  // d(\alpha \psi)/dr = 0
78  // ---------------------
79  cc = 2. * (sqrt(13.) - 1.) / 3. ;
80  }
81  else { // Second condition
82  // d(\alpha \psi)/dr = (\alpha \psi)/(2 rah)
83  // -----------------------------------------
84  cc = 4. / 3. ;
85  }
86  }
87  else { // Dirichlet boundary condition
88  if (first) { // First condition
89  // (\alpha \psi) = 1/2
90  // -------------------
91  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
92  abort() ;
93  }
94  else { // Second condition
95  // (\alpha \psi) = 1/sqrt(2.) \psi_KS
96  // ----------------------------------
97  cout << "!!!!! WARNING: Not yet prepared !!!!!" << endl ;
98  abort() ;
99  }
100  }
101 
102  int nn = 1000 ;
103  double hh = 1./double(nn) ;
104  double integ = 0. ;
105  double rah ; // rah [M]
106 
107  int mm ;
108  double x1, x2, x3, x4, x5 ;
109 
110  // Boole's Rule (Newton-Cotes Integral) for integration
111  // ----------------------------------------------------
112 
113  assert(nn%4 == 0) ;
114  mm = nn/4 ;
115 
116  for (int i=0; i<mm; i++) {
117 
118  x1 = hh * double(4*i) ;
119  x2 = hh * double(4*i+1) ;
120  x3 = hh * double(4*i+2) ;
121  x4 = hh * double(4*i+3) ;
122  x5 = hh * double(4*i+4) ;
123 
124  integ += (hh/45.) * (14.*ff(x1,cc) + 64.*ff(x2,cc)
125  + 24.*ff(x3,cc) + 64.*ff(x4,cc)
126  + 14.*ff(x5,cc)) ;
127 
128  }
129 
130  rah = 2. * exp(integ) ; // rah : normalized by M
131 
132  return rah ;
133 
134 }
135 
136 //*****************************************************************
137 
138 double ff(double xx, const double cc) {
139 
140  double tcc2 = cc*cc/16. ;
141  double tmp = sqrt(1. - xx + tcc2*pow(xx, 4.)) ;
142 
143  double resu = (-1. + tcc2 * pow(xx, 3.)) / tmp / (1. + tmp) ;
144 
145  return resu ;
146 
147 }
148 }
Cmp exp(const Cmp &)
Exponential.
Definition: cmp_math.C:270
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
Lorene prototypes.
Definition: app_hor.h:64
double rah_iso(bool neumann, bool first) const
Computes a radius of apparent horizon (excised surface) in isotropic coordinates. ...
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348