"""Developable sheet geometry: a generalised cylinder.

The sheet bends along one direction: a unit-speed planar profile curve, given by its turning angle
theta(s), is swept along straight rulings. Paper lengths are preserved exactly. theta is a sum of
components: cylinder, wave, ridge, crease, curl (an edge roll) and fold (an accordion step).
"""

import math

import numpy as np

S_MAX = 200.0
DS = 0.02
THETA_MAX = 1.10


def bent_components(rng, s_view=0.0):
    """Random components of a smoothly bent sheet around s_view."""
    comps = []
    radius = math.exp(rng.uniform(math.log(50.0), math.log(900.0)))
    comps.append(dict(type="cylinder", s0=float(s_view + rng.uniform(-30, 30)),
                      R=float(radius * rng.choice([-1.0, 1.0]))))
    for _ in range(int(rng.integers(0, 2))):
        comps.append(dict(type="wave", s0=float(s_view + rng.uniform(-40, 40)),
                          wavelength=float(rng.uniform(40.0, 150.0)),
                          amp=float(rng.uniform(0.08, 0.45))))
    for _ in range(int(rng.integers(0, 3))):
        comps.append(dict(type="ridge", s0=float(s_view + rng.uniform(-35, 35)),
                          sigma=float(rng.uniform(3.0, 16.0)),
                          amp=float(rng.uniform(0.15, 0.80) * rng.choice([-1.0, 1.0]))))
    for _ in range(int(rng.integers(0, 3))):
        comps.append(dict(type="crease", s0=float(s_view + rng.uniform(-35, 35)),
                          width=float(rng.uniform(0.5, 2.6)), half=float(rng.uniform(0.4, 2.5)),
                          amp=float(rng.uniform(0.15, 0.80) * rng.choice([-1.0, 1.0]))))
    if rng.random() < 0.4:
        comps.append(dict(type="curl", s0=float(s_view + rng.uniform(-60, 60)),
                          R=float(rng.uniform(15.0, 60.0) * rng.choice([-1.0, 1.0])),
                          side=float(rng.choice([-1.0, 1.0])), width=float(rng.uniform(2.0, 8.0))))
    return comps


def eval_theta(s, comps):
    """Turning angle at s. Summation order (smooth components, curls, folds) is part of the
    reproducibility of the published datasets."""
    th = np.zeros_like(s)
    for c in comps:
        d = s - c["s0"]
        if c["type"] == "cylinder":
            th = th + d / c["R"]
        elif c["type"] == "wave":
            th = th + c["amp"] * np.sin(2.0 * math.pi * d / c["wavelength"])
        elif c["type"] == "ridge":
            q = d / c["sigma"]
            th = th - c["amp"] * q * np.exp(-0.5 * q * q)
        elif c["type"] == "crease":
            th = th + 0.5 * c["amp"] * (np.tanh((d + c["half"]) / c["width"]) - np.tanh((d - c["half"]) / c["width"]))
    for c in comps:
        if c["type"] == "curl":
            d = c["side"] * (s - c["s0"]) / c["width"]
            th = th + c["width"] * np.logaddexp(0.0, d) / c["R"]
    for c in comps:
        if c["type"] == "fold":
            th = th + 0.5 * c["amp"] * (1.0 + np.tanh((s - c["s0"]) / c["width"]))
    return th


def _cumtrapz(y, dx):
    out = np.zeros_like(y)
    out[1:] = np.cumsum(0.5 * (y[1:] + y[:-1]) * dx)
    return out


class Profile:
    """Profile curve tabulated on s in [-S_MAX, S_MAX]; theta limited per point to +-THETA_MAX."""

    def __init__(self, comps):
        self.s = np.arange(-S_MAX, S_MAX + DS, DS)
        th = eval_theta(self.s, comps)
        self.theta = THETA_MAX * np.tanh(th / THETA_MAX)
        self.f = _cumtrapz(np.cos(self.theta), DS)
        self.zeta = _cumtrapz(np.sin(self.theta), DS)
        self.f -= np.interp(0.0, self.s, self.f)
        self.zeta -= np.interp(0.0, self.s, self.zeta)

    def at(self, s):
        return np.interp(s, self.s, self.f), np.interp(s, self.s, self.zeta)


def sheet_points(shape, U, V, centre_mm):
    """Paper (U, V) mm -> 3-D points; the sheet is near z = 0, +zeta toward the camera = -z."""
    prof = Profile(shape["comps"])
    phi = math.radians(shape["ruling_deg"])
    N = np.array([math.cos(phi), math.sin(phi), 0.0])
    T = np.array([-math.sin(phi), math.cos(phi), 0.0])
    X, Y = U - centre_mm[0], V - centre_mm[1]
    s = X * N[0] + Y * N[1]
    w = X * T[0] + Y * T[1]
    f, zeta = prof.at(s)
    return f[..., None] * N + w[..., None] * T - zeta[..., None] * np.array([0, 0, 1.0])


def rot(axis, a):
    c, s = math.cos(a), math.sin(a)
    if axis == "x":
        return np.array([[1, 0, 0], [0, c, -s], [0, s, c]], float)
    if axis == "y":
        return np.array([[c, 0, s], [0, 1, 0], [-s, 0, c]], float)
    return np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]], float)