% MATLAB Script for Figures 1 and 2 of my response to Curry 12/23/2019.
% Author: Vaughan R. Pratt
% Date: December 25, 2019
%
% Command count: Setup = 12, Download & format = 7, Plot figures = 12
% Total = 31 commands
% (Plus 4 display commands for progress information)

% SETUP (4 + 4 + 4 = 12)

% Define centennial filter and center functions (4 commands)
box = @(w)(ones(w, 1)/w);		    % Box shape, unit area
K = conv(box(11), box(65));		    % K = Centennial convolution kernel
f = @(d)conv(d, K, 'valid');        % f = Centennial filter
c = @(d)(d - mean(d));              % c = Centering function

% Miscellaneous data (4 commands)
Albedo = 0.3;                       % Earth's shortwave reflectivity
Y = (1850:2017)';                   % Time span for "recent" climate
fY = f(Y);                          % Merely truncates each end of Y
wo = weboptions('Timeout', 15);     % The default of 5 seconds is too short

% Support for Kopp vs. Svalgaard (4 commands)
KS = 1;                             % 1 = Kopp, change to 2 for Svalgaard
KSid = {'Kopp' 'Svalgaard'};        % Identity of reconstructor
ab = 'ab';                          % Figure 2 a (Kopp) b (Svalgaard)
CSS = 1.35 + (KS-1)*0.6; % Centennial solar sensitivity: Kopp 1.35 Sval 1.95

% DOWNLOAD CLIMATE DATA FROM WEB; FILTER; CENTER   (2 + 2 + 3 = 7 commands)
% ALL DATA IS FOR 1850-2017.  
% MATLAB downloads are a bit flaky, so if they fail just run curry.m again.

% 1. GMST (Global Mean Surface Temperature) = Annual HadCRUT4 (2 commands)
display(['HadCRUT4 data from www.metoffice.gov.uk/hadobs/' ...
    'hadcrut4/data/current/time_series/ ...']); % Progress view
HC4li = str2num(cell2mat(splitlines(webread(['https://' ...
    'www.metoffice.gov.uk/hadobs/hadcrut4/data/current/time_series/' ...
    'HadCRUT.4.6.0.0.annual_ns_avg.txt'], wo))));
GMST = c(f(HC4li(Y - 1849, 2))); % filter and center annual global HadCRUT4

% 2. ARF (Anthropogenic radiative forcing) = 1.85*log2(CO2) (2 commands)
display('CO2 data from scrippsco2.ucsd.edu/assets/data/atmospheric ...');
CO2cell = webread(['https://scrippsco2.ucsd.edu/assets/data/atmospheric' ...
    '/merged_ice_core_mlo_spo/spline_merged_ice_core_yearly.csv'], wo);
ARF = c(f(1.85*log2(CO2cell{Y, 2}))); % scale CO2 by TCR, filter, center 

% 3. SF (Solar forcing) = SS*(1 - Albedo)*TSI/4 = SS*7/40*TSI (3 commands)
% Recent TSI reconstructions by Svalgaard and Kopp, from an Excel file)
display(['Svalgaard''s and Kopp''s TSI reconstructions ' ...
    'from leif.org/research/ ...']);
urlwrite('https://leif.org/research/Historical-TSI.xls', 'TSI.xls');
TSImat = xlsread('TSI.xls');        % load .xls file into TSImat (411 x 6)
% delete('TSI.xls');                % uncomment if you don't want to keep it
SF = CSS*c(f((1 - Albedo)*TSImat(Y - 1606, KS + 2)/4));% GK,LS in cols C,D

% PLOT FIGURES    (Each figure: 2 commands to plot, 4 to annotate)

% FIGURE 1
% Plot
figure(1)
plot(fY, [GMST ARF])                % Plot climate vs. anthropogenc RF
% Annotate
set(gca, 'XLim', fY([1 end]), 'XTick', 1870:10:2000)
xlabel('\bfFigure 1. GMST vs. anthropogenic RF')
ylabel('Centered filtered temperature (^oC)')
legend('GMST', 'ARF = 1.85*log_2(CO2)')

% FIGURE 2 (a or b if KS = 1 or 2)
% Plot
figure(2)
plot(fY, [GMST-ARF SF]);            % Plot Fig. 1 residual vs. solar forcing
% Annotate
set(gca, 'XLim', fY([1 end]))
xlabel(['\bfFigure 2' ab(KS) '. Fig. 1 residual vs. ' ...
     KSid{KS} '''s reconstruction of TSI'])
ylabel('Temperature (^oC)')
legend('Residual = GMST - ARF', sprintf('SF = %.2f*(1 - Albedo)*TSI/4', CSS))

display(['Std dev of Fig. 2' ab(KS) ' residual: ' ...
    num2str(1000*std(GMST-ARF - SF), 2) ' mK'])
