--- title: "Back to the parametric Hopf torus" author: "Stéphane Laurent" date: '2020-04-12' tags: R, graphics, rgl, geometry, maths rbloggers: yes output: md_document: variant: markdown preserve_yaml: true html_document: highlight: kate keep_md: no highlighter: pandoc-solarized --- In a [previous post](./HopfTorusParametric.html), I explained how to get a parameterization of a Hopf cylinder or torus. There is a clearer way, which I present here. For the Hopf map, the preimage of a point $p=(p_x,p_y,p_z)$ on the unit sphere $S^2$ is the circle on $S^3$ with parametrization: $$ \begin{array}{ccc} \mathcal{C}_p \colon & (0,2\pi[ & \longrightarrow & S^3 \\ & \phi & \longmapsto & \mathcal{C}_p(\phi) \end{array} $$ where $$ \mathcal{C}_p(\phi) = \frac{1}{\sqrt{2(1+p_z)}} \begin{pmatrix} (1+p_z) \cos(\phi) \\ p_x \sin(\phi) - p_y \cos(\phi) \\ p_x \cos(\phi) + p_y \sin(\phi) \\ (1+p_z) \sin(\phi) \end{pmatrix}. $$ Now consider a spherical curve. That is, let $\Gamma$ be a function mapping an interval $I \subset \mathbb{R}$ to the unit sphere $S^2$. Then the Hopf cylinder corresponding to $\Gamma$ has parameterization $$ \begin{array}{ccc} H_\Gamma \colon & I \times (0,2\pi[ & \longrightarrow & S^3 \\ & (t, \phi) & \longmapsto & \mathcal{C}_{\Gamma(t)}(\phi) \end{array}. $$ Recall the [tennis ball curve](./HopfTorusTennisBall.html) example, given for a real constant $A$ and an integer constant $n$ by: $$ \Gamma(t) = \begin{pmatrix} \sin\bigl(\pi/2 - (\pi/2 - A) \cos(nt)\bigr) \cos\bigl(t + A \sin(2nt)\bigr) \\ \sin\bigl(\pi/2 - (\pi/2 - A) \cos(nt)\bigr) \sin\bigl(t + A \sin(2nt)\bigr) \\ \cos\bigl(\pi/2 - (\pi/2 - A) \cos(nt)\bigr) \end{pmatrix}, \quad t \in (0,2\pi[. $$ ``` {.r .numberLines} A <- 0.44 n <- 3 Gamma <- function(t){ alpha <- pi/2 - (pi/2-A)*cos(n*t) beta <- t + A*sin(2*n*t) c( sin(alpha) * cos(beta), sin(alpha) * sin(beta), cos(alpha) ) } HopfInverse <- function(p, phi){ c( (1+p[3])*cos(phi), p[1]*sin(phi) - p[2]*cos(phi), p[1]*cos(phi) + p[2]*sin(phi), (1+p[3])*sin(phi) ) / sqrt(2*(1+p[3])) } Stereo <- function(q){ 2*q[1:3] / (1-q[4]) } F <- function(t, phi){ Stereo(HopfInverse(Gamma(t), phi)) } ``` Now we're ready to plot the stereographic projection of the Hopf torus with `misc3d`: ``` {.r .numberLines} fx <- Vectorize(function(u,v) F(u,v)[1]) fy <- Vectorize(function(u,v) F(u,v)[2]) fz <- Vectorize(function(u,v) F(u,v)[3]) library(misc3d) parametric3d(fx, fy, fz, umin = 0, umax = 2*pi, vmin = 0, vmax = 2*pi, n = 300, smooth = TRUE, color = "#363940") rgl::view3d(90, 0, zoom = 0.65) ``` ![](./figures/HopfTorus_grey.png) A ring cyclide is a Hopf torus. It corresponds to the case when $\Gamma$ describes a circle on the unit sphere $S^2$. Below is a R function to compute such a circle. ``` {.r .numberLines} # helper function: plane passing by points p1, p2, p3 plane3pts <- function(p1,p2,p3){ xcoef <- (p1[2]-p2[2])*(p2[3]-p3[3])-(p1[3]-p2[3])*(p2[2]-p3[2]) ycoef <- (p1[3]-p2[3])*(p2[1]-p3[1])-(p1[1]-p2[1])*(p2[3]-p3[3]) zcoef <- (p1[1]-p2[1])*(p2[2]-p3[2])-(p1[2]-p2[2])*(p2[1]-p3[1]) offset <- p1[1]*xcoef + p1[2]*ycoef + p1[3]*zcoef c(xcoef, ycoef, zcoef, offset) } # helper function: cross product cross <- function(v, w){ c( v[2] * w[3] - v[3] * w[2], v[3] * w[1] - v[1] * w[3], v[1] * w[2] - v[2] * w[1] ) } # circle passing by points three points p1, p2, p3 # given in Cartesian coordinates circle3pts <- function(p1, p2, p3){ p12 <- (p1+p2)/2 p23 <- (p2+p3)/2 v12 <- p2-p1 v23 <- p3-p2 plane <- plane3pts(p1, p2, p3) A <- rbind(plane[1:3], v12, v23) b <- c(plane[4], sum(p12*v12), sum(p23*v23)) center <- c(solve(A) %*% b) r <- sqrt(c(crossprod(p1-center))) i <- (p1-center) / r normal <- plane[1:3] / sqrt(c(crossprod(plane[1:3]))) list(center = center, radius = r, i = i, j = cross(i,normal)) } # circle parameterization: center + radius*(cos(t)*i + sin(t)*j) # circle on unit sphere passing by three points # given in spherical coordinates circleOnUnitSphere <- function(thph1, thph2, thph3){ theta1 <- thph1[1]; phi1 <- thph1[2] theta2 <- thph2[1]; phi2 <- thph2[2] theta3 <- thph3[1]; phi3 <- thph3[2] p1 <- c(sin(theta1)*cos(phi1), sin(theta1)*sin(phi1), cos(theta1)) p2 <- c(sin(theta2)*cos(phi2), sin(theta2)*sin(phi2), cos(theta2)) p3 <- c(sin(theta3)*cos(phi3), sin(theta3)*sin(phi3), cos(theta3)) circle3pts(p1, p2, p3) } ``` The function returns a list with four elements: a point `center`, a number `radius`, and two vectors `i` and `j`. The parameterization of the spherical circle is then `center + radius*(cos(t)*i + sin(t)*j)` for `t` $\in (0, 2\pi[$. Let's try. We enter three pairs of spherical coordinates and we apply the `circleOnUnitSphere` function: ``` {.r .numberLines} thph1 = c(1.3, 1.5) thph2 = c(1.9, 2.8) thph3 = c(1, 2) circ <- circleOnUnitSphere(thph1, thph2, thph3) ``` Then we define the parametrization of the stereographically projected Hopf torus: ``` {.r .numberLines} F <- function(t, phi){ p <- with(circ, center + radius*(cos(t)*i + sin(t)*j)) Stereo(HopfInverse(p, phi)) } ``` And we plot: ``` {.r .numberLines} parametric3d(fx, fy, fz, umin = 0, umax = 2*pi, vmin = 0, vmax = 2*pi, n = 250, smooth = TRUE, color = "#363940") rgl::view3d(90, 0, zoom = 0.65) ``` ![](./figures/cyclide_grey.png) By rotating our spherical circle about the $z$-axis, we can obtain the [linked cyclides](./linkedCyclidesParametric.html). Below is a R function to perform a rotation in spherical coordinates. See [this post](http://stla.github.io/stlapblog/posts/RotationSphericalCoordinates.html) on my former blog for some explanations. ``` {.r .numberLines} # helper functions: basic rotations #### Rx <- function(alpha){ rbind(c(cos(alpha/2), -1i*sin(alpha/2)), c(-1i*sin(alpha/2), cos(alpha/2))) } Ry <- function(alpha){ rbind(c(cos(alpha/2), -sin(alpha/2)), c(sin(alpha/2), cos(alpha/2))) } Rz <- function(alpha){ rbind(c(exp(-1i*alpha/2), 0), c(0, exp(1i*alpha/2))) } # 3D rotation in spherical coordinates #### #' @description Rotation of a vector given in spherical coordinates. #' @param theta_phi spherical coordinates, a vector containing the #' colatitude (or polar angle), between 0 and pi, and the longitude #' (or azimuthal angle), between 0 and 2pi #' @param axis either a letter 'x', 'y' or 'z', a numeric vector of #' length 2 (the spherical coordinates of the axis), or a numeric #' vector of length 3 (the Cartesian coordinates of the axis) #' @param alpha angle of rotation #' @return The spherical coordinates of the transformed vector. rotation <- function(theta_phi, axis="x", alpha){ if(is.character(axis)){ axis <- match.arg(axis, c("x","y","z")) R <- switch(axis, "x" = Rx(alpha), "y" = Ry(alpha), "z" = Rz(alpha)) }else if(length(axis) == 2){ Theta <- axis[1]; Phi <- axis[2] R <- Rz(Phi) %*% Ry(Theta) %*% Rz(alpha) %*% t(Ry(Theta)) %*% t(Conj(Rz(Phi))) }else if(length(axis) == 3){ axis <- axis / sqrt(c(crossprod(axis))) X <- rbind(c(0,1), c(1,0)) Y <- rbind(c(0,-1i), c(1i,0)) Z <- rbind(c(1,0), c(0,-1)) R <- cos(alpha/2)*diag(2) - 1i*sin(alpha/2) * (axis[1]*X + axis[2]*Y + axis[3]*Z) }else{ stop("`axis` must be either: - a letter ('x', 'y' or 'z') - a numeric vector of length two (spherical coordinates) - a numeric vector of length three (Cartesian coordinates)") } theta <- theta_phi[1]; phi <- theta_phi[2] qubit <- c(cos(theta/2), exp(1i*phi)*sin(theta/2)) newqubit <- R %*% qubit z0 <- newqubit[1,1]; z1 <- newqubit[2,1] c(2*atan(Mod(z1)/Mod(z0)), Arg(z1)-Arg(z0)) } ``` Now, let's rotate our spherical circle and plot: ``` {.r .numberLines} thph1 <- rotation(thph1, "z", 2*pi/3) thph2 <- rotation(thph2, "z", 2*pi/3) thph3 <- rotation(thph3, "z", 2*pi/3) circ <- circleOnUnitSphere(thph1, thph2, thph3) parametric3d(fx, fy, fz, umin = 0, umax = 2*pi, vmin = 0, vmax = 2*pi, n = 250, smooth = TRUE, color = "#363940", add = TRUE) thph1 <- rotation(thph1, "z", 2*pi/3) thph2 <- rotation(thph2, "z", 2*pi/3) thph3 <- rotation(thph3, "z", 2*pi/3) circ <- circleOnUnitSphere(thph1, thph2, thph3) parametric3d(fx, fy, fz, umin = 0, umax = 2*pi, vmin = 0, vmax = 2*pi, n = 250, smooth = TRUE, color = "#363940", add = TRUE) ``` ![](./figures/linkedCyclides_grey.png)