% Prediction Evaluator
% MATLAB script to predict centennial climate cc
% as an analytic function cc(y) of year y
% scaled according to data up to Presentyear (provided externally)ed
% CONSTANTS and VECTORS in UPPER case, variables and functions in lower case
% Date: December 6, 2024

if exist('PRSNT') == 0; PRSNT = 2024; end     % If not provided, 2024
Px = PRSNT - 1849; Psig = (PRSNT - 2024)/105; % Index, sigma position                               % Presentyear index

% Externally provided time series for HadCRUT5, CO2, and SF
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:PRSNT)';                   % Common data years
cY  = cf(Y) + mean(Y);                  % Centennial years
HC = HadCRUT5(1:Px) + .5;               % + .5 to undo anomalies
cHC = cf(HC);                           % Centennial HadCRUT5
IRB = ce(log2(CO2(1:Px)/280));          % Convert CO2 to radiative blocking
SR = ce(TSI(1:Px));                     % Convert TSI to solar radiation
mdl = [IRB SR];                         % Model components
cmdl = [cf(IRB) cf(SR)];                % Centennial model components
TR = cmdl\cHC;                          % Transient responses [Rg; Rs]
Rg = TR(1);                             % Transient greenhouse response
Rs = TR(2);                             % Transient solar response
cGW = Rg*cf(IRB);                       % Centennial greenhouse warming
cSW = Rs*cf(SR);                        % Centennial solar warming
cC = cGW + cSW;                         % Centennial climate (observed)

% Create Figure 1 for the 4 plots
if exist('f1') == 0
    f1 = figure('Position', [640 400 800 310*175/Px]); 
end
figure(1); clf; hold on

% Upper half of Figure 1
plotc(cY, [cHC cC - .04] - cC(1))        	% Plot 1
% subtracted .04 to make cHC visible
sigma = 1000*std(cHC - cC);
text(cY(end), 0.65 + Psig, sprintf(' \x03c3 = %.1f mK', sigma), 'Rotation', 15)   	% Std dev of residual 1 (mK)
plotc(Y, [HC mdl*TR + mean(HC)] - .5)       % Plot 2

% ====Formal definition of centennial climate====
irb = @(y) log2(1 + 1.0202.^(y - 2058));    % industrial radiative blocking
fcc = @(y) Rg*irb(y);                      % Formal centennial climate
% NB. TSR only used to improve TRR, not used in formal climate.
cFCC = cf(fcc(Y));                          % Centennially filtered fcc
% NB. cFCC is for comparison with cHC: both need to be equally filtered.

% Lower half of Figure 1
plotc(cY, [cHC  cFCC - .04] - 1 - cC(1))    % Plot 3
% subtracted .04 to make cHC visible
sigmap = 1000*std(cHC - cFCC);              % std dev of residual 3 (mK)
text(cY(end), -.35 + Psig, sprintf(' \x03c3 = %.1f mK', sigmap), 'Rotation', 15)      
plotc(Y, [HC fcc(Y)] - 1.5)                 % Plot 4

% ==== Axes ====
ylabel('Plot #            ')
set(gca, 'XLim', Y([1 end]), 'XTick', 1860:20:2020, ...
    'YLim', [-1.78 1.1], 'YTick', (1:4)/2 - 2, ...
    'YTickLabel', [{'(iv)'}; {'(iii)'}; {'(ii)'}; {'(i)'}])

% ==== Scale ====
line(Y(end - 6)*[1 1], [-1.73 -0.73], 'LineWidth', 5)
text(Y(end - 8), -1.4, '1^oC', 'HorizontalAlignment', 'right')

% ==== Legend ====
legend('HadCRUT5', 'R_g*Greenhouse Blocking', 'Location','northwest', 'AutoUpdate', 'off');
text(Y(end - 14 ), -1.6, sprintf('Present year = %d', PRSNT), ...
    'HorizontalAlignment','right')

% ==== Save plots in Figure 1 as a jpeg ====
saveas(1, sprintf('cent%d.jpg', PRSNT))

% Function to plot data and model
function plotc(t, dm)
h = plot(t, dm);
[h.Color] = deal('b', 'r');
end
