{"spec_id":"stereonet-equal-area","library":"makie","language":"julia","code":"# anyplot.ai\n# stereonet-equal-area: Structural Geology Stereonet (Equal-Area Projection)\n# Library: makie 0.22.10 | Julia 1.11.9\n# Quality: 89/100 | Created: 2026-06-16\n\nusing CairoMakie\nusing Colors\nusing ColorSchemes\nusing Random\nusing Statistics\n\nRandom.seed!(42)\n\n# --- Theme tokens -----------------------------------------------------------\nconst THEME       = get(ENV, \"ANYPLOT_THEME\", \"light\")\nconst PAGE_BG     = THEME == \"light\" ? colorant\"#FAF8F1\" : colorant\"#1A1A17\"\nconst ELEVATED_BG = THEME == \"light\" ? colorant\"#FFFDF6\" : colorant\"#242420\"\nconst INK         = THEME == \"light\" ? colorant\"#1A1A17\" : colorant\"#F0EFE8\"\nconst INK_SOFT    = THEME == \"light\" ? colorant\"#4A4A44\" : colorant\"#B8B7B0\"\nconst IMPRINT     = [\n    colorant\"#009E73\",  # 1 — bedding (anyplot brand green)\n    colorant\"#C475FD\",  # 2 — joint set 1 (lavender)\n    colorant\"#4467A3\",  # 3 — joint set 2 (blue)\n]\n# Kamb density contours render as a single neutral ink-derived layer (varying\n# alpha by level) so the density encoding reads as distinct from the colored\n# feature poles rather than competing with the green/blue categorical hues.\nconst DENSITY_CMAP = cgrad([\n    RGBAf(INK_SOFT.r, INK_SOFT.g, INK_SOFT.b, 0.28f0),\n    RGBAf(INK_SOFT.r, INK_SOFT.g, INK_SOFT.b, 0.92f0),\n])\n\n# --- Equal-area projection (Schmidt net, lower hemisphere) ------------------\n# A downward unit vector (E, N, D) with D >= 0 maps onto the unit disk via the\n# Lambert azimuthal equal-area projection: plot_x = E / sqrt(1+D), etc.\nproj_x(E, N, D) = E / sqrt(1 + D)\nproj_y(E, N, D) = N / sqrt(1 + D)\n\n# Down-positive unit vector for a line given trend (azimuth, CW from N) + plunge\nline_vec(trend, plunge) =\n    (sind(trend) * cosd(plunge), cosd(trend) * cosd(plunge), sind(plunge))\n\n# Great-circle trace for a plane given dip direction + dip (lower hemisphere)\nfunction great_circle(dip_dir, dip)\n    strike = dip_dir - 90.0\n    u1 = (sind(strike), cosd(strike), 0.0)                       # horizontal strike line\n    u2 = line_vec(dip_dir, dip)                                  # down-dip vector\n    φ  = LinRange(0.0, π, 240)\n    xs = [proj_x(cos(t) * u1[1] + sin(t) * u2[1],\n                 cos(t) * u1[2] + sin(t) * u2[2],\n                 cos(t) * u1[3] + sin(t) * u2[3]) for t in φ]\n    ys = [proj_y(cos(t) * u1[1] + sin(t) * u2[1],\n                 cos(t) * u1[2] + sin(t) * u2[2],\n                 cos(t) * u1[3] + sin(t) * u2[3]) for t in φ]\n    return xs, ys\nend\n\n# --- Data: field measurements from a geological mapping campaign ------------\n# Three structural fabric elements, each clustered around a mean orientation.\nset_names = [\"Bedding\", \"Joint set 1\", \"Joint set 2\"]\nmean_dd   = [110.0, 262.0, 18.0]    # mean dip direction (deg)\nmean_dip  = [34.0, 80.0, 61.0]      # mean dip (deg)\nn_obs     = [18, 15, 16]\ndd_spread = [10.0, 9.0, 8.0]\ndip_spread = [5.0, 4.0, 5.0]\n\ndip_dirs = [Float64[] for _ in set_names]\ndips     = [Float64[] for _ in set_names]\nfor s in 1:length(set_names)\n    for _ in 1:n_obs[s]\n        push!(dip_dirs[s], mod(mean_dd[s] + dd_spread[s] * randn(), 360.0))\n        push!(dips[s], clamp(mean_dip[s] + dip_spread[s] * randn(), 2.0, 88.0))\n    end\nend\n\n# Poles to planes (normals): plunge = 90 - dip, trend = dip_dir + 180\npole_xy  = [Tuple{Float64,Float64}[] for _ in set_names]\npole_vec = NTuple{3,Float64}[]      # flat list of all pole vectors for density\nfor s in 1:length(set_names)\n    for k in 1:n_obs[s]\n        v = line_vec(dip_dirs[s][k] + 180.0, 90.0 - dips[s][k])\n        push!(pole_xy[s], (proj_x(v...), proj_y(v...)))\n        push!(pole_vec, v)\n    end\nend\n\n# --- Kamb-style density field over the projection disk ----------------------\n# Each grid node is back-projected to a downward unit vector; density is the\n# sum of a smoothing kernel over all poles (counts preferred orientations).\nng   = 200\ngrid = LinRange(-1.0, 1.0, ng)\nksm  = 45.0                          # counting-cone concentration\ndensity = fill(NaN, ng, ng)\nfor i in 1:ng, j in 1:ng\n    gx, gy = grid[i], grid[j]\n    ρ2 = gx^2 + gy^2\n    ρ2 > 1.0 && continue\n    D = 1.0 - ρ2                      # inverse equal-area projection\n    f = sqrt(2.0 - ρ2)\n    gE, gN, gD = gx * f, gy * f, D\n    acc = 0.0\n    for (pE, pN, pD) in pole_vec\n        acc += exp(ksm * (gE * pE + gN * pN + gD * pD - 1.0))\n    end\n    density[i, j] = acc\nend\ndmax    = maximum(filter(!isnan, density))\ndensity ./= dmax\nlevels  = [0.15, 0.30, 0.50, 0.70, 0.90]\n\n# --- Figure (square — the stereonet has no preferred horizontal axis) -------\nfig = Figure(size = (1200, 1200), fontsize = 14, backgroundcolor = PAGE_BG)\n\nax = Axis(\n    fig[1, 1];\n    title                = \"stereonet-equal-area · julia · makie · anyplot.ai\",\n    titlesize            = 26,\n    titlecolor           = INK,\n    titlegap             = 16,\n    aspect               = DataAspect(),\n    backgroundcolor      = PAGE_BG,\n    xgridvisible         = false,\n    ygridvisible         = false,\n    xticksvisible        = false,\n    yticksvisible        = false,\n    xticklabelsvisible   = false,\n    yticklabelsvisible   = false,\n    leftspinevisible     = false,\n    rightspinevisible    = false,\n    topspinevisible      = false,\n    bottomspinevisible   = false,\n)\nlimits!(ax, -1.20, 1.20, -1.22, 1.20)\n\n# Elevated disk so the net sits on a subtly distinct surface\ndisk_θ = LinRange(0, 2π, 256)\npoly!(ax, Point2f.(cos.(disk_θ), sin.(disk_θ)); color = ELEVATED_BG, strokewidth = 0)\n\n# Net graticule: meridians + small circles every 10° (subtle reference grid)\nnet_col = RGBAf(INK.r, INK.g, INK.b, 0.10f0)\nfor λ in -80.0:10.0:80.0                 # meridians (great circles about N–S axis)\n    φ = LinRange(-90.0, 90.0, 180)\n    xs = [proj_x(cosd(p) * sind(λ), sind(p), cosd(p) * cosd(λ)) for p in φ]\n    ys = [proj_y(cosd(p) * sind(λ), sind(p), cosd(p) * cosd(λ)) for p in φ]\n    lines!(ax, xs, ys; color = net_col, linewidth = 0.8)\nend\nfor φ0 in -80.0:10.0:80.0                # parallels (small circles)\n    λ = LinRange(-90.0, 90.0, 180)\n    xs = [proj_x(cosd(φ0) * sind(l), sind(φ0), cosd(φ0) * cosd(l)) for l in λ]\n    ys = [proj_y(cosd(φ0) * sind(l), sind(φ0), cosd(φ0) * cosd(l)) for l in λ]\n    lines!(ax, xs, ys; color = net_col, linewidth = 0.8)\nend\n\n# Kamb density contours — highlight preferred pole orientations\ncontour!(ax, grid, grid, density;\n    levels = levels, colormap = DENSITY_CMAP, linewidth = 2.2)\n\n# Mean great circle per fabric element (representative plane)\nfor s in 1:length(set_names)\n    xs, ys = great_circle(mean(dip_dirs[s]), mean(dips[s]))\n    lines!(ax, xs, ys; color = IMPRINT[s], linewidth = 3.4)\nend\n\n# Poles to planes, colored by fabric element (one labeled scatter per set)\nfor s in 1:length(set_names)\n    px = [p[1] for p in pole_xy[s]]\n    py = [p[2] for p in pole_xy[s]]\n    scatter!(ax, px, py;\n        color = IMPRINT[s], markersize = 15, strokecolor = PAGE_BG,\n        strokewidth = 1.4, label = set_names[s])\nend\n\n# Primitive circle (horizontal plane / perimeter)\nlines!(ax, cos.(disk_θ), sin.(disk_θ); color = INK_SOFT, linewidth = 2.4)\n\n# Perimeter degree ticks every 10°; cardinal + azimuth labels\nfor az in 0.0:10.0:350.0\n    major = az % 30 == 0\n    r2 = major ? 1.038 : 1.020\n    lines!(ax, [sind(az), sind(az) * r2], [cosd(az), cosd(az) * r2];\n        color = INK_SOFT, linewidth = major ? 1.6 : 0.9)\nend\nfor az in 0.0:30.0:330.0\n    cardinal = az == 0 ? \"N\" : az == 90 ? \"E\" : az == 180 ? \"S\" : az == 270 ? \"W\" : nothing\n    if cardinal === nothing\n        text!(ax, sind(az) * 1.105, cosd(az) * 1.105;\n            text = string(Int(az)), fontsize = 13, color = INK_SOFT,\n            align = (:center, :center))\n    else\n        text!(ax, sind(az) * 1.10, cosd(az) * 1.10;\n            text = cardinal, fontsize = 22, color = INK, font = :bold,\n            align = (:center, :center))\n    end\nend\n\n# North arrow above the primitive circle\npoly!(ax, Point2f[(-0.035, 1.155), (0.035, 1.155), (0.0, 1.235)];\n    color = INK, strokewidth = 0)\n\naxislegend(ax, \"Lower-hemisphere equal area\";\n    position        = :rb,\n    backgroundcolor = ELEVATED_BG,\n    labelcolor      = INK,\n    titlecolor      = INK,\n    framecolor      = INK_SOFT,\n    framewidth      = 0.8,\n    labelsize       = 13,\n    titlesize       = 13,\n    patchsize       = (18, 18),\n    padding         = (8, 8, 6, 6),\n    margin          = (6, 6, 6, 6),\n    rowgap          = 2,\n)\n\n# --- Save -------------------------------------------------------------------\nsave(\"plot-$(THEME).png\", fig; px_per_unit = 2)\n"}