function [y,varargout] = func_cp_2state_rev_eq(p,v,T,varargin)
DH=p(1);Th=p(2);
for i=1:length(T)
    cp = func_cp(T{i},DH,Th);
    y{i} = cp; 
end
varargout = {[]};
if isempty(varargin)
    return;
end
for i=1:length(v)
    K = func_K(T{i},DH,Th);
    X=K./(K+1); 
    x{i}=X;
end
varargout = {x};
f = @func_cp; 
for i=1:length(v)
    s = integral(@(x)f(x,DH,Th),T{i}(1),T{i}(end)); % use row vector for T internally
    str{i} = sprintf('Cp integral = %g', s);
end
f = @func_Tm;
Tm = fzero(@(x)f(x,DH,Th),Th);
str{end+1} = fprintf('Tm = %g\n', Tm);
outfun=varargin{1}{1};
hed=varargin{1}{2};
outfun(hed,str);
end
   

function y = func_cp(T,DH,Th)
R = 1.987e-3;  % kcal/mol*K
K = func_K(T,DH,Th);
tmp = K./(K+1);
cp = tmp./(K+1)*DH^2/R./(273+T)./(273+T);
y = cp*1000;
end

function K = func_K(T,DH,Th)
R = 1.987e-3;  % kcal/mol*K
K = exp(-DH/R*(1./(273+T)-1/(273+Th)));
end

function y = func_Tm(T,DH,Th)
R = 1.987e-3;  % kcal/mol*K
K = func_K(T,DH,Th);
y = (1-K)./(1+K)*DH/R./(273+T) - 2;
end