262 lines
9.0 KiB
C++
262 lines
9.0 KiB
C++
bool polar_mod = true, polar_choose = true;
|
|
|
|
namespace reps {
|
|
|
|
template<class N> static void cyclefix(N& a) {
|
|
while(a > + get_deg<N>(180)) a -= get_deg<N>(360);
|
|
while(a < - get_deg<N>(180)) a += get_deg<N>(360);
|
|
}
|
|
|
|
template<class N> static N cyclefix_on(N a) { cyclefix(a); return a; }
|
|
|
|
/** the Taylor polynomial for 1-sqrt(1-y*y) */
|
|
template<class N> N ssqrt(N y) {
|
|
return y*y/2 + y*y*y*y/8 + y*y*y*y*y*y*y/16 + y*y*y*y*y*y*y*y*y/128;
|
|
}
|
|
|
|
TD struct rep_polar2 {
|
|
|
|
using data = D;
|
|
using N = typename D::Number;
|
|
|
|
struct point { N phi, r; };
|
|
struct isometry { N psi, phi, r; }; // spin by psi first
|
|
|
|
static isometry cspin(int i, int j, N alpha) {
|
|
if(i>j) std::swap(i, j), alpha = -alpha;
|
|
return isometry{.psi = -alpha, .phi = N(0), .r = N(0) };
|
|
}
|
|
static isometry cspin90(int i, int j) { return cspin(i, j, get_deg<N>(90)); }
|
|
static isometry lorentz(int i, int j, N alpha) {
|
|
if(i>j) std::swap(i, j);
|
|
if(i == 0) return isometry{.psi = N(0), .phi = N(0), .r = alpha};
|
|
if(i == 1) return isometry{.psi = N(0), .phi = get_deg<N>(90), .r = alpha};
|
|
throw hr::hr_exception("bad lorentz");
|
|
}
|
|
static isometry id() { return isometry{.psi = N(0), .phi = N(0), .r = N(0)}; };
|
|
static point center() { return point{.phi = N(0), .r = N(0)}; };
|
|
|
|
static std::string print(isometry T) {
|
|
return hr::lalign(0, "{phi=", T.phi, " r=", T.r, " psi=", T.psi, "}");
|
|
}
|
|
|
|
static std::string print(point T) {
|
|
return hr::lalign(0, "{phi=", T.phi, " r=", T.r, "}");
|
|
}
|
|
|
|
static isometry apply(isometry T, isometry U, bool need_psi = true) {
|
|
if(T.r == 0) return isometry {.psi = T.psi+U.psi, .phi = T.psi+U.phi, .r = U.r};
|
|
if(U.r == 0) return isometry {.psi = T.psi+U.psi, .phi = T.phi, .r = T.r};
|
|
N alpha = U.phi + T.psi - T.phi;
|
|
if(polar_mod) cyclefix(alpha);
|
|
isometry res;
|
|
N y1 = sinh(U.r) * sin(alpha);
|
|
|
|
auto ca = cos(alpha);
|
|
auto sa = sin(alpha);
|
|
N x1, x2;
|
|
|
|
// choose the appropriate method
|
|
if(polar_choose && ca >= N(0.5)) {
|
|
N u = ca >= N(.999999) ? ssqrt(sa) : N(1) - ca;
|
|
res.r = cosh(T.r + U.r) - u * sinh(T.r) * sinh(U.r);
|
|
x1 = sinh(T.r + U.r) - u * cosh(T.r) * sinh(U.r);
|
|
if(need_psi) x2 = sinh(T.r + U.r) - u * cosh(U.r) * sinh(T.r);
|
|
}
|
|
else if(polar_choose && ca <= N(-0.5)) {
|
|
N u = ca <= N(-.999999) ? ssqrt(-sa) : ca + N(1);
|
|
res.r = cosh(T.r - U.r) + u * sinh(T.r) * sinh(U.r);
|
|
x1 = sinh(T.r - U.r) + u * cosh(T.r) * sinh(U.r);
|
|
if(need_psi) x2 = sinh(U.r - T.r) + u * cosh(U.r) * sinh(T.r);
|
|
}
|
|
else {
|
|
res.r = sinh(T.r) * sinh(U.r) * ca + cosh(T.r) * cosh(U.r);
|
|
x1 = cosh(T.r) * sinh(U.r) * ca + cosh(U.r) * sinh(T.r);
|
|
if(need_psi) x2 = cosh(U.r) * sinh(T.r) * ca + cosh(T.r) * sinh(U.r);
|
|
}
|
|
|
|
if(res.r < N(1)) res.r = N(0); else res.r = acosh(res.r);
|
|
|
|
N beta = (y1 || x1) ? atan2(y1, x1) : N(0);
|
|
res.phi = T.phi + beta;
|
|
if(polar_mod) cyclefix(res.phi);
|
|
|
|
if(need_psi) {
|
|
N y2 = sinh(T.r) * sin(alpha);
|
|
N gamma = (y2 || x2) ? atan2(y2, x2) : N(0);
|
|
res.psi = T.psi + U.psi + beta + gamma - alpha;
|
|
if(polar_mod) cyclefix(res.psi);
|
|
}
|
|
|
|
return res;
|
|
};
|
|
|
|
static point apply(const isometry& T, const point& x) {
|
|
isometry x1 = apply(T, push(x), false);
|
|
return point { .phi = x1.phi, .r = x1.r};
|
|
};
|
|
|
|
static isometry inverse(isometry T) { return isometry{.psi = -T.psi, .phi = cyclefix_on<N>(get_deg<N>(180)+T.phi-T.psi), .r=T.r }; };
|
|
static isometry push(const point& p) { return isometry{.psi = N(0), .phi = p.phi, .r = p.r}; }
|
|
|
|
static N dist0(const point& x) { return x.r; }
|
|
static N angle(const point& x) { return x.phi; }
|
|
static N get_coord(const point& x, int i) {
|
|
if(i == 0) return cos(x.phi) * sinh(x.r);
|
|
if(i == 1) return sin(x.phi) * sinh(x.r);
|
|
if(i == 2) return cosh(x.r);
|
|
throw hr::hr_exception("bad get_coord");
|
|
}
|
|
};
|
|
|
|
TD struct rep_high_polar {
|
|
|
|
using data = D;
|
|
using N = typename D::Number;
|
|
|
|
struct sphere_data {
|
|
using Number = N;
|
|
static constexpr int Dim = D::Dim-1;
|
|
static constexpr int Flipped = -1;
|
|
};
|
|
|
|
using subsphere = rep_linear_nn<sphere_data>;
|
|
|
|
struct point { typename subsphere::point phi; N r; };
|
|
struct isometry { typename subsphere::isometry psi; typename subsphere::point phi; N r; };
|
|
|
|
static isometry cspin(int i, int j, N alpha) {
|
|
return isometry{.psi = subsphere::cspin(i, j, alpha), .phi = subsphere::center(), .r = N(0) };
|
|
}
|
|
static isometry cspin90(int i, int j) {
|
|
return isometry{.psi = subsphere::cspin90(i, j), .phi = subsphere::center(), .r = N(0) };
|
|
}
|
|
static isometry lorentz(int i, int j, N alpha) {
|
|
if(i>j) std::swap(i, j);
|
|
auto is = isometry{.psi = subsphere::id(), .phi = subsphere::center(), .r = alpha};
|
|
is.phi[D::Dim-2] = N(0);
|
|
is.phi[i] = N(1);
|
|
return is;
|
|
}
|
|
static isometry id() { return isometry{.psi = subsphere::id(), .phi = subsphere::center(), .r = N(0)}; }
|
|
static point center() { return point{.phi = subsphere::center(), .r = N(0)}; };
|
|
|
|
static std::string print(isometry T) {
|
|
return hr::lalign(0, "{phi=", subsphere::print(T.phi), " r=", T.r, " psi=", hr::kz(T.psi.values), "}");
|
|
}
|
|
|
|
static std::string print(point T) {
|
|
return hr::lalign(0, "{phi=", subsphere::print(T.phi), " r=", T.r, "}");
|
|
}
|
|
|
|
static isometry apply(isometry T, isometry U, bool need_psi = true) {
|
|
auto apsi = need_psi ? T.psi * U.psi : subsphere::id();
|
|
if(T.r == 0) return isometry {.psi = apsi, .phi = T.psi*U.phi, .r = U.r};
|
|
if(U.r == 0) return isometry {.psi = apsi, .phi = T.phi, .r = T.r};
|
|
auto aphi = T.psi * U.phi;
|
|
auto cos_alpha = inner<sphere_data>(aphi, T.phi);
|
|
auto& ca = cos_alpha;
|
|
|
|
isometry res;
|
|
|
|
N x1, x2;
|
|
|
|
auto orth = (aphi - T.phi * ca);
|
|
N sin_alpha;
|
|
if(ca > N(0.999999) || ca < N(0.999999))
|
|
sin_alpha = pow(sqnorm<sphere_data>(orth), .5);
|
|
else
|
|
sin_alpha = pow(N(1) - ca * ca, .5);
|
|
|
|
if(sin_alpha == N(0)) {
|
|
if(ca >= N(1)) {
|
|
return isometry{.psi = apsi, .phi = T.phi, .r = T.r + U.r };
|
|
}
|
|
if(ca <= N(-1)) {
|
|
if(T.r >= U.r) {
|
|
return isometry{.psi = apsi, .phi = T.phi, .r = T.r - U.r };
|
|
}
|
|
else {
|
|
return isometry{.psi = apsi, .phi = T.phi*-1, .r = U.r - T.r };
|
|
}
|
|
}
|
|
}
|
|
|
|
orth = orth / sin_alpha;
|
|
N y1 = sinh(U.r) * sin_alpha;
|
|
|
|
// choose the appropriate method
|
|
if(polar_choose && ca >= N(0.5)) {
|
|
N u = ca >= N(.999999) ? ssqrt(sin_alpha) : N(1) - ca;
|
|
res.r = cosh(T.r + U.r) - u * sinh(T.r) * sinh(U.r);
|
|
x1 = sinh(T.r + U.r) - u * cosh(T.r) * sinh(U.r);
|
|
if(need_psi) x2 = sinh(T.r + U.r) - u * cosh(U.r) * sinh(T.r);
|
|
}
|
|
else if(polar_choose && ca <= N(-0.5)) {
|
|
N u = ca <= N(-.999999) ? ssqrt(sin_alpha) : ca + N(1); // ca = u - 1
|
|
res.r = cosh(T.r - U.r) + u * sinh(T.r) * sinh(U.r);
|
|
x1 = sinh(T.r - U.r) + u * cosh(T.r) * sinh(U.r);
|
|
if(need_psi) x2 = sinh(U.r - T.r) + u * cosh(U.r) * sinh(T.r);
|
|
}
|
|
else {
|
|
res.r = sinh(T.r) * sinh(U.r) * ca + cosh(T.r) * cosh(U.r);
|
|
x1 = cosh(T.r) * sinh(U.r) * ca + cosh(U.r) * sinh(T.r);
|
|
if(need_psi) x2 = cosh(U.r) * sinh(T.r) * ca + cosh(T.r) * sinh(U.r);
|
|
}
|
|
|
|
if(res.r < N(1)) res.r = N(0); else res.r = acosh(res.r);
|
|
|
|
auto h1 = pow(x1*x1+y1*y1, -0.5);
|
|
N cos_beta = x1*h1, sin_beta = y1*h1;
|
|
|
|
res.phi = T.phi * cos_beta + orth * sin_beta;
|
|
|
|
if(need_psi) {
|
|
N y2 = sinh(T.r) * sin_alpha;
|
|
auto h2 = pow(x2*x2+y2*y2, -0.5);
|
|
N cos_gamma = x2*h2, sin_gamma = y2*h2;
|
|
|
|
// delta = beta + gamma - alpha
|
|
|
|
auto cos_beta_gamma = cos_beta * cos_gamma - sin_beta * sin_gamma;
|
|
auto sin_beta_gamma = cos_beta * sin_gamma + sin_beta * cos_gamma;
|
|
|
|
auto cos_delta = cos_beta_gamma * cos_alpha + sin_beta_gamma * sin_alpha;
|
|
auto sin_delta = sin_beta_gamma * cos_alpha - cos_beta_gamma * sin_alpha;
|
|
|
|
auto phi1 = T.phi * cos_delta + orth * sin_delta;
|
|
auto orth1 = orth * cos_delta - T.phi * sin_delta;
|
|
|
|
auto phi2 = phi1 - T.phi;
|
|
auto orth2 = orth1 - orth;
|
|
typename subsphere::isometry spinner = subsphere::id();
|
|
|
|
// Tv = v + <v, phi> * (phi'-phi) + <v, orth> * (orth'-orth)
|
|
|
|
for(int i=0; i<D::Dim-1; i++)
|
|
for(int j=0; j<D::Dim-1; j++)
|
|
spinner[i][j] += phi2[i] * T.phi[j] + orth2[i] * orth[j];
|
|
|
|
res.psi = spinner * apsi;
|
|
}
|
|
|
|
return res;
|
|
};
|
|
|
|
static point apply(const isometry& T, const point& x) {
|
|
isometry x1 = apply(T, push(x), false);
|
|
return point { .phi = x1.phi, .r = x1.r};
|
|
};
|
|
|
|
static isometry inverse(isometry T) { return isometry{.psi = subsphere::inverse(T.psi), .phi = subsphere::inverse(T.psi)*T.phi*-1, .r=T.r }; };
|
|
static isometry push(const point& p) { return isometry{.psi = subsphere::id(), .phi = p.phi, .r = p.r}; }
|
|
|
|
static N dist0(const point& x) { return x.r; }
|
|
static N angle(const point& x) { return subsphere::angle(x.phi); }
|
|
static N get_coord(const point& x, int i) { if(i == D::Dim-1) return cosh(x.r); else return x.phi[i] * sinh(x.r); }
|
|
};
|
|
|
|
template<class D> using rep_polar = typename std::conditional<D::Dim==3, rep_polar2<D>, rep_high_polar<D>>::type;
|
|
|
|
}
|