forked from miguelfranca/BHRaytracer
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Spacetime.hpp
80 lines (60 loc) · 1.86 KB
/
Spacetime.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
#pragma once
#include "Matrix.hpp"
class Spacetime
{
public:
virtual MatrixD metric(const VecD &pos) const = 0;
virtual MatrixD imetric(const VecD &pos) const = 0;
virtual Vec<MatrixD> christoffels(const VecD &pos) const = 0;
virtual double
BH_radius() const = 0; // informs of the radius of an existing BH
// find the 0 component of the velocity to satisfy v_\mu v^\mu = V
// V = -1 for massive particles
// V = 0 for null rays
VecD velocity(VecD pos4, VecD vel3, double V) const;
};
class Flat : public Spacetime
{
public:
Flat() {}
MatrixD metric(const VecD &pos) const override;
MatrixD imetric(const VecD &pos) const override;
Vec<MatrixD> christoffels(const VecD &pos) const override;
double BH_radius() const override
{
return 0.;
} // informs of the radius of an existing BH
};
class Schwarzschild : public Spacetime
{
public:
Schwarzschild(double a_M = 1.);
MatrixD metric(const VecD &pos) const override;
MatrixD imetric(const VecD &pos) const override;
Vec<MatrixD> christoffels(const VecD &pos) const override;
double
BH_radius() const override; // informs of the radius of an existing BH
private:
const double M;
inline double f(double r) const { return 1. - 2. * M / r; }
};
class Kerr : public Spacetime
{
public:
Kerr(double a_M = 1., double a_a = 0.);
MatrixD metric(const VecD &pos) const override;
MatrixD imetric(const VecD &pos) const override;
Vec<MatrixD> christoffels(const VecD &pos) const override;
double
BH_radius() const override; // informs of the radius of an existing BH
private:
const double M, a, a2;
inline double f(double r, double sig) const
{
return 1. - 2. * M * r / sig;
}
inline double sigma(double r2, double costheta2) const
{
return r2 + a2 * costheta2;
}
};