-
Notifications
You must be signed in to change notification settings - Fork 1
/
h2tangentvector.cpp
211 lines (165 loc) · 5.18 KB
/
h2tangentvector.cpp
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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
#include "h2tangentvector.h"
#include "h2isometry.h"
H2TangentVector::H2TangentVector()
{
}
H2TangentVector::H2TangentVector(const H2Point &root, const Complex &vector) : root(root), vector(vector)
{
}
H2TangentVector::H2TangentVector(const H2Point &from, const H2Point &to)
{
Complex fromZ, toZ;
fromZ = from.getDiskCoordinate();
toZ = to.getDiskCoordinate();
Complex u = (toZ - fromZ)/(1.0 - conj(fromZ)*toZ);
double distance = H2Point::distance(from,to);
root = from;
vector = distance*u*(1.0-norm(fromZ))/(2.0*std::abs(u));
}
H2TangentVector::H2TangentVector(const H2Point &root) : root(root)
{
vector = Complex(0.0,0.0);
}
H2Point H2TangentVector::exponentiateOld() const
{
double t=length();
Complex z1 = tanh(t/2)*vector/std::abs(vector);
// NB: When translating to the origin via the map (z -z0)/(1 - bar(z0) z), 'vector' is only scaled by a positive factor.
Complex z0 = root.getDiskCoordinate();
Complex z2 = (z1 + z0)/(1.0 + conj(z0)*z1);
assert (norm(z2)<1.0);
return H2Point::fromDiskCoordinate(z2);
}
H2Point H2TangentVector::exponentiate() const
{
return exponentiate(1.0);
}
H2Point H2TangentVector::exponentiate(const double & t) const
{
double L,lambda,absZ,absV;
Complex z,v,outputV;
H2Point output;
z = root.getDiskCoordinate();
v = vector;
absZ = std::abs(z);
absV = std::abs(v);
L = std::abs(t)*length();
Complex u = t >= 0 ? v/absV : -1.0*v/absV;
//lambda = (((1.0-absZ)*exp(L/2.0)) - (1.0+absZ)*exp(-L/2.0))/(((1.0-absZ)*exp(L/2.0)) + (1.0+absZ)*exp(-L/2.0));
lambda = ((1.0-absZ) - ((1.0+absZ)*exp(-1.0*L)))/((1.0-absZ) + ((1.0+absZ)*exp(-1.0*L)));
outputV = ((z*((lambda*absZ) + 1.0))+(u*(lambda+absZ)))/((lambda*absZ+1.0)+(u*conj(z)*(lambda+absZ)));
output.setDiskCoordinate(outputV);
return output;
}
H2Point H2TangentVector::exponentiate(const double &t, const H2TangentVector &v)
{
return v.exponentiate(t);
}
std::vector<H2Point> H2TangentVector::exponentiate(const double &t, const std::vector<H2TangentVector> &V)
{
std::vector<H2Point> out;
out.reserve(V.size());
for (const auto & v : V)
{
out.push_back(v.exponentiate(t));
}
return out;
}
H2TangentVector H2TangentVector::parallelTransportOld(const double &t) const
{
H2Point y0 = root;
H2Point yt = exponentiate(t);
H2Isometry f;
f.setDiskCoordinates(Complex(1.0,0.0), y0.getDiskCoordinate());
H2TangentVector u(f*yt, (1.0 - norm((f*yt).getDiskCoordinate())) * ((f*(*this)).vector));
return (f.inverse())*u;
}
H2TangentVector H2TangentVector::parallelTransport(const double &t) const
{
H2Point yt = exponentiate(t);
Complex v,outV;
double L,absV;
Complex z = root.getDiskCoordinate();
v = vector;
absV = std::abs(v);
// Absolute value t? NO.
L = t*length();
Complex u = v/absV;
Complex thing = cosh(L/2)+u*conj(z)*sinh(L/2);
outV = v/(thing*thing);
H2TangentVector output(yt,outV);
return output;
}
double H2TangentVector::length() const
{
return (2.0/(1.0- norm(root.getDiskCoordinate())))*std::abs(vector);
}
double H2TangentVector::lengthSquared() const
{
double d = (1.0- norm(root.getDiskCoordinate()));
return (4.0/(d*d))*norm(vector);
}
double H2TangentVector::scalProd(const H2TangentVector &v1, const H2TangentVector &v2)
{
assert(v1.root == v2.root);
double d = (1.0- norm(v1.root.getDiskCoordinate()));
return (4.0/(d*d))*real(v1.vector*conj(v2.vector));
}
double H2TangentVector::scalProd(const std::vector<H2TangentVector> &V1, const std::vector<H2TangentVector> &V2)
{
assert(V1.size() == V2.size());
double out = 0.0;
for (uint j=0;j<V1.size();++j)
{
out += scalProd(V1[j],V2[j]);
}
return out;
}
H2Point H2TangentVector::getRoot() const
{
return root;
}
Complex H2TangentVector::getVector() const
{
return vector;
}
std::ostream & operator<<(std::ostream & out, const H2TangentVector & vec)
{
out << "H2TangentVector has root " << vec.getRoot() << " and tangent vector " << vec.getVector();
return out;
}
H2TangentVector operator*(const double &t, const H2TangentVector &v)
{
return H2TangentVector(v.root,t*v.vector);
}
H2TangentVector operator*(const H2Isometry &f, const H2TangentVector &v)
{
Complex u, a;
f.getDiskCoordinates(u, a);
Complex z = v.root.getDiskCoordinate();
Complex df = u*(1 - norm(a))/((1.0 - conj(a)*z)*(1.0 - conj(a)*z));
return H2TangentVector(f*v.root, df*v.vector);
}
H2TangentVector operator+(const H2TangentVector & v1, const H2TangentVector & v2)
{
if (!(v1.root == v2.root))
{
std::cout << "v1.root = " << v1.root << ", v2.root = " << v2.root << std::endl;
}
//assert (v1.root == v2.root);
return H2TangentVector(v1.root,v1.vector+v2.vector);
}
H2TangentVector H2TangentVector::parallelTransport(const double &t, const H2TangentVector &v)
{
return v.parallelTransport(t);
}
std::vector<H2TangentVector> H2TangentVector::parallelTransport(const double &t, const std::vector<H2TangentVector> &V)
{
std::vector<H2TangentVector> out;
out.reserve(V.size());
for (const auto & v : V)
{
out.push_back(v.parallelTransport(t));
}
return out;
}