%% Some example plotting code to show the ouput of the FDA and the fossil
% SI inversions. These are based on those submitted with the Geology
% paper but are not formatted specifically for publication. 

clear all
close all

% Read in the initial data used for the FDA
a = csvread('ginkgo_data.csv');
CO2      = a(:,1);
CO2lo   = a(:,2);
SI      = a(:,3);
SIlo = a(:,4);

% Read in the output from the FDA. In this case it contains CO2fine and
% SI_outfine

load simple_out_g_030409.mat

%%
% Calculate the median, 2.5 and 97.5 percentiles of the FDA functions.
% Becuase fineCO2 varies for each function, we have to first bin the data
% into equally spaced bins.

numbbreak=1000;
mmin=min(min(CO2fine));
mmax=max(max(CO2fine));
mmin=min(CO2); mmax=max(CO2);
bbreak=linspace(mmin,mmax,numbbreak);

for i = 1:numbbreak-1
    found=SI_outfine(find(bbreak(i)<CO2fine & CO2fine <bbreak(i+1)));
    true_quart(i,:) = quantile(found,[0.025,0.5,0.975]);
end

true_quart(numbbreak,:)=true_quart(numbbreak-1,:);

%% Plot the FDA functions - median and 95% CI

ylab='SI (%)';
xlab='CO_2 (ppmv)';

figure; % dotted black line CI
[bb]=gcf;
plot(CO2,SI,'ro','markersize',4);hold on
plot(bbreak,true_quart(:,1),':k');hold on
plot(bbreak,true_quart(:,2),'k');hold on
plot(bbreak,true_quart(:,3),':k');hold on
ylabel(ylab);xlabel(xlab);

figure; % solid grey CI
errorbargrey(bbreak,true_quart(:,2),true_quart(:,2)-true_quart(:,1),true_quart(:,2)-true_quart(:,3));hold on;
errorbarxy(CO2,SI,CO2lo,SIlo,[],[],'or','k');hold on
plot(bbreak,true_quart(:,2),'k');hold on
ylabel(ylab);xlabel(xlab);

%%
load CO2_SI_statistics_060409

foss_obs=load('ginkgo_foss_data_22.csv');

[numobs nn] = size(foss_obs);
complete_foss_obs=foss_obs;

%%
gam=4.3;
ws=7.4;
sunef=ones(numobs,6);
rco=mean_std./280;
wst=complete_foss_obs(:,5)./570;
rco(:,6)=complete_foss_obs(:,4)./280;
deltatt=gam.*(log(rco));
for i=1:6;sunef(:,i)=sunef(:,i).*(ws.*wst);end
deltat=deltatt-sunef;


%%
select=[ 1 2 3 4 5 6 7 8];

ixx=kernel_1; fff=kernel_2;

figure;
for j = 1:8
i=select(j);
subplot(4,2,j)
plot(ixx(i,:),fff(i,:),'k')
firstline=['SI to invert: ',num2str(complete_foss_obs(i,2)),' %'];
medianround=round(mean_std(i,4));
secondline=['Median: ',num2str(medianround),'shelagh mcdonald CO_2 (ppmv)'];
title({firstline;secondline}); xlabel(xlab); ylabel('Frequency')
end

figure
subplot(1,2,1)
errorbar(complete_foss_obs(:,5),mean_std(:,4),mean_std(:,4)-mean_std(:,3),mean_std(:,4)-mean_std(:,5),'bo');set(gca,'xdir','reverse');hold on
plot(complete_foss_obs(:,5),complete_foss_obs(:,4),'r*');
ylabel('CO_2 (ppmv)');xlabel('Time (Ma)');
subplot(1,2,2)
errorbar(complete_foss_obs(:,5),deltat(:,4),deltat(:,4)-deltat(:,3),deltat(:,4)-deltat(:,5),'bo');set(gca,'xdir','reverse');hold on
plot(complete_foss_obs(:,5),deltat(:,6),'r*');
ylabel('\Delta T (^oC)');xlabel('Time (Ma)');
