% J.C. Doyle (doyle@cds.caltech.edu)
% code for rank vs freq plot comparison

% define scaling rank function
Xmax=1e4;P=1:1e3;b=1.5;
[x,I,J]=unique(-floor(Xmax*P.^(-1/b))); 
x=-x;nx=length(x);P=P(I);
% create initial frequency data
p=[ 1 diff(P)];
% define exponential freq function
 a=2.13;x1=1000;
 pp=zeros(1,nx-1); ptest=1:200;
 pp(1)=P(1);ppinc=0;
% shift data to fit rank and freq functions
for k=2:nx-1
    ptest=p(k)+ppinc;
    xtest=x1*(1-log10(ptest)/a);
        if x(k)<xtest;
            pp(k)=0;
            ppinc=ppinc+p(k);
        else
            pp(k)=ptest;
            ppinc=0;
      end
end
pp(nx)=p(nx)+ppinc;
i=find(pp>0);pp=pp(i);x=x(i);
%Frequency plots
figure(4);hold off;semilogy(x,pp,'sb');hold on
axis([0 1200 0 2e2]);
figure(3);hold off; loglog(x,pp,'sb');hold on
axis([30 2e3 1e0 2e2]);
%rank plots
PP=cumsum(pp);
figure(1); hold off; loglog(x,PP,'sb'); 
axis([30 1e4 1e0 1e3]); hold on
figure(2); hold off;semilogy(x,PP,'sb'); 
axis([0 2e3 1e0 1e3]); hold on

% define exponential rank function
% then run the same algorithm
Xmax=1000;bb=3.1;b=bb/(Xmax);
P=1:1:1e3;
[x,I,J]=unique(-ceil(Xmax-log10(P)/b)); 
x=-x; nx=length(x);P=P(I);
p=[ 1 diff(P)];
% define scaling freq function
 a=1.8;x1=700;
 pp=zeros(1,nx-1);pp(1)=P(1);ppinc=0;
for k=2:nx-1
    ptest=p(k)+ppinc;
    xtest=x1*ptest^(-1/a);
        if x(k)<xtest;
            pp(k)=0;
            ppinc=ppinc+p(k);
        else
            pp(k)=ptest;
            ppinc=0;
      end
end
pp(nx)=p(nx)+ppinc;
i=find(pp>0); pp=pp(i); x=x(i);
% frequency plots
figure(4);semilogy(x,pp,'ok');hold on
figure(3);loglog(x,pp,'ok');hold on
%rank plots
PP=cumsum(pp);
figure(1); loglog(x,PP,'ok');
figure(2); semilogy(x,PP,'ok');

