function [R,navg1,navg2] = twomolmodelanalytic_integrate_hyperbolic(td)
%[R,navg1,navg2] = twomolmodelanalytic_integrate_hyperbolic(td)
%100325
%100616

s1=1;
J1=1.5e21;
n01=2.6e18;
tau1=7.2e-4;
f1=1.7e24;
V1=3.8e-29;
sigma1=4.95e-21;

s2=1;
J2=1.5e19/2;
n02=1.1e18;
tau2=.00009;
f2=f1;
V2=6e-28;
sigma2=2.3e-18;

kd1=s1*J1/n01+1/tau1+sigma1*f1;
kr1=s1*J1/n01+1/tau1;
kd2=s2*J2/n02+1/tau2+sigma2*f2;
kr2=s2*J2/n02+1/tau2;
kc1=s1*J1/n02;
kc2=s2*J2/n01;

nr1=(s1*J1*kr2-s2*J2*kc1)/(kr1*kr2-kc1*kc2);
nr2=(s2*J2*kr1-s1*J1*kc2)/(kr1*kr2-kc1*kc2);
nd1=(s1*J1*kd2-s2*J2*kc1)/(kd1*kd2-kc1*kc2);
nd2=(s2*J2*kd1-s1*J1*kc2)/(kd1*kd2-kc1*kc2);

Dn1=nr1-nd1;
Dn2=nr2-nd2;
kappa=sqrt((kd1-kd2)^2/4+kc1*kc2);
kda=(kd1+kd2)/2;
N1=((kda-kd1)*Dn1-kc1*Dn2)/kappa;
N2=((kda-kd2)*Dn2-kc2*Dn1)/kappa;


for i=1:max(size(td))
navg1(i)=quad(@integrnd, 0, td(i))/td(i);
navg2(i)=quad(@integrnd2, 0, td(i))/td(i);
%ntilde1(i)=navg1(i)/n01;
%ntilde2(i)=navg2(i)/n02;
R1(i)=V1*sigma1*f1*navg1(i);
R2(i)=V2*sigma2*f2*navg2(i);
R(i)=R1(i)+R2(i);
end

function n1=integrnd(ta)
%n1=X(1)*exp(r1*ta)+X(2)*exp(r2*ta)+nd1;
%n1=nd1+exp(-kda*ta).*(Dn1*cosh(kappa*ta)+N1*sinh(kappa*ta));
n1=nd1+.5*((Dn1+N1)*exp((kappa-kda)*ta)+(Dn1-N1)*exp((-kappa-kda)*ta));
end

function n2=integrnd2(ta)
%n2=X(3)*exp(r1*ta)+X(4)*exp(r2*ta)+nd2;
%n2=nd2+exp(-kda*ta).*(Dn2*cosh(kappa*ta)+N2*sinh(kappa*ta));
n2=nd2+.5*((Dn2+N2)*exp((kappa-kda)*ta)+(Dn2-N2)*exp((-kappa-kda)*ta));
end

end

