%this script drives the cpd heat model
%last updated 12/2/2014 by Troy Holland, BYU
clear all;
close all;
clc;

%put in the NMR parameters
p0=0.479;%		!p0		
c0=0.0692;%		!c0
sigp1=4.87;%		!sig+1
mw1=430.3;%		!mw
mdel=43.0;%		!mdel (7 will be subtracted internally to the CPD model)

%Enter the heating rate in K/s. You will likely only want to change "heating0",
%though you are free to change other values if you like.

heating0=5e3;%(K/s)
Tf0=1600;%K
hold0=0;%s

heating1=0;
Tf1=1600;
hold1=.001;


%Input the kinetic parameters. The parameter n is simply a multiplier that
%allows you to change all the values at once by a factor. I suspect you
%won't need it, so it should be left at a value of 1, but it was convenient
%for me.

n=1;

ab=2.6e15*n;%	!ab
eb0=55400*n;%		!eb
ebsig=1800*n;% 		!ebsig
ac=0.9*n;%		!ac=rho
ec0=0*n;%		!ec0
ag=3.e15*n;%		!ag
eg0=60000*n;%		!eg
egsig=8100*n;%     	!egsig
acr=3.e15*n;%		!Acr (pre-exponential factor for crosslinking rate)
ecr=65000*n;%		!Ecr (Activation energy for crosslinking rate)
an=9e17*n;%		!an (Pre-exponential factor for HCN release)
en0=100000*n;%		!en (Activation energy for HCN release)
ensig=17000*n;%		!ensig (deviation bound for dist en for HCN release)


%Build the input matrix. Note that there are additional inputs that
%shouldn't need to be changed in most cases.

Input_matrix=[p0,c0,sigp1,mw1,mdel,heating0,Tf0,hold0,heating1,Tf1,hold1,...
              ab,eb0,ebsig,ac,ec0,ag,eg0,egsig,acr,ecr,an,en0,ensig];


output=cpdheat(Input_matrix);

tar=output(:,1);%unitless fraction
total_volatiles=output(:,2);%unitless fraction
time=output(:,3);%seconds
temperature=output(:,4);%Kelvin



