{"spec_id":"contour-3d","library":"ggplot2","language":"r","code":"#' anyplot.ai\n#' contour-3d: 3D Contour Plot\n#' Library: ggplot2 3.5.1 | R 4.4.1\n#' Quality: 93/100 | Created: 2026-09-10\n\nlibrary(ggplot2)\nlibrary(dplyr)\nlibrary(ragg)\n\nset.seed(42)\n\n# --- Theme tokens -----------------------------------------------------------\nTHEME    <- Sys.getenv(\"ANYPLOT_THEME\", \"light\")\nPAGE_BG  <- if (THEME == \"light\") \"#FAF8F1\" else \"#1A1A17\"\nINK      <- if (THEME == \"light\") \"#1A1A17\" else \"#F0EFE8\"\nINK_SOFT <- if (THEME == \"light\") \"#4A4A44\" else \"#B8B7B0\"\n\n# --- Camera: orthographic projection (elevation 28, azimuth 40) -------------\n# ggplot2 has no 3D grammar, so the surface mesh and its contour isolines are\n# projected to 2D screen coordinates ourselves (the same technique any static\n# 3D renderer uses under the hood), then drawn with geom_polygon/geom_path.\nelev <- 28 * pi / 180\nazim <- 40 * pi / 180\n\nview_dir <- c(cos(elev) * cos(azim), cos(elev) * sin(azim), sin(elev))\nworld_up <- c(0, 0, 1)\nright_axis <- c(\n  view_dir[2] * world_up[3] - view_dir[3] * world_up[2],\n  view_dir[3] * world_up[1] - view_dir[1] * world_up[3],\n  view_dir[1] * world_up[2] - view_dir[2] * world_up[1]\n)\nright_axis <- right_axis / sqrt(sum(right_axis^2))\nup_axis <- c(\n  right_axis[2] * view_dir[3] - right_axis[3] * view_dir[2],\n  right_axis[3] * view_dir[1] - right_axis[1] * view_dir[3],\n  right_axis[1] * view_dir[2] - right_axis[2] * view_dir[1]\n)\n\nZ_LIFT <- 1.3  # visual height exaggeration so the landscape reads clearly\nproject_x <- function(x, y, z) x * right_axis[1] + y * right_axis[2] + z * Z_LIFT * right_axis[3]\nproject_y <- function(x, y, z) x * up_axis[1]    + y * up_axis[2]    + z * Z_LIFT * up_axis[3]\n\n# --- Data: electric field intensity around a point charge, 40x40 grid -------\n# A single smooth peak keeps the equipotential isolines nested and legible in\n# projection (a multi-peak landscape produces overlapping loops per level).\ngrid_n <- 40\nx_vals <- seq(-4, 4, length.out = grid_n)\ny_vals <- seq(-4, 4, length.out = grid_n)\n\nfield_intensity <- function(x1, x2) {\n  3.4 * exp(-((x1 - 0.4)^2 / 5.0 + (x2 + 0.3)^2 / 7.2))\n}\nz_mat <- outer(x_vals, y_vals, field_intensity)\n\nz_min <- min(z_mat)\nz_max <- max(z_mat)\nz_span <- z_max - z_min\nfloor_z <- z_min - 0.35 * z_span\nceil_z  <- z_max + 0.15 * z_span\n\nx_min <- min(x_vals); x_max <- max(x_vals)\ny_min <- min(y_vals); y_max <- max(y_vals)\n\n# --- Bilinear interpolation so contour isolines get a surface height --------\ninterp_z <- function(x, y) {\n  ix <- findInterval(x, x_vals, all.inside = TRUE)\n  iy <- findInterval(y, y_vals, all.inside = TRUE)\n  x0 <- x_vals[ix]; x1 <- x_vals[ix + 1]\n  y0 <- y_vals[iy]; y1 <- y_vals[iy + 1]\n  tx <- (x - x0) / (x1 - x0)\n  ty <- (y - y0) / (y1 - y0)\n  z00 <- z_mat[cbind(ix, iy)]\n  z10 <- z_mat[cbind(ix + 1, iy)]\n  z01 <- z_mat[cbind(ix, iy + 1)]\n  z11 <- z_mat[cbind(ix + 1, iy + 1)]\n  (z00 * (1 - tx) + z10 * tx) * (1 - ty) + (z01 * (1 - tx) + z11 * tx) * ty\n}\n\n# --- Surface as filled quads, painter's algorithm (far cells drawn first) ---\nsurface_quads <- vector(\"list\", (grid_n - 1) * (grid_n - 1))\nslot <- 0\nfor (i in seq_len(grid_n - 1)) {\n  for (j in seq_len(grid_n - 1)) {\n    slot <- slot + 1\n    xs <- x_vals[c(i, i + 1, i + 1, i)]\n    ys <- y_vals[c(j, j, j + 1, j + 1)]\n    zs <- z_mat[cbind(c(i, i + 1, i + 1, i), c(j, j, j + 1, j + 1))]\n    surface_quads[[slot]] <- data.frame(\n      poly_id = slot,\n      depth   = i + j,\n      px      = project_x(xs, ys, zs),\n      py      = project_y(xs, ys, zs),\n      z_mid   = mean(zs)\n    )\n  }\n}\nsurface_df <- bind_rows(surface_quads) |> arrange(desc(depth), poly_id)\n\n# --- Contour isolines (equipotential lines): on the surface, and projected\n# onto the base plane for reference ------------------------------------------\ncontour_levels <- pretty(c(z_min, z_max), n = 6)\ncontour_levels <- contour_levels[contour_levels > z_min & contour_levels < z_max]\nraw_lines <- grDevices::contourLines(x_vals, y_vals, z_mat, levels = contour_levels)\n\nsurface_lines <- bind_rows(lapply(seq_along(raw_lines), function(k) {\n  ln <- raw_lines[[k]]\n  z_line <- interp_z(ln$x, ln$y)\n  data.frame(line_id = k, px = project_x(ln$x, ln$y, z_line), py = project_y(ln$x, ln$y, z_line))\n}))\nbase_lines <- bind_rows(lapply(seq_along(raw_lines), function(k) {\n  ln <- raw_lines[[k]]\n  data.frame(line_id = k, px = project_x(ln$x, ln$y, floor_z), py = project_y(ln$x, ln$y, floor_z))\n}))\n\n# --- Axis box: three edges meeting at the front-left-bottom corner ----------\naxis_lines <- data.frame(\n  x    = c(x_min, x_min, x_min),\n  y    = c(y_min, y_min, y_min),\n  z    = c(floor_z, floor_z, floor_z),\n  xend = c(x_max, x_min, x_min),\n  yend = c(y_min, y_max, y_min),\n  zend = c(floor_z, floor_z, ceil_z)\n)\naxis_lines$px    <- project_x(axis_lines$x, axis_lines$y, axis_lines$z)\naxis_lines$py    <- project_y(axis_lines$x, axis_lines$y, axis_lines$z)\naxis_lines$pxend <- project_x(axis_lines$xend, axis_lines$yend, axis_lines$zend)\naxis_lines$pyend <- project_y(axis_lines$xend, axis_lines$yend, axis_lines$zend)\n\nx_breaks <- pretty(x_vals, n = 4); x_breaks <- x_breaks[x_breaks >= x_min & x_breaks <= x_max]\ny_breaks <- pretty(y_vals, n = 4); y_breaks <- y_breaks[y_breaks >= y_min & y_breaks <= y_max]\nz_breaks <- pretty(c(floor_z, ceil_z), n = 4); z_breaks <- z_breaks[z_breaks >= floor_z & z_breaks <= ceil_z]\n\n# Large offsets keep the x-tick and y-tick label clusters from colliding near\n# the shared corner — both grow away from the box along their own axis line.\n# 0.45x the axis range clears the surface's projected silhouette even for the\n# negative-side ticks nearest the tall part of the mesh (0.3x left \"-4\"/\"-2\"\n# sitting on top of the surface after the Z_LIFT projection); the matching\n# Z_TICK_OFFSET below is widened in step so the two label columns don't swap\n# one collision (with the mesh) for another (with each other).\ntick_gap <- 0.45 * (x_max - x_min)\nticks <- rbind(\n  data.frame(x = x_breaks, y = y_min - tick_gap, z = floor_z, label = x_breaks),\n  data.frame(x = x_min - tick_gap, y = y_breaks, z = floor_z, label = y_breaks)\n)\nticks$px <- project_x(ticks$x, ticks$y, ticks$z)\nticks$py <- project_y(ticks$x, ticks$y, ticks$z)\n\n# Z ticks sit on the vertical axis line itself; nudge the label text\n# (not the axis line) sideways into the open gap left of the mesh. Widened\n# alongside tick_gap above so the z-tick column stays clear of both the mesh\n# and the (now farther-out) y-tick label column.\nZ_TICK_OFFSET <- 4.5\nz_ticks <- data.frame(x = x_min, y = y_min, z = z_breaks, label = z_breaks)\nz_ticks$px <- project_x(z_ticks$x, z_ticks$y, z_ticks$z) - Z_TICK_OFFSET\nz_ticks$py <- project_y(z_ticks$x, z_ticks$y, z_ticks$z)\n\naxis_labels <- data.frame(\n  x     = c(x_max + 0.7, x_min, x_min),\n  y     = c(y_min, y_max + 0.7, y_min),\n  z     = c(floor_z, floor_z, ceil_z + 0.6),\n  label = c(\"x (m)\", \"y (m)\", \"E (kV/m)\")\n)\naxis_labels$px <- project_x(axis_labels$x, axis_labels$y, axis_labels$z)\naxis_labels$py <- project_y(axis_labels$x, axis_labels$y, axis_labels$z)\n\n# --- Plot ---------------------------------------------------------------\np <- ggplot() +\n  geom_polygon(\n    data = surface_df, aes(px, py, group = poly_id, fill = z_mid),\n    color = PAGE_BG, linewidth = 0.05\n  ) +\n  geom_path(\n    data = base_lines, aes(px, py, group = line_id),\n    color = INK_SOFT, linewidth = 0.35, alpha = 0.45, linetype = \"22\"\n  ) +\n  geom_path(\n    data = surface_lines, aes(px, py, group = line_id),\n    color = INK, linewidth = 0.3, alpha = 0.7\n  ) +\n  geom_segment(\n    data = axis_lines, aes(x = px, y = py, xend = pxend, yend = pyend),\n    color = INK_SOFT, linewidth = 0.6\n  ) +\n  geom_text(data = ticks, aes(px, py, label = label), color = INK_SOFT, size = 2.9) +\n  geom_text(data = z_ticks, aes(px, py, label = label), color = INK_SOFT, size = 2.9) +\n  geom_text(data = axis_labels, aes(px, py, label = label), color = INK, size = 3.2, fontface = \"bold\") +\n  scale_fill_steps(\n    low = \"#009E73\", high = \"#4467A3\", n.breaks = 7, name = \"E (kV/m)\"\n  ) +\n  labs(title = \"contour-3d · r · ggplot2 · anyplot.ai\") +\n  coord_fixed(ratio = 1, clip = \"off\") +\n  theme_void(base_size = 8) +\n  theme(\n    plot.background   = element_rect(fill = PAGE_BG, color = PAGE_BG),\n    panel.background  = element_rect(fill = PAGE_BG, color = NA),\n    plot.title        = element_text(color = INK, size = 12, hjust = 0.5, margin = margin(b = 14)),\n    plot.margin       = margin(t = 20, r = 30, b = 10, l = 20),\n    legend.background = element_rect(fill = PAGE_BG, color = NA),\n    legend.text       = element_text(color = INK_SOFT, size = 7.5),\n    legend.title      = element_text(color = INK, size = 9)\n  )\n\n# --- Save -----------------------------------------------------------------\nggsave(\n  filename = sprintf(\"plot-%s.png\", THEME),\n  plot     = p,\n  device   = ragg::agg_png,\n  width    = 8,\n  height   = 4.5,\n  units    = \"in\",\n  dpi      = 400\n)\n"}