Back to the parametric Hopf torus
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
In a previous post, 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=(px,py,pz) on the unit sphere
S2 is the circle on
S3 with parametrization:
Cp:(0,2π[⟶S3ϕ⟼Cp(ϕ)
Now consider a spherical curve. That is, let
Γ be a function mapping an
interval I⊂R to
the unit sphere S2. Then the Hopf
cylinder corresponding to
Γ has parameterization
HΓ:I×(0,2π[⟶S3(t,ϕ)⟼CΓ(t)(ϕ).
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
:
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)
A ring cyclide is a Hopf torus. It corresponds to the case when Γ describes a circle on the unit sphere S2. Below is a R function to compute such a circle.
# 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
∈(0,2π[.
Let’s try. We enter three pairs of spherical coordinates and we apply
the circleOnUnitSphere
function:
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:
F <- function(t, phi){ p <- with(circ, center + radius*(cos(t)*i + sin(t)*j)) Stereo(HopfInverse(p, phi)) }
And we plot:
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)
By rotating our spherical circle about the z-axis, we can obtain the linked cyclides. Below is a R function to perform a rotation in spherical coordinates. See this post on my former blog for some explanations.
# 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:
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)
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.