clear all
D=20; % number of subinterval for alpha in [ainit,aend]
ainit=1.2;aend=3.5;ainc=(aend-ainit)/D;
alpha=ainit:ainc:aend;
tau_1=1;
R=[];
N=10; % largest number of bases
figure(201);
hold off
for n=2:2:N % even base number only
    RR=[];
    HO=ceil((1:n)/2); % H-O bond strength with H-O = 1 normalized
    E=2+HO; % hydrogen bond strength including H-N bond with 2 = H-N bond strength. 
    P=E/sum(E); % base distribtuion
    for i=1:(D+1) % To generate the rate as a function of alpha=tau_{3,4}/tau_{1,2}. 
        H=-P*(log(P)/log(2))';
        tau=1+(alpha(i)-1)*(HO-1); % with tau_1=1 normalized, HO happsens to be the sequence 1 1 2 2 3 3, etc
        T=tau*P';
        RR=[RR,H/T]; % replication rate as function of alpha
    end
    R(n,:)=RR; 
end
plot(alpha,R(2,:),'r-',alpha,R(4,:),'b-',alpha,R(6,:),'y-',alpha,R(8,:),'m-',alpha,R(10,:),'g-')
legend('\itn = 2','\itn = 4','\itn = 6','\itn = 8','\itn = 10')
xlabel('\it\alpha = \tau_{3,4} / \tau_{1,2}')
ylabel('Normalized Rate \it\tau_{1,2} R')


    