Visualizing Vector Fields and Approximating Curvature on Surfaces

R code to visualize vector fields and approximate curvature

library(ggplot2)
library(dplyr)
library(tidyr)
library(pracma) # for numerical derivatives     

Function to compute metric-induced vector fields on sphere & torus

vector_field_sphere <- function(theta, phi, R) {
  A <- R * sin(phi)
  B <- R
  P <- (pracma::grad(function(x) R * sin(x[2]), c(theta, phi))[2]) / B
  Q <- -(pracma::grad(function(x) R, c(theta, phi))[1]) / A
  return(c(P, Q))
}   

# Function to compute vector field on torus
vector_field_torus <- function(u, v, R, r) {
  A <- R + r * cos(v)
  B <- r
  P <- (pracma::grad(function(x) R + r * cos(x[2]), c(u, v))[2]) / B
  Q <- -(pracma::grad(function(x) r, c(u, v))[1]) / A
  return(c(P, Q))
}

Generate grid points and compute vector fields

R <- 1
r <- 0.3
num_points <- 20
theta_seq <- seq(0, 2 * pi, length.out = num_points)
phi_seq <- seq(0, pi, length.out = num_points)
u_seq <- seq(0, 2 * pi, length.out = num_points)
v_seq <- seq(0, 2 * pi, length.out = num_points)

Create data frames for vector fields

sphere_data <- expand.grid(theta = theta_seq, phi = phi_seq) %>%
  rowwise() %>%
  mutate(vector = list(vector_field_sphere(theta, phi, R))) %>%
  unnest_wider(vector, names_sep = "_") %>%
  rename(P = vector_1, Q = vector_2) |>
  mutate(Q = case_when(is.na(Q) ~ 0, TRUE ~ Q))

torus_data <- expand.grid(u = u_seq, v = v_seq) %>%
  rowwise() %>%
  mutate(vector = list(vector_field_torus(u, v, R, r))) %>%
  unnest_wider(vector, names_sep = "_") %>%
  rename(P = vector_1, Q = vector_2) |>
  mutate(Q = case_when(is.na(Q) ~ 0, TRUE ~ Q))

Select a random point on the torus to illustrate curvature approximation

u_ix <- runif(1,min = 1, max=num_points-1) |> round()
v_ix <- runif(1,min = 1, max=num_points-1) |> round()
# the random point
u <- u_seq[u_ix]
v <- v_seq[v_ix]
# the next points to form a square
u_du <- u_seq[u_ix+1]
v_dv <- v_seq[v_ix+1]
# the length of the sides of the square
du <- u_du - u
dv <- v_dv - v
# compute area of square_data
area_square <- du * dv

# data frame for the square
square_data <- tibble(
  u = c(u, u_du, u_du, u),
  v = c(v, v, v_dv, v_dv)
)

function to compute curl of vector field at (u,v)

curl_fn <- function(.data){
  n = names(.data)
  
# note: !!u is a way to unquote the variable u inside dplyr filter
  filter(.data, .data[[n[1]]] == !!u, .data[[n[2]]] == !!v)  |> glimpse()

  P1 <- .data |> filter(.data[[n[1]]] == !!u,    .data[[n[2]]] == !!v) |> select(P) |> pull()
  Q1 <- .data |> filter(.data[[n[1]]] == !!u,    .data[[n[2]]] == !!v) |> select(Q) |> pull()
  P2 <- .data |> filter(.data[[n[1]]] == !!u,    .data[[n[2]]] == !!v_dv) |> select(P) |> pull()
  Q2 <- .data |> filter(.data[[n[1]]] == !!u_du, .data[[n[2]]] == !!v) |> select(Q) |> pull()
  
  print(paste('- P1 = ', P1, '- Q1 = ',Q1, '- P2 = ',P2, '- Q2 = ',Q2))
  print(Q2)
  curl_value <- ( (Q2 - Q1)/du - (P2 - P1)/dv )
  return(curl_value)
}

– torus case – – – – – – – – – – – – – – – – – – –

curl <- curl_fn(torus_data)

Rows: 1
Columns: 4
$ u <dbl> 1.322776
$ v <dbl> 1.984164
$ P <dbl> -0.9157733
$ Q <dbl> 0
[1] "- P1 =  -0.9157733266467 - Q1 =  0 - P2 =  -0.735723910662144 - Q2 =  0"
[1] 0

# Approximate curvature as curl / area
curvature_approx <- curl / area_square
k_torus <- curl / area_square

print(paste("Approximate curvature at (u,v)=(", 
  round(u,2), ",", round(v,2), "): ", round(curvature_approx, 4)))

[1] "Approximate curvature at (u,v)=( 1.32 , 1.98 ):  -4.9787"

Plot the vector field on the torus with the square highlighted

p_torus <- ggplot(torus_data, aes(x = u, y = v)) +
  geom_rect(xmin = -Inf, xmax = Inf, ymin = pi/2, ymax = 3*pi/2,
            fill='#C2A5CF', color="violet", linetype="dashed", linewidth=1, alpha = 0.1) +
  geom_rect(xmin = -Inf, xmax = Inf, ymin = 0, ymax = pi/2,
            fill='#ACD39E', color="green", linetype="dashed", linewidth=1, alpha = 0.1) +
  geom_rect(xmin = -Inf, xmax = Inf, ymin = 3*pi/2, ymax = 2*pi,
            fill='#ACD39E', color="green", linetype="dashed", linewidth=1, alpha = 0.1) +
  geom_segment(aes(xend = u + P * 0.1, yend = v + Q * 0.1), 
    arrow = arrow(length = unit(0.1, "inches"))) +
  geom_point(data = tibble(u=u, v=v), aes(x=u, y=v), color="blue", size=3) +
  geom_rect(data = square_data, aes(xmin=u, xmax=u_du, ymin=v, ymax=v_dv),
            fill=NA, color="black", linewidth=1) +
  labs(title = "Metric-Induced Vector Field on Torus", x = "u", y = "v",
    caption = paste0('Curl of Loop  = ', round(curl,3),
                     '; Area = ', round(area_square,3),
                     '; Curvature = ', round(curvature_approx,2)))

p_torus

– sphere case – – – – – – – – – – – – – – – – – – –

theta_ix <- runif(1,min = 1, max=num_points-1) |> round()
# pick a point near the transition of signs of curvature
# phi_ix <- round(num_points/2) 
phi_ix <- runif(1,min = 1, max=num_points-1) |> round()
# the random point
u <- theta_seq[theta_ix]
v <- phi_seq[phi_ix]
# the next points to form a square
u_du <- theta_seq[theta_ix+1]
v_dv <- phi_seq[phi_ix+1]
# the length of the sides of the square
du <- u_du - u
dv <- v_dv - v

# data frame for the square
square_data_sphere <- tibble(
  theta = c(u, u_du, u_du, u),
  phi = c(v, v, v_dv, v_dv)
)

# compute area of square_data
area_square_sphere <- du * dv

compute the curl of metric-vector field at (theta,phi)

curl <- curl_fn(sphere_data)

Rows: 1
Columns: 4
$ theta <dbl> 0
$ phi   <dbl> 0
$ P     <dbl> 1
$ Q     <dbl> 0
[1] "- P1 =  0.999999999993889 - Q1 =  0 - P2 =  0.986361303396612 - Q2 =  0"
[1] 0

# Approximate curvature as curl / area
curvature_approx <- curl / area_square_sphere
# print(paste("Approximate curvature at (theta,phi)=(", 
#   round(theta,2), ",", round(phi,2), "): ", round(curvature_approx, 4)))

Plot the vector field on the sphere with the square highlighted

p_sphere <- ggplot(sphere_data, aes(x = theta, y = phi)) +
  geom_rect(xmin = -Inf, xmax = Inf, ymin = pi/2, ymax = pi,
            fill='#ACD39E', color="violet", linetype="dashed", linewidth=1, alpha = 0.1) +
  geom_rect(xmin = -Inf, xmax = Inf, ymin = 0, ymax = pi/2,
            fill='#ACD39E', color="green", linetype="dashed", linewidth=1, alpha = 0.1) +
  geom_segment(aes(xend = theta + P * 0.1, yend = phi + Q * 0.1), 
    arrow = arrow(length = unit(0.1, "inches"))) +
  geom_rect(data = square_data_sphere, aes(xmin=min(theta), xmax=max(theta), ymin=min(phi), ymax=max(phi)),
            fill=NA, color="black", linewidth=1) +
  geom_point(data = tibble(), 
             aes(x=square_data_sphere$theta[1], 
                 y=square_data_sphere$phi[1]), color="blue", size=3) +
  labs(title = "Metric-Induced Vector Field on Sphere", x = "Theta", y = "Phi",
    caption = paste0('Curl = ', round(curl,3),
                     '; Area = ', round(area_square_sphere,3),
                     '; Curvature = ', round(curvature_approx,2))) +
    theme(plot.caption = element_text(hjust = 0))


 p_sphere

Print area of square

#print(paste("Area of square at (theta,phi)=(", 
#  round(theta,2), ",", round(phi,2), "): ", round(area_square_sphere, 4)))

Approximate curvature as curl / area

curvature_approx <- curl / area_square_sphere
#print(paste("Curl at (theta,phi)=(", 
#  round(theta,2), ",", round(phi,2), "): ", round(curl, 4)))

Plot the square used for curl approximation on the torus

# draw a torus using plotly
library(plotly)

# Parameters for the torus
u <- theta_seq
v <- phi_seq

torus_surface_data <- expand.grid(theta = theta_seq, phi = phi_seq) %>% mutate(
    x = (R + r * cos(phi)) * cos(theta),
    y = (R + r * cos(phi)) * sin(theta),
    z = r * sin(phi)
  )

x_torus <- matrix(torus_surface_data$x, nrow = length(theta_seq), ncol = length(phi_seq))
y_torus <- matrix(torus_surface_data$y, nrow = length(theta_seq), ncol = length(phi_seq))
z_torus <- matrix(torus_surface_data$z, nrow = length(theta_seq), ncol = length(phi_seq))

# Generate data for a subset of the torus (e.g., a specific angle range)
theta_subset <- seq(square_data$u[1], square_data$u[2], length.out = 25) # small rectangle on the torus
phi_subset <- seq(square_data$v[1], square_data$v[3], length.out = 25)
# theta_subset <- seq(0, 0.4, length.out = 25) # small rectangle on the torus
# phi_subset <- seq(0, 0.4, length.out = 25)
torus_subset_data <- expand.grid(theta = theta_subset, phi = phi_subset) %>%
  mutate(
    x = (R + r * cos(phi)) * cos(theta),
    y = (R + r * cos(phi)) * sin(theta),
    z = r * sin(phi)
  )

x_subset <- matrix(torus_subset_data$x, nrow = length(theta_subset), ncol = length(phi_subset))
y_subset <- matrix(torus_subset_data$y, nrow = length(theta_subset), ncol = length(phi_subset))
z_subset <- matrix(torus_subset_data$z, nrow = length(theta_subset), ncol = length(phi_subset))


# Define axis limits with a buffer for consistent scaling
# Calculate min/max values manually
xmin <- min(x_torus)
xmax <- max(x_torus)
ymin <- min(y_torus)
ymax <- max(y_torus)
zmin <- min(z_torus)
zmax <- max(z_torus)

# Add a small buffer and ensure plot is centered
buffer <- 0.1
# A more robust approach for equal aspect ratio is to use the max dimension for all axes
max_dim <- max(c(xmax - xmin, ymax - ymin, zmax - zmin)) / 2
center_x <- (xmin + xmax) / 2
center_y <- (ymin + ymax) / 2
center_z <- (zmin + zmax) / 2
axis_range <- c(min(center_x - max_dim - buffer, center_y - max_dim - buffer, center_z - max_dim - buffer),
                max(center_x + max_dim + buffer, center_y + max_dim + buffer, center_z + max_dim + buffer))

# Create the 3D plot
p <- plot_ly(showscale = FALSE) %>%
  add_trace(
    x = x_torus,
    y = y_torus,
    z = z_torus,
    type = "surface",
    name = "Full Torus",
    colorscale = "Viridis",
    opacity = 0.5 # Make it semi-transparent
  )

# Add the subset trace to the existing plot 'p'
p <- p %>%
  add_trace(
    x = x_subset,
    y = y_subset,
    z = z_subset,
    type = "surface",
    name = "Small Rectangle",
    text = paste('Curvature Approx =', round(k_torus, 3)),
    colorscale = "Reds",
    size = 20,
    opacity = 1.0 # Make it opaque
  ) 

p <- p %>%
  layout(
    title = paste0('Torus + square w/ Approximate Curvature = ', round(k_torus, 3)),
    scene = list(
      xaxis = list(title = 'X', range = R*axis_range),
      yaxis = list(title = 'Y', range = axis_range),
      zaxis = list(title = 'Z', range = R*axis_range),
      # set aspect mode to 'cube' for equal spacing across all axes
      aspectmode = 'cube'
    )
  ) 

# Display the plot
p