From 81ed6121fa1a6ea7a44bd57513d72947cba8323c Mon Sep 17 00:00:00 2001
From: stla <laurent_step@yahoo.fr>
Date: Fri, 6 Oct 2023 18:50:09 +0200
Subject: [PATCH] fixed conformal torus

---
 R/someMeshes.R                     | 37 +++++++-------
 inst/essais/essai_conformalTorus.R | 77 ++++++++++++++++++++++++++++++
 2 files changed, 95 insertions(+), 19 deletions(-)
 create mode 100644 inst/essais/essai_conformalTorus.R

diff --git a/R/someMeshes.R b/R/someMeshes.R
index b4dbf15..c96b0cb 100644
--- a/R/someMeshes.R
+++ b/R/someMeshes.R
@@ -101,29 +101,28 @@ torusMesh2 <- function(R, r, nu, nv) {
   vs    <- matrix(NA_real_, nrow = 3L, ncol = nunv)
   tris1 <- matrix(NA_integer_, nrow = 3L, ncol = nunv)
   tris2 <- matrix(NA_integer_, nrow = 3L, ncol = nunv)
-  u_ <- seq(0, 2*pi, length.out = nu + 1L)[-1L]
-  cosu_ <- cos(u_)
-  sinu_ <- sin(u_)
-  v_ <- seq(0, 2*pi, length.out = nv + 1L)[-1L]
-  cosv_ <- cos(v_)
-  sinv_ <- sin(v_)
-  kxy <- R*R - r*r
-  kz <- sqrt(kxy) * r
-  kcosu_ <- kxy * cosu_
-  ksinu_ <- kxy * sinu_
-  w_ <- R - r * cosv_
-  h_ <- kz * sinv_ / w_
+  rho <- R/r
+  s <- sqrt(rho*rho - 1)
+  r <- s*r
+  u_ <- seq(0, 2, length.out = nu + 1L)[-1L]
+  ccu_ <- s * cospi(u_)
+  ssu_ <- s * sinpi(u_)
+  v_ <- seq(0, 2, length.out = nv + 1L)[-1L]
+  cospiv_ <- cospi(v_)
+  sinpiv_ <- sinpi(v_)
+  w_ <- rho - cospiv_
+  h_ <- sinpiv_ / w_
   jp1_ <- c(2L:nv, 1L)
   j_ <- 1L:nv
   for(i in 1L:(nu-1L)){
     i_nv <- i*nv
     k1 <- i_nv - nv
     rg <- (k1 + 1L):i_nv
-    kcosu_i <- kcosu_[i]
-    ksinu_i <- ksinu_[i]
-    vs[, rg] <- rbind(
-      kcosu_i / w_,
-      ksinu_i / w_,
+    ccu_i <- ccu_[i]
+    ssu_i <- ssu_[i]
+    vs[, rg] <- r * rbind(
+      ccu_i / w_,
+      ssu_i / w_,
       h_
     )
     k_ <- k1 + j_
@@ -135,8 +134,8 @@ torusMesh2 <- function(R, r, nu, nv) {
   i_nv <- nunv
   k1 <- i_nv - nv
   rg <- (k1 + 1L):i_nv
-  vs[, rg] <- rbind(
-    kxy / w_,
+  vs[, rg] <- r * rbind(
+    s / w_,
     0,
     h_
   )
diff --git a/inst/essais/essai_conformalTorus.R b/inst/essais/essai_conformalTorus.R
new file mode 100644
index 0000000..7062a73
--- /dev/null
+++ b/inst/essais/essai_conformalTorus.R
@@ -0,0 +1,77 @@
+library(rgl)
+library(cgalMeshes)
+
+N <- 1024
+
+g <- function(phi) {
+  r <- sin(phi) / sqrt(1 - sin(phi)^2)
+  R <- cos(phi) / (1 - sin(phi)) - r
+  h <- R/r
+  sqrt(h*h - 1)
+}
+
+phi <- pi/4 # c'est la root de g=1
+r <- sin(phi) / sqrt(1 - sin(phi)^2)
+R <- cos(phi) / (1 - sin(phi)) - r
+
+torus <- function(R, r) {
+  h <- R/r
+  s <- sqrt(h*h - 1)
+  r <- 1/s/r
+  f <- function(u, v) {
+    w <- h - cospi(2*v)
+    rbind(
+      s * cospi(2*u/s) / w,
+      s * sinpi(2*u/s) / w,
+      sinpi(2*v) / w
+    ) / r
+  }
+  parametricMesh(
+    f, c(0, s), c(0, 1), c(TRUE, TRUE),
+    nu = N, nv = N
+  )
+}
+
+mesh <- torus(7, 3)
+summary(t(mesh$vb))
+# r = 1/sqrt((h-1)*(h+1))
+
+
+u_ <- seq(0, 1, length.out = N)
+v_ <- seq(0, 1, length.out = N)
+UV <- as.matrix(
+  expand.grid(V = v_, U = u_)
+)
+
+rot <- function(alpha, uv) {
+  t(rbind(
+    c(cos(alpha), -sin(alpha)),
+    c(sin(alpha),  cos(alpha))
+  ) %*% t(uv))
+}
+
+UVrot <- rot(pi/4, UV)
+
+clrs1 <- ifelse(
+  (floor(4*sqrt(2)*UVrot[, 1L]) %% 2) == (floor(4*sqrt(2)*UVrot[, 2L]) %% 2),
+  "yellow", "navy"
+)
+plot(UV, col = clrs1, asp = 1, pch = ".")
+
+mesh$material <- list(color = clrs1)
+shade3d(mesh, polygon_offset = 1) 
+#bbox3d()
+
+Villarceau <- function(beta, theta0 = 0) {
+  d <- (1 - sin(beta) * sin(phi))
+  cbind(
+    cos(theta0 + beta) * cos(phi) / d,
+    sin(theta0 + beta) * cos(phi) / d,
+    cos(beta) * sin(phi) / d
+  ) 
+}
+
+beta_ <- seq(0, 2*pi, len = 400)
+pts <- Villarceau(beta_)
+
+points3d(pts, size = 7)
\ No newline at end of file