function [y,varargout] = func_cp_2state_irrev(p,v,T,varargin)
DHu=p(1);Eu=p(2);Tu=p(3);
for i=1:length(v)
    cp = func_cp(T{i},DHu,Eu,Tu,v(i),T{i}(1));
    y{i} = cp; 
end
varargout = {[]};
if isempty(varargin)
    return;
end
for i=1:length(v)
    f = @myIntegrand;
    for j=1:length(T{i})
        tmp1(j) = integral(@(x)f(x,Eu,Tu),T{i}(1),T{i}(j));
    end
    tmp2 = exp(-tmp1/v(i));   
    x{i}=1-tmp2;
end
varargout = {x};    
f = @func_cp;
for i=1:length(v)
    s = integral(@(x)f(x,DHu,Eu,Tu,v(i),T{i}(1)),T{i}(1),T{i}(end)); 
    str{i} = sprintf('Cp integral = %g', s);
end
str{end+1}=sprintf('\n');
outfun=varargin{1}{1};
hed=varargin{1}{2};
outfun(hed,str);
end
   

function y = func_cp(T,DHu,Eu,Tu,v,T0)
R = 1.987e-3;  % kcal/mol*K
f = @myIntegrand;
for i=1:length(T)
    tmp1(i) = integral(@(x)f(x,Eu,Tu),T0,T(i));
end
tmp2 = exp(-tmp1/v);    
ku = func_k(T,Eu,Tu);
cp = DHu*ku/v.*tmp2;    
y = cp*1000;   
end


function y = myIntegrand(T,E,Ts)
k = func_k(T,E,Ts);
y = k;
end

function k = func_k(T,E,Ts)
R = 1.987e-3;  % kcal/mol*K
k = exp(-E/R*(1./(T+273)-1/(Ts+273)));
end

