% zero-dimensional, steady state energy balance model 
% for a blackbody planet 

clear

sigma=5.67e-8;      % Stefan-Boltzmann constant, W/m^2/K^4
epsilon=1;          % emissivity for blackbody object
Se=1370;            % solar constant W/m^2, Earth=1370

re=1;               % distance from sun to Earth in AU
rp=[0.387 0.723 re 1.524 9.529 30.087];  % list of planets

albedo=[0.12 0.65 0.31 0.15 0.47 0.41];  % mean albedo of planet surface

rps=['ro';'go';'co';'bo';'mo';'r+'];
rls=['Mercury';'Venus  ';'Earth  ';'Mars   ';'Saturn ';'Neptune'];

Tnot=273.15;        % convert between K and degrees C 

% zero-D steady-state blackbody temperature 
T=(0.25*(Se*re^2./rp.^2).*(1-albedo)/epsilon/sigma).^(1/4)


figure(1)
clf
axis([0 max(rp) min(T)-5 max(T)+5])
hold on
for n=1:length(T)
    plot(rp(n), T(n), rps(n,:))
end
ylabel('temperature (K)')
xlabel('distance from Sun (AU)')
title('blackbody temperatures of solar system objects')
legend(rls)