Modern Thermodynamics with Julia in CoCalc

import Pkg; Pkg.add("PlotThemes")

using Plots, Printf, Statistics, LinearAlgebra
using PlotThemes
theme(:bright)

# Stefan–Boltzmann constant (W·m⁻²·K⁻⁴)
const σ_SB = 5.670374419e-8

println("🌍 MODERN THERMODYNAMICS: ENERGY & CLIMATE")
println(repeat("=", 55))
println("Exploring renewable energy and climate thermodynamics")
println()
println("Energy Constants:")
println("• Stefan-Boltzmann: σ = $(σ_SB) W/(m²·K⁴)")

using Printf

println("=== RENEWABLE ENERGY THERMODYNAMICS ===")

# Solar energy analysis
function solar_cell_power(area, irradiance, efficiency)
    # Solar cell power output: P = Area × Irradiance × Efficiency
    return area * irradiance * efficiency
end

function wind_power(air_density, area, wind_speed, efficiency)
    # Wind turbine power: P = 0.5 ρ A v^3 η
    return 0.5 * air_density * area * wind_speed^3 * efficiency
end

# Solar photovoltaic analysis
println("Solar Photovoltaic Systems:")
solar_data = [
    ("Residential (5 kW)", 25, 1000, 0.20),
    ("Commercial (100 kW)", 500, 1000, 0.22),
    ("Utility scale (10 MW)", 50000, 1000, 0.24),
    ("Concentrated PV", 100, 2000, 0.40)
]

println(repeat("-", 70))
@printf "%-20s %10s %10s %10s %12s\n" "System Type" "Area (m²)" "Irrad. (W/m²)" "Efficiency" "Power (kW)"
println(repeat("-", 70))

for (system, area, irradiance, efficiency) in solar_data
    power = solar_cell_power(area, irradiance, efficiency) / 1000  # Convert to kW
    @printf "%-20s %10.0f %10.0f %10.1f%% %12.1f\n" system area irradiance (efficiency*100) power
end

println()

# Wind energy analysis
println("Wind Energy Systems:")
wind_data = [
    ("Small residential", 1.2, 5, 0.35),
    ("Medium commercial", 1.2, 8, 0.45),
    ("Large offshore", 1.2, 12, 0.50),
    ("Offshore wind farm", 1.2, 15, 0.50)
]

println(repeat("-", 80))
@printf "%-20s %12s %12s %10s %12s\n" "System Type" "Air Density" "Wind Speed" "Efficiency" "Power/m²"
println(repeat("-", 80))

for (system, air_density, wind_speed, efficiency) in wind_data
    power_per_area = wind_power(air_density, 1.0, wind_speed, efficiency)
    @printf "%-20s %12.1f %12.0f %10.1f%% %12.0f\n" system air_density wind_speed (efficiency*100) power_per_area
end

println()
println("Renewable Energy Insights:")
println("• Solar: Power scales with area and efficiency")
println("• Wind: Power scales with wind speed cubed (v^3)!")
println("• Efficiency improvements have dramatic impact")
println("• Geographic location crucial (sun/wind resources)")
println("• Storage needed for intermittent sources")
println()

# Betz limit demonstration
betz_limit = 16/27  # ≈ 0.593
println("Wind Energy Theoretical Limits:")
@printf "• Betz limit: %.1f%% (maximum extractable kinetic energy)\n" (betz_limit * 100)
println("• Modern turbines achieve ~75-85% of Betz limit")
println("• Demonstrates fundamental thermodynamic constraints")
println()

# julia
using Printf

# Stefan–Boltzmann constant (W·m⁻²·K⁻⁴)
const σ_SB = 5.670374419e-8

# Climate thermodynamics
"""Earth's equilibrium temperature: T = [(1-α)S/(4σ)]^(1/4)"""
function earth_temperature(solar_constant, albedo)
    return ((1 - albedo) * solar_constant / (4 * σ_SB))^(1/4)
end

"""Greenhouse warming effect in Kelvin"""
function greenhouse_effect(T_no_atmosphere, T_with_atmosphere)
    return T_with_atmosphere - T_no_atmosphere
end

println("=== CLIMATE THERMODYNAMICS ===")
println("Earth's Energy Balance and Greenhouse Effect")
println()

# Earth energy balance
S = 1361  # W/m² (solar constant)
albedo_earth = 0.30  # Earth's average albedo
albedo_ice = 0.90    # Ice/snow albedo
albedo_ocean = 0.06  # Ocean albedo
albedo_forest = 0.15 # Forest albedo

T_no_atmosphere = earth_temperature(S, albedo_earth)
T_current = 288.15  # K (current global average)
greenhouse_warming = greenhouse_effect(T_no_atmosphere, T_current)

println("Earth Energy Balance Analysis:")
@printf "• Solar constant: S = %.0f W/m²\n" S
@printf "• Earth's albedo: α = %.2f\n" albedo_earth
@printf "• Temperature without atmosphere: %.1f K = %.1f°C\n" T_no_atmosphere (T_no_atmosphere - 273.15)
@printf "• Current global average: %.1f K = %.1f°C\n" T_current (T_current - 273.15)
@printf "• Greenhouse effect: +%.1f K (+%.1f°C)\n" greenhouse_warming greenhouse_warming
println()

# Albedo feedback effects
println("Albedo Feedback in Climate System:")
surfaces = [
    ("Fresh snow", 0.85),
    ("Sea ice", 0.70),
    ("Desert sand", 0.35),
    ("Forest", 0.15),
    ("Grassland", 0.20),
    ("Ocean", 0.06)
]

println(repeat("-", 60))
@printf "%-15s %10s %15s %15s\n" "Surface Type" "Albedo" "T_equilibrium (K)" "T_equilibrium (°C)"
println(repeat("-", 60))

for (surface, albedo) in surfaces
    T_eq = earth_temperature(S, albedo)
    T_eq_celsius = T_eq - 273.15
    @printf "%-15s %10.2f %15.1f %15.1f\n" surface albedo T_eq T_eq_celsius
end

println()
println("Climate Feedback Insights:")
println("• Ice-albedo feedback: Melting ice → lower albedo → more warming")
println("• Forest fires increase albedo temporarily")
println("• Ocean absorption crucial for energy balance")
println("• Clouds have complex albedo and greenhouse effects")
println("• Small albedo changes cause significant temperature changes")
println()

# CO₂ and temperature relationship (simplified)
println("Greenhouse Gas Effects:")
co2_preindustrial = 280  # ppm
co2_current = 420       # ppm
co2_doubled = 560       # ppm

# Rough climate sensitivity: ΔT ≈ 3°C per CO₂ doubling
climate_sensitivity = 3.0  # K per CO₂ doubling
warming_current = climate_sensitivity * log(co2_current / co2_preindustrial) / log(2)
warming_doubled = climate_sensitivity * log(co2_doubled / co2_preindustrial) / log(2)

@printf "• Pre-industrial CO₂: %.0f ppm\n" co2_preindustrial
@printf "• Current CO₂: %.0f ppm\n" co2_current
@printf "• Current warming: +%.1f°C (from CO₂ alone)\n" warming_current
@printf "• Doubled CO₂ warming: +%.1f°C\n" warming_doubled
println("• Climate sensitivity: fundamental thermodynamic property")
println("• Logarithmic relationship: equal % increases → equal warming")
println()