%% 14.14 Snowfall

depth = load('snowd.dat');
nWeeks = length(depth);

% fit a quadratic to the snowfall data
weeks = 1:nDays;
coeff = polyfit(weeks, depth, 2);
r = roots(coeff);

% the snow gone day is the largest root rounded up
gone = ceil(max(r));

% plot the data and quadratic fit
weeksGone = 1:gone;
plot(weeks, depth, 'o', weeksGone, polyval(coeff, weeksGone), '-');
axis([0 gone 0 max(depth)]);
xlabel('# Weeks');
ylabel('Depth in millimeters');
s = sprintf('Snow gone week %d', gone);
title(s);




%% Vostok data
[age, dt, co2, dust] = textread('vostok.txt', '%d %f %f %f', ...
                                'headerlines', 2, ...
                                'emptyvalue', NaN);
                               
[hAx, hLine1, hLine2] = plotyy(age, dt, age, co2);
title('temperature difference and CO2 levels versus year from present')
xlabel('years from present')
ylabel(hAx(1), 'temperature difference (deg Celcius)');
ylabel(hAx(2), 'C02 level in parts per million');


figure

scatter(co2, dt)
title('temperature difference versus CO2 level')
xlabel('CO2 level in parts per million')
ylabel('temperature difference (deg Celcius)');


coeffs = polyfit(co2, dt, 1);
x = linspace(150, 350, 20);
yfit = polyval(coeffs, x);
hold on
plot(x, yfit, 'r-')
legend('measured values', 'best fit line', 'Location', 'NorthWest')

polyval(coeffs, 400)

y = polyval(coeffs, co2);
res = dt - y;
[X, I] = sort(co2);
plot(X, res(I), '-')
xlabel('CO2 level in parts per million')
ylabel('residual error (deg Celcius)');
title('residual plot after line fitting')

