function Ts = temp_gradient2(So, obliquity, eccentricity, perihelion, albedo, opticalthickness)

% Creates a table of values for mean temperature by latitude and plots same.

%

% Ts = temp_gradient2(So, obliquity, eccentricity, perihelion, albedo, opticalthickness)

%

% Where:

% So is the solar constant for the body

% obliquity is the axial tilt of the body

% eccentricity is the eccentricity of the body’s orbit

% perihelion is the longitude of periapsis of the body’s orbit

%

% Ts is an array containing the surface temperatures for each 1 degree

%      latitude band from the south pole to the north pole

% orbital parameters

%eccentricity = 0.0172; % eccentricity of Earth’s orbit

%obliquity = 23.5; % tilt of axis rotation (degrees)

%So = 1365; % solar constant (W/m^2)

%perihelion = 281.37; % longitude of perihelion (degrees longitude)

% other parameters

%albedo = 0.3; % fraction of incident radiation reflected

%frac_IR_absorbed = 0.8; % fraction of IR absorbed

% grid

latitude = -90:90; % 90 South to 90 North

[day lat] = meshgrid(1:5:365, latitude);

% calculate annual mean insolation vs latitude

% see function below

insol = daily_insolation(So,eccentricity,obliquity,perihelion,lat,day);

insol_mean = mean(insol,2);

% emission temperature

sigma = 5.67e-8; % Stefan-Boltzmann constant

Te = (insol_mean*(1-albedo)/sigma).^(1/4); % equation 2-4

% surface temperature

% from leaky greenhouse model

%Ts = (2/(2-frac_IR_absorbed))^(1/4)*Te; % equation 2-13

convection = 30/67;

Ts = ((1 + opticalthickness- convection)*bsxfun(@power,Te,4));

Ts = bsxfun(@power, Ts, .25);

% plot it

figure;

plot(latitude, Ts);

line1 = strcat(‘obliquity – \epsilon: ‘, num2str(obliquity), ‘ eccentricity – e: ‘, num2str(eccentricity));

line2 = strcat(‘Solar constant – So: ‘, num2str(So), ‘W/m^2’, ‘ Longitude of perihelion – \omega: ‘, num2str(perihelion));

title({line1 ; line2});

xlim([-90 90]);

set(gca,’XTick’,-90:15:90)

xlabel(‘latitude’);

ylabel(‘Surface temperature (K)’);

function Fsw = daily_insolation(So,eccentricity,obliquity,perihelion,lat,day)

% based on the code by I. Eisenman

% [day lat]=meshgrid(1:5:365, -90:90);

% Fsw=daily_insolation(eccentricity,obliquity,perihelion,omega,lat,day);

% plot(-90:90,mean(Fsw”))

epsilon=obliquity * pi/180;

omega=perihelion * pi/180;

% === Calculate insolation ===

lat=lat*pi/180; % latitude

% lambda (or solar longitude) is the angular distance along Earth’s orbit measured from spring equinox (21 March)

% estimate lambda from calendar day using an approximation from Berger 1978 section 3

delta_lambda_m=(day-80)*2*pi/365.2422;

beta=(1-eccentricity.^2).^(1/2);

lambda_m0=-2*( (1/2*eccentricity+1/8*eccentricity.^3).*(1+beta).*sin(-omega)-…

1/4*eccentricity.^2.*(1/2+beta).*sin(-2*omega)+1/8*eccentricity.^3.*(1/3+beta).*(sin(-3*omega)) );

lambda_m=lambda_m0+delta_lambda_m;

lambda=lambda_m+(2*eccentricity-1/4*eccentricity.^3).*sin(lambda_m-omega)+…

(5/4)*eccentricity.^2.*sin(2*(lambda_m-omega))+(13/12)*eccentricity.^3.*sin(3*(lambda_m-omega));

%So=1365; % solar constant (W/m^2)

delta=asin(sin(epsilon).*sin(lambda)); % declination of the sun

Ho=acos(-tan(lat).*tan(delta)); % hour angle at sunrise/sunset

% no sunrise or no sunset: Berger 1978 eqn (8),(9)

Ho( ( abs(lat) >= pi/2 – abs(delta) ) & ( lat.*delta > 0 ) )=pi;

Ho( ( abs(lat) >= pi/2 – abs(delta) ) & ( lat.*delta <= 0 ) )=0;

% Insolation: Berger 1978 eq (10)

Fsw=So/pi*(1+eccentricity.*cos(lambda-omega)).^2 ./ …

(1-eccentricity.^2).^2 .* …

( Ho.*sin(lat).*sin(delta) + cos(lat).*cos(delta).*sin(Ho) );