31 lines
1.5 KiB
Julia
31 lines
1.5 KiB
Julia
# Calculate the temperature of the earth using the simplest model
|
||
# @jake
|
||
# https://discourse.julialang.org/t/seven-lines-of-julia-examples-sought/50416/121
|
||
|
||
using Unitful, Plots
|
||
p_sun = 386e24u"W" # power output of the sun
|
||
radius_a = 6378u"km" # semi-major axis of the earth
|
||
radius_b = 6357u"km" # semi-minor axis of the earth
|
||
orbit_a = 149.6e6u"km" # distance from the sun to earth
|
||
orbit_e = 0.017 # eccentricity of r = a(1-e^2)/(1+ecos(θ)) & time ≈ 365.25 * θ / 360 where θ is in degrees
|
||
a = 0.75 # absorptivity of the sun's radiation
|
||
e = 0.6 # emmissivity of the earth (very dependent on cloud cover)
|
||
σ = 5.6703e-8u"W*m^-2*K^-4" # Stefan-Boltzman constant
|
||
temp_sky = 3u"K" # sky temperature
|
||
|
||
|
||
|
||
t = (0:0.25:365.25)u"d" # day of year in 1/4 day increments
|
||
θ = 2*π/365.25u"d" .* t # approximate angle around the sun
|
||
r = orbit_a * (1-orbit_e^2) ./ (1 .+ orbit_e .* cos.(θ)) # distance from sun to earth
|
||
area_projected = π * radius_a * radius_b # area of earth facing the sun
|
||
ec = sqrt(1-radius_b^2/radius_a^2) # eccentricity of earth
|
||
|
||
area_surface = 2*π*radius_a^2*(1 + radius_b^2/(ec*radius_b^2)*atanh(ec)) # surface area of the earth
|
||
|
||
q_in = p_sun * a * area_projected ./ (4 * π .* r.^2) # total heat impacting the earth
|
||
|
||
temp_earth = (q_in ./ (e*σ*area_surface) .+ temp_sky^4).^0.25 # temperature of the earth
|
||
|
||
plot(t*u"d^-1", temp_earth*u"K^-1" .- 273.15, label = false, title = "Temperature of Earth", xlabel = "Day", ylabel = "Temperature [C]")
|