{"spec_id":"skewt-logp-atmospheric","library":"makie","language":"julia","code":"# anyplot.ai\n# skewt-logp-atmospheric: Skew-T Log-P Atmospheric Diagram\n# Library: makie 0.22.10 | Julia 1.11.9\n# Quality: 85/100 | Created: 2026-05-22\n\nusing CairoMakie\nusing Colors\nusing Random\n\nRandom.seed!(42)\n\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\",\n    colorant\"#C475FD\",\n    colorant\"#4467A3\",\n    colorant\"#BD8233\",\n    colorant\"#AE3030\",\n    colorant\"#2ABCCD\",\n    colorant\"#954477\",\n]\n\n# Skew-T parameters and thermodynamic constants\nconst SKEW     = 45.0\nconst Lv_CONST = 2.501e6\nconst Rd_CONST = 287.0\nconst Rv_CONST = 461.5\nconst Cp_CONST = 1005.0\n\n# Coordinate helpers\nskew_x(T_C, P_hPa) = T_C + SKEW * log10(1000.0 / P_hPa)\ny_p(P_hPa)         = log10(P_hPa)\n\n# Saturation vapor pressure via Bolton (1980), hPa\nsat_es(T_C) = 6.112 * exp(17.67 * T_C / (T_C + 243.5))\n\n# Saturation mixing ratio, g/kg\nws_sat(T_C, P_hPa) = 622.0 * sat_es(T_C) / (P_hPa - sat_es(T_C))\n\n# Log-spaced pressure array from 1000 → 100 hPa\nconst P_FINE = collect(exp10.(LinRange(log10(1000.0), log10(100.0), 300)))\n\n# Tropical convective sounding (standard atmosphere with instability)\nconst P_OBS = [1000.0, 950.0, 925.0, 900.0, 850.0, 800.0, 750.0, 700.0,\n               650.0, 600.0, 550.0, 500.0, 450.0, 400.0, 350.0, 300.0,\n               250.0, 200.0, 150.0, 100.0]\n\nconst T_OBS = [28.4,  25.6,  23.8,  21.6,  17.2,  12.4,   8.0,   3.2,\n               -1.8,  -7.2, -12.8, -19.2, -25.6, -33.2, -41.6, -50.4,\n               -59.2, -65.8, -68.4, -72.6]\n\nconst TD_OBS = [24.8,  21.4,  19.6,  16.8,  10.4,   5.2,  -0.8,  -8.4,\n               -15.2, -22.6, -30.2, -38.4, -46.8, -52.6, -58.8, -63.2,\n               -68.4, -72.0, -74.2, -77.6]\n\n# Reference line helpers\n\nfunction isotherm_line(T0, p_arr)\n    return [skew_x(T0, P) for P in p_arr], y_p.(p_arr)\nend\n\nfunction dry_adiabat_line(theta_K, p_arr)\n    Rcp = Rd_CONST / Cp_CONST\n    return [skew_x(theta_K * (P / 1000.0)^Rcp - 273.15, P) for P in p_arr], y_p.(p_arr)\nend\n\nfunction mixing_ratio_line(ws_gkg, p_arr)\n    xs = map(p_arr) do P\n        es_t = ws_gkg * P / (622.0 + ws_gkg)\n        es_t <= 0.0 && return NaN\n        a = log(es_t / 6.112)\n        skew_x(243.5 * a / (17.67 - a), P)\n    end\n    return xs, y_p.(p_arr)\nend\n\nfunction moist_adiabat_line(T0_C, P0_hPa, p_arr)\n    xs  = Float64[]\n    ys  = Float64[]\n    T_K = T0_C + 273.15\n    for i in 1:length(p_arr)\n        P = p_arr[i]\n        P > P0_hPa + 0.5 && continue\n        push!(xs, skew_x(T_K - 273.15, P))\n        push!(ys, y_p(P))\n        if i < length(p_arr)\n            ws    = max(0.0, ws_sat(T_K - 273.15, P) / 1000.0)\n            dT_dp = (Rd_CONST * T_K + Lv_CONST * ws) /\n                    (P * (Cp_CONST + Lv_CONST^2 * ws / (Rv_CONST * T_K^2)))\n            T_K   = T_K + dT_dp * (p_arr[i + 1] - P)\n        end\n    end\n    return xs, ys\nend\n\n# Compute lifted parcel temperatures (Bolton 1980) and return LCL metadata.\n# p_arr must be sorted descending (1000 → 100 hPa).\nfunction lifted_parcel_temps(T_surf_C, Td_surf_C, P_surf_hPa, p_arr)\n    T_K  = T_surf_C + 273.15\n    Td_K = Td_surf_C + 273.15\n    Rcp  = Rd_CONST / Cp_CONST\n\n    # LCL temperature and pressure (Bolton 1980)\n    T_LCL_K = 56.0 + 1.0 / (1.0 / (Td_K - 56.0) + log(T_K / Td_K) / 800.0)\n    P_LCL   = P_surf_hPa * (T_LCL_K / T_K)^(Cp_CONST / Rd_CONST)\n\n    T_moist = T_LCL_K\n    prev_P  = P_LCL\n    T_out   = Float64[]\n\n    for P in p_arr\n        if P >= P_LCL\n            # Dry adiabatic lifting\n            push!(T_out, T_K * (P / P_surf_hPa)^Rcp - 273.15)\n        else\n            # Moist adiabatic lifting — one Euler step from prev level\n            ws      = max(0.0, ws_sat(T_moist - 273.15, prev_P) / 1000.0)\n            dT_dp   = (Rd_CONST * T_moist + Lv_CONST * ws) /\n                      (prev_P * (Cp_CONST + Lv_CONST^2 * ws / (Rv_CONST * T_moist^2)))\n            T_moist = T_moist + dT_dp * (P - prev_P)\n            prev_P  = P\n            push!(T_out, T_moist - 273.15)\n        end\n    end\n    return T_out, T_LCL_K - 273.15, P_LCL\nend\n\n# Parcel path at both resolutions\nparcel_Ts_obs, T_LCL_C, P_LCL_hPa = lifted_parcel_temps(\n    T_OBS[1], TD_OBS[1], P_OBS[1], P_OBS)\nparcel_Ts_fine, _, _ = lifted_parcel_temps(\n    T_OBS[1], TD_OBS[1], P_OBS[1], P_FINE)\n\n# Figure\nfig = Figure(\n    size            = (1600, 900),\n    fontsize        = 14,\n    backgroundcolor = PAGE_BG,\n)\n\nax = Axis(\n    fig[1, 1];\n    title              = \"skewt-logp-atmospheric · julia · makie · anyplot.ai\",\n    titlesize          = 20,\n    titlecolor         = INK,\n    xlabel             = \"Temperature (°C)\",\n    ylabel             = \"Pressure (hPa)\",\n    xlabelsize         = 14,\n    ylabelsize         = 14,\n    xticklabelsize     = 12,\n    yticklabelsize     = 12,\n    xlabelcolor        = INK,\n    ylabelcolor        = INK,\n    xticklabelcolor    = INK_SOFT,\n    yticklabelcolor    = INK_SOFT,\n    xtickcolor         = INK_SOFT,\n    ytickcolor         = INK_SOFT,\n    backgroundcolor    = PAGE_BG,\n    leftspinecolor     = INK_SOFT,\n    bottomspinecolor   = INK_SOFT,\n    topspinevisible    = false,\n    rightspinevisible  = false,\n    xgridvisible       = false,\n    ygridvisible       = false,\n    yreversed          = true,\n)\n\n# Y-axis: log10(P) with 1000 hPa at bottom, 100 hPa at top\nylims!(ax, y_p(100.0) - 0.02, y_p(1000.0) + 0.02)\nconst P_YTICKS = [1000.0, 925.0, 850.0, 700.0, 500.0, 400.0, 300.0, 200.0, 100.0]\nax.yticks = (y_p.(P_YTICKS), string.(Int.(P_YTICKS)))\n\n# X-axis: skewed temperature coordinate, ticks at surface temps\nxlims!(ax, -47.0, 88.0)\nconst T_XTICKS = Float64[-40, -30, -20, -10, 0, 10, 20, 30, 40]\nax.xticks = (T_XTICKS, string.(Int.(T_XTICKS)) .* \"°\")\n\n# Background isobars\nisobar_col = RGBAf(INK.r, INK.g, INK.b, 0.08f0)\nfor P in [950.0, 925.0, 850.0, 800.0, 750.0, 700.0, 650.0, 600.0, 550.0,\n          500.0, 450.0, 400.0, 350.0, 300.0, 250.0, 200.0, 150.0]\n    hlines!(ax, y_p(P); color = isobar_col, linewidth = 0.6)\nend\n\n# Isotherms (every 10 °C)\niso_col = THEME == \"light\" ?\n    RGBAf(0.55f0, 0.55f0, 0.55f0, 0.28f0) :\n    RGBAf(0.62f0, 0.62f0, 0.62f0, 0.22f0)\n\nfor T0 in -60.0:10.0:60.0\n    xs, ys = isotherm_line(T0, P_FINE)\n    lines!(ax, xs, ys; color = iso_col, linewidth = 0.7)\n    x_bot = skew_x(T0, 1000.0)\n    if -47.0 <= x_bot <= 88.0\n        text!(ax, x_bot, y_p(1000.0) - 0.014;\n              text    = \"$(Int(T0))°\",\n              fontsize = 9,\n              color   = INK_SOFT,\n              align   = (:center, :top))\n    end\nend\n\n# Dry adiabats (potential temperature from -40 to 80 °C, every 10 °C)\ndry_col = THEME == \"light\" ?\n    RGBAf(0.84f0, 0.37f0, 0.0f0, 0.40f0) :\n    RGBAf(0.90f0, 0.55f0, 0.20f0, 0.33f0)\n\nfor (i, theta_C) in enumerate(-40.0:10.0:80.0)\n    xs, ys = dry_adiabat_line(theta_C + 273.15, P_FINE)\n    lbl    = i == 1 ? \"Dry adiabat\" : nothing\n    if isnothing(lbl)\n        lines!(ax, xs, ys; color = dry_col, linewidth = 0.8, linestyle = :dash)\n    else\n        lines!(ax, xs, ys; color = dry_col, linewidth = 0.8, linestyle = :dash, label = lbl)\n    end\nend\n\n# Moist adiabats (surface start temps -5 to 35 °C, every 5 °C)\nmoist_col = THEME == \"light\" ?\n    RGBAf(0.0f0, 0.447f0, 0.698f0, 0.42f0) :\n    RGBAf(0.25f0, 0.62f0, 0.87f0, 0.35f0)\n\nfor (i, T_start) in enumerate(-5.0:5.0:35.0)\n    xs, ys = moist_adiabat_line(T_start, 1000.0, P_FINE)\n    lbl    = i == 1 ? \"Moist adiabat\" : nothing\n    if isnothing(lbl)\n        lines!(ax, xs, ys; color = moist_col, linewidth = 0.8, linestyle = :dashdot)\n    else\n        lines!(ax, xs, ys; color = moist_col, linewidth = 0.8, linestyle = :dashdot, label = lbl)\n    end\nend\n\n# Mixing ratio lines (g/kg), displayed only below 600 hPa\np_low   = P_FINE[P_FINE .>= 599.0]\nmix_col = THEME == \"light\" ?\n    RGBAf(0.0f0, 0.62f0, 0.45f0, 0.48f0) :\n    RGBAf(0.0f0, 0.75f0, 0.55f0, 0.40f0)\n\nfor (i, ws_gkg) in enumerate([2.0, 4.0, 8.0, 12.0, 20.0])\n    xs, ys = mixing_ratio_line(ws_gkg, p_low)\n    valid  = .!isnan.(xs) .& isfinite.(xs)\n    lbl    = i == 1 ? \"Mixing ratio\" : nothing\n    if any(valid)\n        if isnothing(lbl)\n            lines!(ax, xs[valid], ys[valid]; color = mix_col, linewidth = 0.7, linestyle = :dot)\n        else\n            lines!(ax, xs[valid], ys[valid]; color = mix_col, linewidth = 0.7, linestyle = :dot, label = lbl)\n        end\n        # Label at top of visible section (600 hPa)\n        x_top = xs[valid][end]\n        y_top = ys[valid][end]\n        if -47.0 <= x_top <= 88.0\n            text!(ax, x_top, y_top;\n                  text    = \"$(Int(ws_gkg))\",\n                  fontsize = 8,\n                  color   = mix_col,\n                  align   = (:center, :bottom))\n        end\n    end\nend\n\n# CAPE region — poly! fills the closed polygon between parcel and environment.\n# This uses Makie's native polygon primitive for efficient filled-area rendering.\ncape_mask = parcel_Ts_obs .> T_OBS\nif any(cape_mask)\n    ci          = findall(cape_mask)\n    env_xs_c    = [skew_x(T_OBS[i],        P_OBS[i]) for i in ci]\n    par_xs_c    = [skew_x(parcel_Ts_obs[i], P_OBS[i]) for i in ci]\n    ys_c        = y_p.(P_OBS[ci])\n    cape_pts    = Point2f.(\n        vcat(par_xs_c, reverse(env_xs_c)),\n        vcat(ys_c,     reverse(ys_c)),\n    )\n    cape_fill = RGBAf(IMPRINT[5].r, IMPRINT[5].g, IMPRINT[5].b, 0.22f0)\n    poly!(ax, cape_pts; color = cape_fill, strokewidth = 0, label = \"CAPE\")\nend\n\n# Lifted parcel trace (smooth, fine resolution)\nparcel_xs_fine = [skew_x(T, P) for (T, P) in zip(parcel_Ts_fine, P_FINE)]\nlines!(ax, parcel_xs_fine, y_p.(P_FINE);\n       color = IMPRINT[4], linewidth = 2.0, linestyle = :dashdotdot,\n       label = \"Lifted parcel\")\n\n# LCL marker\nlcl_x = skew_x(T_LCL_C, P_LCL_hPa)\nlcl_y = y_p(P_LCL_hPa)\nscatter!(ax, [lcl_x], [lcl_y];\n         color = IMPRINT[4], markersize = 10, marker = :diamond, strokewidth = 0)\ntext!(ax, lcl_x + 1.5, lcl_y;\n      text    = \"LCL\",\n      fontsize = 9,\n      color   = IMPRINT[4],\n      align   = (:left, :center))\n\n# Sounding profiles\nT_xs   = [skew_x(T, P) for (T, P) in zip(T_OBS, P_OBS)]\nTD_xs  = [skew_x(Td, P) for (Td, P) in zip(TD_OBS, P_OBS)]\nobs_ys = y_p.(P_OBS)\n\nlines!(ax, T_xs, obs_ys;\n       color = IMPRINT[1], linewidth = 3.5, label = \"Temperature\")\nscatter!(ax, T_xs, obs_ys;\n         color = IMPRINT[1], markersize = 7, strokewidth = 0)\n\nlines!(ax, TD_xs, obs_ys;\n       color = IMPRINT[2], linewidth = 3.5, linestyle = :dash, label = \"Dewpoint\")\nscatter!(ax, TD_xs, obs_ys;\n         color = IMPRINT[2], markersize = 7, strokewidth = 0)\n\naxislegend(ax;\n    position        = :rt,\n    backgroundcolor = ELEVATED_BG,\n    labelcolor      = INK,\n    framecolor      = INK_SOFT,\n    framewidth      = 0.8,\n    labelsize       = 12,\n)\n\nsave(\"plot-$(THEME).png\", fig; px_per_unit = 2)\n"}