%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%
%%%%% Even-Based Channel Capacity Rate --- Linear Time Scales:
%%%%%     tau_2n-1=tau_2n=tau_1(1+(alpha-1)(n-1))
%%%%%     p_(2k-1)=p_(2k)=p^(1+(alpha-1)(k-1)); k=1 to n
%%%%%     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
alpha=2; % 
p_A=[];
p_C=[];
CRate=[];
Rmd=[];
solup=[];
alphaa=1.5;
alphab=3;
alphainc=.25;
%for alpha=1.65:.5:3
for alp=alphaa:alphainc:alphab
      p=.1:.00001:.9;
      n=2;   
      y=p+p.^alp-0.5;
      %y=2*p.^((alp-1)*n+1)-p.^(alp-1)-2*p+1;
      [rmd,ind]=min(abs(y));
      Rmd=[Rmd rmd];
      p_A=[p_A p(ind)];
      p_C=[p_C p(ind)^alp];
end
Rmd 
CRate=-log(p_A);
figure(1)
hold off
alpha=alphaa:alphainc:alphab;
plot(alpha,2*p_A,'k-', 'linewidth', 3)
hold on
plot(alpha,2*p_C,'b-', 'linewidth', 1)
%plot(alpha, CRate, 'k-', 'linewidth', 3)
%plot([alphaa alphab],[.54 .54],'b-.',[alphaa alphab],[.54-.09 .54-.09],'b-.',[alphaa alphab],[.54+.09 .54+.09],'b-.',...
%    [alphaa alphab],[.38 .38],'g-.',[alphaa alphab],[.38-.09 .38-.09],'g-.',[alphaa alphab],[.38+.09 .38+.09],'g-.')
%plot([alphaa alphab],[.54+0.09*54/(54+38) .54+0.09*54/(54+38)],'g-',[alphaa alphab],[.38+0.09*38/(54+38) .38+0.09*38/(54+38)],'g-')
plot([alphaa alphab],[.54+0.09*54/(54+38) .54+0.09*54/(54+38)],'k--','linewidth',3)
plot([alphaa alphab],[.38+0.09*38/(54+38) .38+0.09*38/(54+38)],'b--','linewidth',1)

legend('Capacity p_A + p_T','Capacity p_C + p_G','Human p_A + p_T','Humam p_C + p_G','fontsize',16)
%legend('Capacity p_A + p_T','Capacity p_C + p_G','\tau_A K','Human p_A + p_T','Humam p_C + p_G','fontsize',16)
plot([5/3 5/3],[0 1.4],'g--')
axis([alphaa alphab 0 1.2])
%title('Replication Time Scales:  \tau_{2k-1} = \tau_{2k} = \tau_1 + \Delta\tau(k-1)','fontsize',18)
xlabel('\alpha = \tau_{3,4} / \tau_{1,2}','fontsize',18)
ylabel('Distributions','fontsize',16)

%plot([5/3 5/3],[.54+0.09 .54+0.09],'mo',[5/3 5/3],[.38 .38],'go','markersize',25)

