% MATLAB script to plot 65-year HadCRUT4 with and without TSI
% In EOS October 23, Shaun Lovejoy plotted GISS against CO2 forcing
% Here we plot HadCRUT3 against CO2 forcing after removing essentially all
% other 20th century influences on global climate.  The "fast" contributors
% are the 65-year Atlantic Multidecadal Oscillation and all faster
% contributors, which are removed using a 65-year moving average filter.
% This is the blue curve.  A slower contributor remains, namely the Sun,
% whose Total Solar Irradiance (TSI) rises between 1900 and 1950.
% Subtracting the forcing contribution of TSI from HadCRUT4 yields the red
% curve.  Its straightness indicates that what now remains of HadCRUT4
% grows linearly with CO2 forcing.  This is consistent with the theory that
% rising CO2 forcing results in a rising global climate.  Moreover it
% demonstrates that whatever is causing global temperature to rise is doing
% so in proportion to CO2 forcing.
%
% NB. CO2 is used here only as a proxy for all greenhouse forcing, to which
% CO2 itself many only be contributing about 80%.

% SET UP FIGURE
figure(1); clf; grid on; box on; hold on
w = 65;                                 % Width of filter.
title(sprintf('%d-year HadCRUT4 with and without TSI', w))
xlabel('Forcing = log_2(CO2/280)')
ylabel('Anomaly (^oC)')
set(gca, 'XLim', [.08 .31], 'Ylim', [-.53 .15])

% USE forcing = log2(CO2/280) TO DEFINE THE X-AXIS
load CO2.csv                            % CO2 since 1868
Forcing = log2(CO2/280);                % CO2 forcing in units of doublings
Fw = ma(Forcing, w);                   % w-year CO2 forcing
yrs = 1900:10:1980;                     % Auxiliary x-axis index: years
text(Fw(yrs-1867-floor(w/2)), yrs*0+.149, num2cell(yrs), 'Rotation', -90)

% GET HadCRUT4, SHOW IT AT 5-YR SMOOTHING
load HadCRUT4.csv                       % HadCRUT4 since 1868
h = 11;
Hh = ma(HadCRUT4, h);                   % h-year HadCRUT4
rh = (1:length(Fw))+floor((w - h)/2);      % Index range of Fw in Hh
plot(Fw, Hh(rh), 'Color', .7*[1 1 1])  % Plot it less w-h years

% SMOOTH HadCRUT4, PLOT IT, AND SHOW ITS R2 (as (1-R2)*100%)
Hw = ma(HadCRUT4, w);                  % w-year HadCRUT4 (1900-1982)
plot(Fw, Hw, 'LineWidth', 4, 'Color', 'b') % Plot it
F1 = [Fw Fw*0+1];                     % Append column of 1's to Fw
TrendParamH = polyfit(Fw, Hw, 1)';    % Trend parameters defining
Htrend = F1*TrendParamH;                % the trend line for Hw
H1_R2 = var(Hw - Htrend)/var(Hw);     % Compute 1-R2 for Hw
text(.132, -.21, sprintf('1-R2 = %.2f%%', H1_R2*100), ...
    'Color', 'b');                      %   and show it

% GET TSI, ANOMALIZE, SCALE, AND CONVERT TO DEGREES
load TSI.csv	% Total Solar Irradiance since 1868 (reconstructed)
A = 0.306;      % Earth's bond albedo (NASA Earth fact sheet, DRW, 2004)
L = 1.08;       % Net solar climate sensitivity (best-fit estimate)
TSIanom = TSI - mean(TSI((1961:1990) - 1867));	% TSI anomaly ref 1961-1990
SolForc = TSIanom/4*(1 - A)*L;   % Area of sphere = 4x disk; albedo, sol.c.s.
s = 21;
SolFs = ma(SolForc, s);
rs = (1:length(Fw))+floor((w-s)/2);       % Index range of Fw in SolFs
plot(Fw, SolFs(rs), 'Color', [.9 .6 0])  % Plot it less w-s years

% SUBTRACT SOLAR FORCING FROM 
Nw = ma(HadCRUT4 - SolForc, w);            % w-year HadCRUT4 less TSI
plot(Fw, Nw, 'LineWidth', 4, 'Color', 'r')%   Plot it
TrendParamN = polyfit(Fw, Nw, 1)';        % Trend parameters defining
Ntrend = F1*TrendParamN;                    %   the trend line  for Nw
plot(Fw, Ntrend, 'LineWidth', 2, 'Color', 'g');  % Plot trend
text(.121, -.16, sprintf('1-R2 = %.2f%%', var(Nw-Ntrend)/var(Nw)*100), ...
    'HorizontalAlignment', 'right', 'Color', 'r');% ... and label it
N1_R2 = var(Nw - Ntrend)/var(Nw);                  % Compute 1-R2 for Nw
text(.134, -.13, sprintf('Clim.Sens. = %.2f', TrendParamN(1)), ...
    'HorizontalAlignment', 'right', 'Color', [0 .7 0]);% and show it

% LEGEND: GIVE NAMES TO CURVES
legend(sprintf('HadCRUT4, %d-yr', h), sprintf('HadCRUT4, %d-yr', w), ...
    sprintf('SolForc, %d-yr', s), ...
    sprintf('HadCRUT4 - SolForc, %d-yr', w), ...
    'Trend Line', 'Location', 'SouthEast')

% Parse the H5 rise from 1910 to 1943
r = [1910 1940] - 1867;
rh = r - floor(h/2);
rs = r - floor(s/2);
rw = r - floor(w/2);
riseh = diff(Hh(rh))
rises = diff(SolFs(rs))
risew = diff(Hw(rw))
risec = TrendParamN(1)*diff(Fw(rw))
risec + rises