{"spec_id":"skewt-logp-atmospheric","library":"ggplot2","language":"r","code":"#' anyplot.ai\n#' skewt-logp-atmospheric: Skew-T Log-P Atmospheric Diagram\n#' Library: ggplot2 3.5.1 | R 4.4.1\n#' Quality: 90/100 | Created: 2026-05-21\n\nlibrary(ggplot2)\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\"\nIMPRINT   <- c(\"#009E73\", \"#C475FD\", \"#4467A3\", \"#BD8233\",\n                 \"#AE3030\", \"#2ABCCD\", \"#954477\")\n\n# --- Coordinate Helpers ------------------------------------------------------\n# Y-axis: negative log10 of pressure so that high pressure (surface) sits at\n# the bottom and low pressure (upper atm) sits at the top.\nlp <- function(P_hpa) -log10(P_hpa)\n\n# X-axis: temperature + skew offset to tilt isotherms 45°\nSKEW <- 45   # °C per log10(P) decade\nskewx <- function(T_c, P_hpa) T_c + SKEW * log10(1000 / P_hpa)\n\n# --- Atmospheric Physics (Bolton 1980 approximations) -----------------------\nes_hpa <- function(T_c) 6.112 * exp(17.67 * T_c / (T_c + 243.5))\n\nws_kgkg <- function(T_c, P_hpa) {\n  es <- pmin(es_hpa(T_c), P_hpa * 0.99)\n  0.622 * es / (P_hpa - es)\n}\n\n# Pseudo-adiabatic lapse rate (Euler integration, P decreasing = rising air)\nmoist_adiabat <- function(T0_c, P0 = 1000, Pend = 100, n = 200) {\n  Lv <- 2.501e6; Rd <- 287.0; Cp <- 1004.0\n  P_lev <- exp(seq(log(P0), log(Pend), length.out = n))\n  T_c <- numeric(n)\n  T_c[1] <- T0_c\n  for (i in 2:n) {\n    T_K  <- T_c[i - 1] + 273.15\n    P    <- P_lev[i - 1]\n    ws   <- ws_kgkg(T_c[i - 1], P)\n    num  <- Rd * T_K / P + Lv * ws / P\n    den  <- Cp + Lv^2 * ws * 0.622 / (Rd * T_K^2)\n    T_c[i] <- T_c[i - 1] + (num / den) * (P_lev[i] - P)\n  }\n  data.frame(P = P_lev, T_c = T_c)\n}\n\n# --- Radiosonde Sounding (tropical pre-storm profile) -----------------------\npres_obs <- c(1000, 925, 850, 750, 700, 650, 600, 550, 500, 450,\n              400, 350, 300, 250, 200, 150, 100)\ntemp_obs <- c(30, 24, 17, 9, 4, -2, -8, -15, -22, -29,\n              -37, -44, -52, -58, -62, -66, -70)\ndew_obs  <- c(24, 20, 13, 2, -6, -14, -22, -33, -43, -55,\n              -62, -68, -73, -77, -80, -83, -86)\n\n# Smooth interpolation in log-P space\np_fine  <- exp(seq(log(max(pres_obs)), log(min(pres_obs)), length.out = 120))\nT_fine  <- approx(log(pres_obs), temp_obs, xout = log(p_fine))$y\nTd_fine <- approx(log(pres_obs), dew_obs,  xout = log(p_fine))$y\n\nsounding <- data.frame(\n  y   = lp(p_fine),\n  xT  = skewx(T_fine,  p_fine),\n  xTd = skewx(Td_fine, p_fine)\n)\n\n# --- Reference Lines --------------------------------------------------------\np_ref <- exp(seq(log(1050), log(95), length.out = 400))\ny_ref <- lp(p_ref)\n\n# Temperature isotherms — mapped to series for legend\niso_df <- do.call(rbind, lapply(seq(-80, 60, by = 10), function(T0) {\n  data.frame(group = T0, y = y_ref, x = skewx(T0, p_ref), series = \"Isotherms\")\n}))\n\n# Dry adiabats (Poisson: T_K = theta * (P/1000)^0.286)\ndry_df <- do.call(rbind, lapply(seq(265, 390, by = 10), function(th_K) {\n  T_c <- th_K * (p_ref / 1000)^0.286 - 273.15\n  data.frame(group = th_K, y = y_ref, x = skewx(T_c, p_ref), series = \"Dry Adiabats\")\n}))\n\n# Moist pseudo-adiabats starting from the surface\nmoist_df <- do.call(rbind, lapply(c(4, 10, 16, 22, 28, 34), function(T0) {\n  df       <- moist_adiabat(T0, P0 = 1000, Pend = 100, n = 200)\n  df$group <- T0\n  df$y     <- lp(df$P)\n  df$x     <- skewx(df$T_c, df$P)\n  df$series <- \"Moist Adiabats\"\n  df\n}))\n\n# Saturation mixing ratio lines (lower troposphere only)\np_mix  <- exp(seq(log(1050), log(500), length.out = 200))\nmix_df <- do.call(rbind, lapply(c(1, 2, 4, 8, 16), function(w_gpkg) {\n  w_kg <- w_gpkg / 1000\n  es   <- w_kg * p_mix / (0.622 + w_kg)\n  T_c  <- 243.5 * log(es / 6.112) / (17.67 - log(es / 6.112))\n  data.frame(group = w_gpkg, y = lp(p_mix), x = skewx(T_c, p_mix),\n             series = \"Mixing Ratios\")\n}))\n\n# --- Y-axis Setup (log-pressure labels) -------------------------------------\np_labeled <- c(1000, 850, 700, 500, 400, 300, 200, 100)\ny_breaks  <- lp(p_labeled)\ny_labels  <- as.character(p_labeled)\n\n# Isobar positions for horizontal reference lines\ny_isobars <- lp(c(1000, 850, 700, 500, 400, 300, 200, 100))\n\n# Plot y limits (from 1050 hPa at bottom to 95 hPa at top)\ny_lo <- lp(1050)  # ≈ -3.021\ny_hi <- lp(95)    # ≈ -1.978\n\n# --- Long-form sounding for colour + linetype legend -----------------------\nsnd_long <- rbind(\n  data.frame(y = sounding$y, x = sounding$xT,  series = \"Temperature\"),\n  data.frame(y = sounding$y, x = sounding$xTd, series = \"Dewpoint\")\n)\n\n# --- Combined scale values (order controls legend display) ------------------\nseries_order <- c(\"Temperature\", \"Dewpoint\", \"Isotherms\",\n                  \"Dry Adiabats\", \"Moist Adiabats\", \"Mixing Ratios\")\n\ncolor_values <- c(\n  \"Temperature\"   = IMPRINT[1],   # #009E73 green\n  \"Dewpoint\"      = IMPRINT[2],   # #C475FD orange-red\n  \"Isotherms\"     = INK_MUTED,      # gray, theme-adaptive\n  \"Dry Adiabats\"  = IMPRINT[5],   # #AE3030 amber\n  \"Moist Adiabats\"= IMPRINT[4],   # #BD8233 purple (distinct from sky blue)\n  \"Mixing Ratios\" = IMPRINT[6]    # #2ABCCD sky blue\n)\n\nlinetype_values <- c(\n  \"Temperature\"   = \"solid\",\n  \"Dewpoint\"      = \"dashed\",\n  \"Isotherms\"     = \"solid\",\n  \"Dry Adiabats\"  = \"dashed\",\n  \"Moist Adiabats\"= \"dotted\",\n  \"Mixing Ratios\" = \"longdash\"\n)\n\n# --- Plot -------------------------------------------------------------------\np <- ggplot() +\n  # Horizontal isobars at standard levels\n  geom_hline(\n    yintercept = y_isobars,\n    color = INK_MUTED, linewidth = 0.2, alpha = 0.40\n  ) +\n  # Reference lines — mapped to series aesthetic for unified legend\n  geom_path(\n    data = iso_df,\n    aes(x = x, y = y, group = group, color = series, linetype = series),\n    linewidth = 0.28, alpha = 0.55\n  ) +\n  geom_path(\n    data = dry_df,\n    aes(x = x, y = y, group = group, color = series, linetype = series),\n    linewidth = 0.28, alpha = 0.55\n  ) +\n  geom_path(\n    data = moist_df,\n    aes(x = x, y = y, group = group, color = series, linetype = series),\n    linewidth = 0.28, alpha = 0.55\n  ) +\n  geom_path(\n    data = mix_df,\n    aes(x = x, y = y, group = group, color = series, linetype = series),\n    linewidth = 0.28, alpha = 0.55\n  ) +\n  # Main sounding profiles with legend (Temperature = first/brand colour)\n  geom_path(\n    data = snd_long,\n    aes(x = x, y = y, color = series, linetype = series),\n    linewidth = 1.5\n  ) +\n  scale_color_manual(\n    name   = NULL,\n    values = color_values,\n    breaks = series_order\n  ) +\n  scale_linetype_manual(\n    name   = NULL,\n    values = linetype_values,\n    breaks = series_order\n  ) +\n  scale_y_continuous(\n    breaks = y_breaks,\n    labels = y_labels,\n    limits = c(y_lo, y_hi)\n  ) +\n  scale_x_continuous(breaks = seq(-30, 60, by = 10)) +\n  coord_cartesian(xlim = c(-42, 62)) +\n  labs(\n    x     = \"Temperature (°C)\",\n    y     = \"Pressure (hPa)\",\n    title = paste0(\n      \"Tropical Radiosonde Sounding · \",\n      \"skewt-logp-atmospheric · r · ggplot2 · anyplot.ai\"\n    )\n  ) +\n  theme_minimal(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    panel.grid.major  = element_blank(),\n    panel.grid.minor  = element_blank(),\n    panel.border      = element_rect(color = INK_SOFT, fill = NA, linewidth = 0.5),\n    axis.title        = element_text(color = INK,      size = 10),\n    axis.text         = element_text(color = INK_SOFT, size = 8),\n    plot.title        = element_text(color = INK,      size = 11, face = \"bold\"),\n    legend.background = element_rect(fill = ELEVATED_BG, color = INK_SOFT,\n                                     linewidth = 0.3),\n    legend.text       = element_text(color = INK_SOFT, size = 8),\n    legend.position.inside = c(0.88, 0.05),\n    legend.justification   = c(1, 0),\n    plot.margin       = margin(20, 30, 15, 20)\n  ) +\n  guides(\n    color = guide_legend(\n      override.aes = list(\n        linewidth = c(1.5, 1.5, 0.5, 0.5, 0.5, 0.5),\n        alpha     = c(1.0, 1.0, 0.8, 0.8, 0.8, 0.8)\n      )\n    ),\n    linetype = \"none\"\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"}