for n = 1:10
lwabsorbedatm = n/2;
[Tatm(n),Tsurf(n)]=modelwithatm(lwabsorbedatm);
for n = 1:10
surface_albedo = n/2;
[Tatm(n),Tsurf(n)]=modelwithatm(surface_albedo);
end
end
plot(.2:.2:2, Tatm)
hold on
plot(.2:.2:2, Tsurf)
hold off
legend('ta','ts')
axis([0 2 100 400])
function [Tatm, Tsurf] = modelwithatm(lwabsorbedatm,surface_albedo)
albedo =.31;
Ks = 1361;
Re = 6173;
Sigma = 5.67E-8;
Epsilon = 1;
atm_albedo = .35;
swradiation_sun = .9;
lwradiation_sun = .1;
Epsilon_atm = .9;
syms Tsurf
syms Tatm
a = Ks*.9*pi*Re^2*0.02;
d = Ks*.9*pi*Re^2*(1-.02 -.35)*.2*.02;
c = Ks*.9*pi*Re^2*(1-.02 -.35)*(1-surface_albedo);
q =Ks*.1*pi*Re^2*lwabsorbedatm;
f = Ks*.1*pi*Re^2*(1-.95);
y = 4*pi*Re^2*Epsilon*Sigma*Tsurf^4*lwabsorbedatm;
z = 4*pi*Re^2*Epsilon_atm*Sigma*Tatm^4;
p = 4*pi*Re^2*Epsilon*Sigma*Tsurf^4;
Energy_absorbed_atm = a + d + q + (p*.95);
Energy_absorbed_surf = c + f + z;
Energy_emitted_atm = 2*z;
Energy_emitted_surface = p;
equationa = (Energy_absorbed_atm == Energy_emitted_atm);
equations = (Energy_absorbed_surf == Energy_emitted_surface);
[tempa, temps] = solve([equationa, equations],[Tatm,Tsurf]);
Tsurf = abs(real(double(temps(1))))
Tatm = abs(real(double(tempa(1))))
end