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
29char 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
62namespace Lorene {
63double ff(double, const double) ;
64
65 //----------------------------------------------------------//
66 // Radius of the apparent horizon (excised surface) //
67 //----------------------------------------------------------//
68
69double 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
138double 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}
double rah_iso(bool neumann, bool first) const
Computes a radius of apparent horizon (excised surface) in isotropic coordinates.
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Lorene prototypes.
Definition app_hor.h:64