% MATLAB script to plot CHC - cGW vs. cSW
% Author: Vaughan Pratt
% Date: November, 2024

load HadCRUT5.csv
load CO2.csv
load TSI.csv

% Filtered HadCRUT5 and forcings
ce = @(D) D - mean(D);                          % Center data D _____
TRPZD = conv(ones(11, 1)/11, ones(65, 1)/65);   % Trapezoid    /_____\
cf = @(d) ce(conv(d, TRPZD, 'Valid'));  % Centennial filter: 11 & 65 yrs

% Transient responses to centennial radiative and solar forcings
Y  = (1850:2024)';           	    % Common data years
cY  = cf(Y) + mean(Y);              % Centennial years
HC = HadCRUT5 + .5;           	    % + .5 to undo anomalies
cHC = cf(HC);                       % Centennial HadCRUT5
GB = log2(CO2/280);           	    % Convert CO2 to Greenhouse Blocking
SR = ce(TSI);                 	    % Convert TSI to Solar Radiation
mdl = [GB SR];                      % Model components
cGB = cf(GB);                       % Centennial greenhouse blocking
cSR = cf(SR);                       % Centennial solar radiation
cmdl = [cGB cSR];                   % Centennial model components
TR = cmdl\cHC;                      % Transient responses [R_g; R_s]
cGW = TR(1)*cGB;                    % Centennial greenhouse warming
cSW = TR(2)*cSR;                    % Centennial solar warming

figure(1); clf
h = plot(cY, [cHC - cGW cSW]);
[h.Color] = deal('b', 'r');
legend('cHC - cGW', sprintf('cSW \x03c3 = %.1f mK', ...
    1000*std((cHC - cGW) - cSW)), 'Location', 'northwest')
ylabel('cSW (^oC)')
set(gca, 'XLim', cY([1 end]), 'YLim', [-1 1]*.03)
saveas(1, 'TSIresid.jpg')
