% MATLAB script to compute Table 1 and Figure 1 for the paper
% "A one-page account of man-made global warming"
% Vaughan Pratt
% Stanford University
% May 11, 2020

% Function to fit a trend line to a given region of data
trend = @(data, region)[region region.^0]*polyfit(region, data(region), 1)';

% Script to load all data
loadYCRUTEM4CO2EMIT1851to2019

% Some useful constants
Y0 = Y(1) - 1;                              % Year = index + Y0
split8 =  Y(1:(Y(end)-Y(1))/8:end);         % Split period into 8
midy = split8(1+(end-1)*3/4);               % 3/4 year (1977)
AtmInExamoles = 5140/28.97;                 % Earth's atmosphere (examoles)
ppmtoeg = (AtmInExamoles/10^6)*44;          % 44 g/mol of CO2

% 168 years of annual data for emissions, GHGs, and land temp. (^oF)
emeg = emi*44/12/10^6;                      % CO2 emissions in exagrams
ghg = CO2*ppmtoeg + 12.9;                   % Convert CO2 ppm to Eg, add H2O
landTF = CRUTEM4*9/5;                       % Convert Celsius to Fahrenheit

spl = 168*3/4 + 1;                          % 3/4 split index (1977)
r0 = (1:spl)';                              % 1851-1977
r1 = (spl+1:169)';                          % 1978-2019
hiatr = (2001:2012)' - Y0;                  % hiatus range
tr0 = trend(landTF, r0);
sl0 = slope(landTF, r0);
tr1 = trend(landTF, r1);
sl1 = slope(landTF, r1);
trh = trend(landTF, hiatr);

% =====COMPUTE TABLE 1=====

dat = [Y cumsum(emeg) ghg [tr0;tr1]];           % 168x4 table
dat3 = round(dat([1 1977-Y0 end], :)',2);       % 3x4 table (1851, 1977, 2019)
tabl = zeros(4, 7);                             % Prepare 4x7 table
for i = 1:4                                     % Compute table
    tabl(i, 1) = dat3(i, 1);
    for j = 1:2
        Delta(i, j) = dat3(i, j+1) - dat3(i, j);
        tabl(i, j*3-1) = Delta(i, j);
        tabl(i, j*3)   = Delta(i, j)/Delta(1, j)*100;
        tabl(i, j*3+1) = dat3(i, j+1);
    end
end
tabl2 = round(tabl, 2);
display(tabl2)                                          % Show table

% =====PLOT FIGURE 1=====
% For paint order, assumes ax(i).SortMethod == 'childorder' (not 'depth')
% Important to preserve order of painting, and note line before legend.

close all                                       % Start fresh with Fig 1
figure('Position', [10 10 900 435]); clf
hold on
[ax, h1, h2] = plotyy(Y, ghg, Y, landTF);       % Plot GHGs and land temp.
set(h1, 'LineWidth', 3)                         % Thicken GHG for visibility

% All annotations except text boxes and legend
set(ax, 'XTick', split8, 'XLim', [1851 2019])   % Common x-axis
xlabel('\bfFigure 1. Greenhouse gases and land temperatures: 1851-1977-2019');
set(ax(1), 'YLim', [15 16.2], 'YTick', 15:0.3:16.2)
ylabel(ax(1), 'Greenhouse gases: H2O + CO2 (Eg)')
set(ax(2), 'YLim', [-1.5 2.5], 'YTick', -1.5:2.5, 'Color', 'w')
ylabel(ax(2), 'Land temperature anomalies (^oF)')

text(1855, 15.69, [                             % State sources in figure
    {['                                                              ' ...
    'SOURCES: Emissions, CO2, Temp.']}
    {'https://cdiac.ess-dive.lbl.gov/ftp/ndp030/global.1751\_2014.ems'}
    {'http://scrippsco2.ucsd.edu/assets/data/atmospheric/mer'}
    {'ged\_ice\_core\_mlo\_spo/spline\_merged\_ice\_core\_yearly.csv'}
    {'https://crudata.uea.ac.uk/cru/data/temperature/CRUTEM4-gl.dat3'}
], 'FontSize', 12, 'BackgroundColor', 'w')
text(2017, 15.24, 'CO2 CONVERSIONS:', 'FontSize', 12, ...
    'HorizontalAlignment', 'right')
text(2017, 15.14, [
    {'1 Eg CO2 = 28.97/44 * 1000/5.15 ppm = 128 ppm'}
    {'Hence 280 ppm CO2 = 280/128 = 2.19 Eg CO2'}
    ], 'FontSize', 12, 'HorizontalAlignment', 'right', 'BackgroundColor', 'w')

% Three vertical boundaries (light gray)
for y = split8(3:2:end-2)
    line(ax(2), [y y] + .5, [-1.5 2.5], 'Color', 0.85*[1 1 1], ...
        'LineWidth', 2, 'HandleVisibility', 'Off')
end
% Redraw temp. over these boundaries
line(ax(2), Y, landTF, 'Color', 'r', 'HandleVisibility', 'Off')

% Plot the three trend lines
line(ax(2), r0+Y0, tr0, 'LineStyle', '-', 'LineWidth', 3)
line(ax(2), r1+Y0, tr1, 'LineStyle', '-', 'LineWidth', 7)
line(ax(2), hiatr+Y0, trh, 'LineStyle', '-', 'LineWidth', 3, 'Color', 'k')

set(ax(1), 'Color', 'none');   % Tell legend to bring left axis to top (!)
legend('GHGs: H2O + CO2 (Eg)', 'Land temp: CRUTEM4 (^oF)', ...
    sprintf('Trend 1851-1977: %.2f ^oF/cy', tabl2(4, 3)), ...
    sprintf('Trend 1977-2019: %.2f ^oF/cy', tabl2(4, 6)), ...
    'Trend 2001-2012 (Hiatus)', 'Location', 'NorthWest', 'Color', 'w')

% Uncomment if you want to save a new version
% saveas(1, 'Fig1.jpg')
