clear;

hh2o=-285.8;
gh2o=-237.2;
hco2=-393.5;
gco2=-394.4;
%hh2so4=-801.14;
%gh2so4=-793.06;

hso2=-296.81;
gso2=-414.44;

% methane, ethane, propane, butane, pentane, hexane, heptane, octane,
% nonane, decane, undecane, dodecane
h=[-74.8 -84.7 -104.5 -126.5 -173.2 -198.6 -224 -250 -274 -300.9 -327.2 -350.9];

% 2-methyl propane, 2-methyl butane, 2-methyl pentane, 2-methyl hexane,
% 2-methyl heptane, 2,2-dimethyl propane
h=[h -134.5 -178.9 -204.6 -229.5 -255 -189.8];

% cyclopropane, cyclopentane, cyclohexane
h=[h 53.3 -107.1 -156.3];

% ethene, propene, but-1-ene, trans-but-2-ene, cis-but-2-ene,
% buta-1,2-diene, buta-1,3-diene, phenylethene
h=[h 52.5 20.2 -.4 -12.2 -7.8 162.3 109.9 103.8];

% ethyne, propyne
h=[h 228 186.6];

% benzene, methylbenzene, ethylbenzene, propylbenzene, 1,2-dimethylbenzene,
% 1,3-dimethylbenzene, 1,4-dimethylbenzene, ethenylbenzene
h=[h 49 12.1 -13.1 -38.3 -24.4 -25.4 -24.4 103.8];

% methyl amine, dimethyl amine, trimethyl amine, 1-aminobutane
h=[h -23 -18.5 -23.7 -127.6];

% methanol, ethanol, propan-1-ol, propan-2-ol, pentan-1-ol, hexan-1-ol,
% heptan-1-ol, octan-1-ol, ethan-1,2-diol, cyclohexanol
h=[h -239.1 -277.1 -302.7 -317.9 -327.4 -353.6 -377.8 -403.3 -426.6 -454.8 -348.8];

% methoxymethane, ethoxyethane
h=[h -184 -279];

% methanal, ethanal, propanal, butanal
h=[h -108.7 -191.5 -217.1 -241.2];

% acetone, butanone
h=[h -248 -275.7];
h=[h -425 -484.5 -510.7]; % methanoic acid, ethanoic acid, Propanoic acid
h=[h -1003 -527.5 -733]; % glutamic acid, glycine, serine
h=[h -393.5]; % CO2
h=[h -332.9]; % urea
h=[h 254.5 12.4 -165 101.2]; % cyclooctratetraene, nitrobenzene, phenol, pyridine
h=[h -829.94]; % oxalic acid
h=[h -384.9]; % benzoic acid
h=[h -618 -446.4 -776.3 -646.8 -640.6]; % valine, phenylalanine, threonine, leucine, isoleucine
h=[h -678.7 -507.6 -685.6]; % lysine, proline, tyrosine
h=[h -577.5 -534.08 -973.32 -789.02]; % methionine, cysteine, aspartic acid, asparagine
h=[h -1273.3 -2221.2 -1265.6 -2233.08]; % glucose, sucrose, fructose, lactose


g=[-50.8 -32.9 -23.4 -15.6 -9.2 -4.2 1.3 6.4 11.9 17.4 22.8 28.4];
g=[g -17.8 -14.5 -8.1 -2 3.8 -15.2];
g=[g 104.1 36.5 26.8];
g=[g 68.2 74.7 72 62.9 65.9 201.5 151.9 202.5];
g=[g 209.2 194.2];
g=[g 124.5 110.6 19.7 123.8 110.6 107.8 110.3 202.5];
g=[g 32.1 59.2 76.7 -81.8];
g=[g -166.4 -174.9 -171.3 -180.3 -168.9 -163.3 -160.0 -150.0 -136.4 -323.2 -1134.2];
g=[g -114.1 -122.7];
g=[g -113 -128.2 -142.1 -306.4];
g=[g -154 -156];
g=[g -361.4 -389.9 -383.5]; % methanoic acid, ethanoic acid, Propanoic acid
g=[g -476.25 -233.33 -364.17]; % glutamic acid, glycine, serine
g=[g -394.4]; % CO2
g=[g -196.8]; % urea
g=[g 358.1 141.6 -47.5 181.2]; % cyclooctratetraene, nitrobenzene, phenol, pyridine
g=[g -565.52]; % oxalic acid
g=[g -245.1]; % benzoic acid
g=[g -212.95 -64.42 -358.19 -204.53 -204.53]; % valine, phenylalanine, threonine, leucine, isoleucine
g=[g -135.64 -150.26 -219.04]; % lysine, proline, tyrosine
g=[g -177.39 -197.96 -484.67 -281.4]; % methionine, cysteine, aspartic acid, asparagine
g=[g -793.74 -1320.1 -820.7 -1334.42]; % glucose, sucrose, fructose, lactose



n=[1 2 3 4 5 6 7  8 9 10 11 12 4 5 6 7 8 5 3 5 6 2 3 4 4 4 4 4 8 2 3 6 7 8 9 8 8 8 8 1 2 3 4];
n=[n 1 2 3 3 4 5 6 7 8 2 6 ];
n=[n 2 4 ];
n=[n 1 2 3 4];
n=[n 3 4];
n=[n 1 2 3]; 
n=[n 5 2 3 1]; % glutamic acid, glycine, serine, CO2
n=[n 1]; % urea
n=[n 8 6  6  5]; % cyclooctratetraene, nitrobenzene, phenol, pyridine
n=[n 2]; % oxalic acid
n=[n 7]; % benzoic acid
n=[n 5 9 4 6 6]; %valine, phenylalanine, threonine, leucine, isoleucine
n=[n 6 5 9]; % lysine, proline, tyrosine
n=[n 5 3 4 4]; % methionine, cysteine, aspartic acid, asparagine
n=[n 6 12 6 12]; % glucose, sucrose, fructose, lactose

m=[2 3 4 5 6 7 8 9 10 11 12 13 5 6 7 8 9 6 3 5 6 2 3 4 4 4 3 3 4 1 2 3 4 5 6 5 5 5 4 2.5 3.5 4.5 5.5];
m=[m 2 3 4 4 5 6 7 8 9 3 6];
m=[m 3 5];
m=[m 1 2 3 4];
m=[m 3 4];
m=[m 1 2 3];
m=[m 4.5];
m=[m 2.5 3.5 0];
m=[m 2]; % urea
m=[m 4 2.5 3 2.5]; % cyclooctratetraene, nitrobenzene, phenol, pyridine
m=[m 1]; % oxalic acid
m=[m 3]; % benzoic acid
m=[m 5.5 5.5 4.5 6.5 6.5]; % valine, phenylalanine, threonine, leucine, isoleucine
m=[m 7 4.5 5.5]; % lysine, proline, tyrosine
m=[m 5.5 3.5 3.5 4]; % methionine, cysteine, aspartic acid, asparagine
m=[m 6 11 6 11]; % glucose, sucrose, fructose, lactose

N=[zeros(1,39) 1 1 1 1 zeros(1,22) 1 1 1 0];
N=[N 2]; % urea
N=[N 0 1 0 1]; % cyclooctratetraene, nitrobenzene, phenol, pyridine
N=[N 0]; % oxalic acid
N=[N 0]; % benzoic acid
N=[N 1 1 1 1 1]; % valine, phenylalanine, threonine, leucine, isoleucine 
N=[N 2 1 1]; % lysine, proline, tyrosine
N=[N 1 1 1 2]; % methionine, cysteine, aspartic acid, asparagine
N=[N 0 0 0 0]; % glucose, sucrose, fructose, lactose

O=[zeros(1,47) 1 1 1 1 2 1 1 1 1]; 
O=[O 1 1 1 1];
O=[O 1 1];
O=[O 2 2 2];
O=[O 4];
O=[O 2 3 2];
O=[O 1]; % urea
O=[O 0 2 1 0]; % cyclooctratetraene, nitrobenzene, phenol, pyridine
O=[O 4]; % oxalic acid
O=[O 2]; % benzoic acid
O=[O 2 2 3 2 2]; % valine, phenylalanine, threonine, leucine, isoleucine
O=[O 2 2 3]; % lysine, proline, tyrosine
O=[O 1 2 4 3]; % methionine, cysteine, aspartic acid, asparagine
O=[O 6 11 6 11]; % glucose, sucrose, fructose, lactose

S=zeros(1,84);
S=[S 1 1 0 0]; % methionine, cysteine, aspartic acid, asparagine
S=[S 0 0 0 0]; % glucose, sucrose, fructose, lactose


on=(-2*m+N*3+2*O+2*S)./n;


hc=n*hco2+m*hh2o-h+S*hso2;
gc=n*gco2+m*gh2o-g+S*gso2;
format short g;
disp([[1:92]'  hc' gc' on' (hc./gc)' (hc./(12*n))']);


p=polyfit(on,gc./(12*n),1);
yhat=polyval(p,[-4 4]);

%subplot (1,3,1);
plot(on,gc./(12*n),'k^',[-4 4],yhat,'k');
set(gca,'xlim',[-4 4.2],'ylim',[-90 5]);
xlabel('Carbon oxidation number','fontsize',12);
ylabel('Gibbs free energy of combustion (J/mgC)','fontsize',12);

disp('correlation between oxidation number and Gibbs energy');
[r,p]=corr([on' (gc./(12*n))'])

pause;

%subplot(1,3,2);


loglog(-hc./(12*n),-gc./(12*n),'^k',[9 100],[9 100],'k',[9 100],0.9*[9 100],'k--',[9 100],1.1*[9 100],'k--');
set(gca,'fontsize',12,'xlim',[9 100],'ylim',[9 100]);
xlabel('enthalpy of combustion (J/mgC)','fontsize',12);
ylabel('Gibbs free energy of combustion (J/mgC)','fontsize',12);

pause;

T=O+N+n+m;
T=n;
y=gc./(12*n);
%y=hc./(12*n);
X=[ones(1,length(on)); m./T; N./T; O./T; S./T]';
[B,BINT,R,RINT,STATS]=regress(y',X);
yhat=X*B;
disp('statistics for figure 3');
disp([B BINT]);
disp(STATS(1));

%subplot(1,3,3);


plot(y',yhat,'k^',[-100 0],[-100 0],'k');
xlabel('Gibbs free energy of combustion (J/mgC)','fontsize',12);
ylabel('Calculated Gibbs free energy of combustion (J/mgC)','fontsize',12);
pause;
clf;


chon=[2.76	6.33	4.78	0.58
2.36	5.61	5.40	0.51
2.39	5.49	3.94	0.50
2.31	5.26	3.75	0.46
1.92	4.48	3.52	0.38
1.80	4.16	3.17	0.35
1.90	4.45	3.39	0.36
1.46	3.37	2.18	0.26
1.48	3.39	2.38	0.27
0.99	2.35	1.74	0.18
1.03	2.45	1.91	0.17
0.69	1.69	1.71	0.11];

H=chon(:,2)./chon(:,1)/2;
O=chon(:,3)./chon(:,1);
N=chon(:,4)./chon(:,1);
S=zeros(12,1);

X=[ones(12,1) H N O S];


yhat=X*B;
disp(yhat);
disp(length(m));

xxx=abs(2*(gc-hc)./(gc+hc));
I=find(xxx<=0.1);
disp('fraction of H and G data within 10% of each other');
disp(length(I)/length(m));


% Gibbs energies for particulate matter in traps

depth=[100 110 120:10:160 175 200 250 300 500];
G=[-31 -24 -32 -33 -30 -31 -31 -35 -33 -31 -30 -21];
xmat=[100 141  176 301 501];
xmat=[100 141  176];
%xmat=[100 176];
k=length(xmat)-1;
for j=1:k;
    I=find(depth>=xmat(j) & depth<xmat(j+1));
    Z(j)=mean(depth(I));
    Y(j)=mean(G(I));
    YS(j)=std(G(I))/sqrt(length(I));
end;
clf;
errorbar(Z,Y,YS,'k^');
hold on;
plot(depth(9:12),G(9:12),'k^');
set(gca,'fontsize',12);
xlabel('Depth (m)','fontsize',12);
ylabel('standard Gibbs energy of combustion (J/mgC)','fontsize',12);
hold off;

[r,p]=corr([depth(9:12)' G(9:12)'],'type','pearson')
