{"spec_id":"stereonet-equal-area","library":"ggplot2","language":"r","code":"#' anyplot.ai\n#' stereonet-equal-area: Structural Geology Stereonet (Equal-Area Projection)\n#' Library: ggplot2 3.5.1 | R 4.4.1\n#' Quality: 91/100 | Created: 2026-06-16\n\nlibrary(ggplot2)\nlibrary(dplyr)\nlibrary(scales)\nlibrary(ragg)\n\nset.seed(42)\n\n# --- Theme tokens -----------------------------------------------------------\nTHEME       <- Sys.getenv(\"ANYPLOT_THEME\", \"light\")\nPAGE_BG     <- if (THEME == \"light\") \"#FAF8F1\" else \"#1A1A17\"\nELEVATED_BG <- if (THEME == \"light\") \"#FFFDF6\" else \"#242420\"\nINK         <- if (THEME == \"light\") \"#1A1A17\" else \"#F0EFE8\"\nINK_SOFT    <- if (THEME == \"light\") \"#4A4A44\" else \"#B8B7B0\"\nINK_MUTED   <- if (THEME == \"light\") \"#6B6A63\" else \"#A8A79F\"\n\n# Imprint palette — first categorical series is always brand green\nIMPRINT_PALETTE <- c(\"#009E73\", \"#C475FD\", \"#4467A3\", \"#BD8233\",\n                     \"#AE3030\", \"#2ABCCD\", \"#954477\", \"#99B314\")\n\nNET_COL  <- alpha(INK, 0.10)   # faint equal-area net grid\nPRIM_COL <- INK_SOFT           # primitive circle + ticks\nCONT_COL <- alpha(INK, 0.50)   # density contour overlay\n\n# --- Projection helpers (lower-hemisphere Schmidt / equal-area) -------------\n# Equal-area radius of a downward line: r = sqrt(1 - sin(plunge)). North is +y,\n# East is +x, azimuth measured clockwise from North.\ndeg2rad <- function(d) d * pi / 180\n\nproject_line <- function(trend, plunge) {\n  z <- sin(deg2rad(plunge))            # downward direction cosine\n  r <- sqrt(pmax(0, 1 - z))            # equal-area radial distance (primitive = 1)\n  tr <- deg2rad(trend)\n  data.frame(x = r * sin(tr), y = r * cos(tr))\n}\n\n# Project a downward unit vector given as (North, East, Down) components.\nproject_vec <- function(N, E, D) {\n  r <- sqrt(pmax(0, 1 - D))\n  h <- sqrt(N^2 + E^2)\n  s <- ifelse(h > 0, r / h, 0)\n  data.frame(x = E * s, y = N * s)\n}\n\n# Great circle of a plane (right-hand-rule strike, dip), lower hemisphere.\ngreat_circle <- function(strike, dip, n = 160) {\n  S  <- deg2rad(strike)\n  Dr <- deg2rad(dip)\n  # u = horizontal strike line, d = down-dip line; both unit and orthogonal.\n  u <- c(cos(S), sin(S), 0)\n  d <- c(-cos(Dr) * sin(S), cos(Dr) * cos(S), sin(Dr))\n  phi <- seq(0, pi, length.out = n)   # lower-hemisphere half (Down >= 0)\n  N <- cos(phi) * u[1] + sin(phi) * d[1]\n  E <- cos(phi) * u[2] + sin(phi) * d[2]\n  D <- cos(phi) * u[3] + sin(phi) * d[3]\n  project_vec(N, E, D)\n}\n\n# --- Equal-area net grid (meridians + small circles about the E-W axis) -----\nnet_arcs <- list()\ngi <- 0\nlambda <- seq(0, pi, length.out = 120)\nfor (g in seq(-80, 80, 10)) {           # meridians: great circles through E & W\n  gr <- deg2rad(g)\n  N <- sin(lambda) * -sin(gr)\n  E <- cos(lambda)\n  D <- sin(lambda) * cos(gr)\n  gi <- gi + 1\n  net_arcs[[gi]] <- transform(project_vec(N, E, D), grp = gi)\n}\nfor (b in seq(10, 170, 10)) {           # small circles: cones about the E-W axis\n  br <- deg2rad(b)\n  N <- sin(br) * cos(lambda)\n  E <- rep(cos(br), length(lambda))\n  D <- sin(br) * sin(lambda)\n  gi <- gi + 1\n  net_arcs[[gi]] <- transform(project_vec(N, E, D), grp = gi)\n}\nnet <- bind_rows(net_arcs)\n\n# --- Primitive circle, perimeter ticks, cardinal labels ---------------------\nct <- seq(0, 2 * pi, length.out = 400)\nprimitive <- data.frame(x = cos(ct), y = sin(ct))\n\nang   <- seq(0, 350, 10)\nmajor <- ang %% 90 == 0\nticks <- data.frame(\n  x0 = sin(deg2rad(ang)),\n  y0 = cos(deg2rad(ang)),\n  x1 = sin(deg2rad(ang)) * ifelse(major, 1.06, 1.035),\n  y1 = cos(deg2rad(ang)) * ifelse(major, 1.06, 1.035)\n)\n\ncardinal <- data.frame(\n  lab = c(\"E\", \"S\", \"W\"),\n  ang = c(90, 180, 270)\n)\ncardinal$x <- sin(deg2rad(cardinal$ang)) * 1.13\ncardinal$y <- cos(deg2rad(cardinal$ang)) * 1.13\n\n# --- Structural data: bedding, two joint sets, a fault set ------------------\nmake_set <- function(n, strike_mu, strike_sd, dip_mu, dip_sd, label) {\n  data.frame(\n    strike       = rnorm(n, strike_mu, strike_sd) %% 360,\n    dip          = pmin(89, pmax(2, rnorm(n, dip_mu, dip_sd))),\n    feature_type = label\n  )\n}\n\nmeasurements <- bind_rows(\n  make_set(22, 42,  12, 24, 6, \"Bedding\"),\n  make_set(18, 118, 9,  80, 6, \"Joint Set 1\"),\n  make_set(16, 205, 10, 74, 7, \"Joint Set 2\"),\n  make_set(12, 312, 11, 58, 8, \"Fault\")\n)\nfeature_levels <- c(\"Bedding\", \"Joint Set 1\", \"Joint Set 2\", \"Fault\")\nmeasurements$feature_type <- factor(measurements$feature_type, levels = feature_levels)\n\n# Poles to planes: plunge 90 - dip toward (strike - 90)\npoles <- measurements\npole_xy <- project_line((poles$strike + 270) %% 360, 90 - poles$dip)\npoles$x <- pole_xy$x\npoles$y <- pole_xy$y\n\n# Great circles for every plane (coloured by feature type)\nplanes_list <- vector(\"list\", nrow(measurements))\nfor (i in seq_len(nrow(measurements))) {\n  arc <- great_circle(measurements$strike[i], measurements$dip[i])\n  arc$grp <- i\n  arc$feature_type <- measurements$feature_type[i]\n  planes_list[[i]] <- arc\n}\nplanes <- bind_rows(planes_list)\n\n# --- Density field over the projected poles (Kamb-style clustering) ---------\ngx <- seq(-1, 1, length.out = 90)\ngrid <- expand.grid(x = gx, y = gx)\nband <- 0.13\ndens <- numeric(nrow(grid))\nfor (i in seq_len(nrow(poles))) {\n  dens <- dens + exp(-((grid$x - poles$x[i])^2 + (grid$y - poles$y[i])^2) / (2 * band^2))\n}\ngrid$z <- dens\ngrid$z[sqrt(grid$x^2 + grid$y^2) > 1] <- NA   # clip to the primitive circle\n\n# --- Plot -------------------------------------------------------------------\np <- ggplot() +\n  geom_path(data = net, aes(x, y, group = grp),\n            color = NET_COL, linewidth = 0.3) +\n  geom_path(data = primitive, aes(x, y),\n            color = PRIM_COL, linewidth = 0.8) +\n  geom_segment(data = ticks, aes(x = x0, y = y0, xend = x1, yend = y1),\n               color = PRIM_COL, linewidth = 0.5) +\n  geom_path(data = planes, aes(x, y, group = grp, color = feature_type),\n            linewidth = 0.45, alpha = 0.30) +\n  geom_contour(data = grid, aes(x, y, z = z),\n               color = CONT_COL, linewidth = 0.5, bins = 6) +\n  geom_point(data = poles, aes(x, y, color = feature_type),\n             size = 2.8, alpha = 0.95, stroke = 0.4) +\n  annotate(\"text\", x = 0, y = 1.155, label = \"N\",\n           color = INK, fontface = \"bold\", size = 6.5) +\n  geom_text(data = cardinal, aes(x, y, label = lab),\n            color = INK_SOFT, size = 4.6) +\n  scale_color_manual(values = IMPRINT_PALETTE[1:4], name = \"Feature type\") +\n  coord_fixed(xlim = c(-1.2, 1.2), ylim = c(-1.2, 1.24),\n              expand = FALSE, clip = \"off\") +\n  labs(\n    title    = \"stereonet-equal-area · r · ggplot2 · anyplot.ai\",\n    subtitle = \"Lower-hemisphere Schmidt net — poles, planes & pole-density contours\"\n  ) +\n  guides(color = guide_legend(override.aes = list(linewidth = 0, size = 5, alpha = 1))) +\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 = 14, hjust = 0.5,\n                                    margin = margin(b = 4)),\n    plot.subtitle    = element_text(color = INK_SOFT, size = 10, hjust = 0.5,\n                                    margin = margin(b = 6)),\n    legend.position  = \"bottom\",\n    legend.title     = element_text(color = INK, size = 11),\n    legend.text      = element_text(color = INK_SOFT, size = 10),\n    legend.key       = element_rect(fill = PAGE_BG, color = NA),\n    plot.margin      = margin(14, 14, 10, 14)\n  )\n\n# --- Save -------------------------------------------------------------------\nggsave(\n  filename = sprintf(\"plot-%s.png\", THEME),\n  plot     = p,\n  device   = ragg::agg_png,\n  width    = 6,\n  height   = 6,\n  units    = \"in\",\n  dpi      = 400\n)\n"}