M-FILES USED FOR EN#214 ARTICLE ON SCALES
 
CODE IS NOT OPTIMIZED OR REFINED - JUST THE CODE ACTUALLY USED



PROGRAM 1 - THIS PROGRAM SCALE1 MAKES THE TABLES ON PAGE 5
*************BEGIN PROGRAM 1*******************************************************************
% scale1.m

tt=1;
ttt=1;
tttt=1;
ttttt=1
for k=0:12
   t5(k+1)=tt;
   t4(k+1)=tttt;
   t52(k+1)=ttt;
   t42(k+1)=ttttt;
   tt=1.5*tt;
   ttt=1.5*ttt;
   if ttt>2;ttt=ttt/2;end
   tttt=(4/3)*tttt;
   ttttt=(4/3)*ttttt;
   if ttttt>2;ttttt=ttttt/2;end
   
end
t5
t4
t52
t52(13)=2*t52(13);
t52s=sort(t52) 
t42s=sort(t42)
teqtemp=2.^([0:12]/12)
%tdisp=[0:12;t5;t52;t52s;t4;t42;t42s;teqtemp]'
tdisp5=[0:12;t5;t52;t52s;teqtemp]'
tdisp4=[0:12;t4;t42;t42s;teqtemp]'
 
**************END PROGRAM 1********************************************************************




THIS PROGRAM DOES THE MAIN CALCULATIONS/PLOTS OF THE ERROR FITTINGS
**************************BEGIN PROGRAM 2***********************************************************
% scale3.m
sqerr=[]
maxerr=[]

for N=5:50
   
n=[0:N];
r=2^(1/N);
tn=r.^n;

favrat=5/4;
me1=100;
for m=1:N
   e=abs(tn(m)-favrat);
   if e<me1;me1=e;end
end
favrat=4/3;
me2=100;
for m=1:N
   e=abs(tn(m)-favrat);
   if e<me2;me2=e;end
end
favrat=3/2;
me3=100;
for m=1:N
   e=abs(tn(m)-favrat);
   if e<me3;me3=e;end
end
favrat=5/3;
me4=100;
for m=1:N
   e=abs(tn(m)-favrat);
   if e<me4;me4=e;end
end

sqerr(N)=me1^2+me2^2+me3^2+me4^2;
maxerr(N)=max([me1 me2 me3 me4]);
toterr(N)=me1+me2+me3+me4;
% record error magnitudes
me1r(N)=me1;
me2r(N)=me2;
me3r(N)=me3;
me4r(N)=me4;


end

figure(1)
stem([0:50],[0 sqerr])
hold on
plot([-5 55],[0 0],'k')
hold off
axis([-2 55 -.002 0.022])
title('Total Square Error')

figure(2)
stem([0:50],[0 maxerr])
hold on
plot([-5 55],[0 0],'k')
hold off
axis([-2 55 -.01 0.1])
title('Max Abs Error')

figure(3)
stem([0:50],[0 toterr])
hold on
plot([-5 55],[0 0],'k')
hold off
axis([-2 55 -.01 0.28])
title('Total Abs Error')

figure(4)
stem([0:50],[0 me1r])
hold on
plot([-5 55],[0 0],'k')
hold off
axis([-2 55 -.005 0.075])
title('Abs Error at Third')

figure(5)
stem([0:50],[0 me2r])
hold on
plot([-5 55],[0 0],'k')
hold off
axis([-2 55 -.005 0.080])
title('Abs Error at Fourth')

figure(6)
stem([0:50],[0 me3r])
hold on
plot([-5 55],[0 0],'k')
hold off
axis([-2 55 -.005 0.090])
title('Abs Error at Fifth')

figure(7)
stem([0:50],[0 me4r])
hold on
plot([-5 55],[0 0],'k')
hold off
axis([-2 55 -.005 0.085])
title('Abs Error at Sixth')

***************************END PROGRAM 2***************************************************



PROGRAM 3 - THIS PROGRAM PRODUCES FIG. 4
********************BEGIN PROGRAM #******************************************
% scale4.m
clear
ra=1
for k=2:6
   for m=2:6
      rat(ra)=k/m;
      ra=ra+1;
   end
end
figure(1)
for k=1:length(rat)
   if rat(k)>0.99999
      if rat(k)<2.00001
         plot(rat(k),0.6,'r*')
         hold on
      end
   end
end

ra=1
for k=2:12
   for m=2:12
      rat(ra)=k/m;
      ra=ra+1;
   end
end
figure(1)
for k=1:length(rat)
   if rat(k)>0.99999
      if rat(k)<2.00001
         plot(rat(k),0.4,'g*')
         hold on
      end
   end
end

ra=1
for k=2:24
   for m=2:24
      rat(ra)=k/m;
      ra=ra+1;
   end
end
figure(1)
for k=1:length(rat)
   if rat(k)>0.99999
      if rat(k)<2.00001
         plot(rat(k),0.2,'b*')
         hold on
      end
   end
end

ra=1
for k=2:48
   for m=2:48
      rat(ra)=k/m;
      ra=ra+1;
   end
end
figure(1)
for k=1:length(rat)
   if rat(k)>0.99999
      if rat(k)<2.00001
         plot(rat(k),0,'k*')
         hold on
      end
   end
end

hold off
axis([0.9 2.1 -.3 .7])
figure(1)
**************************END PROGRAM 3*****************************************



PROGRAM 4 - THIS PROGRAM COMPARES FITS TO HARMONICS (Fig 10-12, etc)
************************BEGIN PROGRAM 4*****************************************
% scale11.m

clear
n=0:12;
s=2.^(n/12);
 
 
for n=1:32
   d=16;
    
   if n<16.0001;  d=8;  end
   if n< 8.0001;  d=4;  end
   if n< 4.0001;  d=2;  end
   if n< 2.0001;  d=1;  end
    
   hb(n)=n/d;
end
hb
hbs=sort(hb)
s

figure(1)
plot([0 3],[0 0],'k')
hold on
plot(s,0,'r*')
stem(hb,ones(1,32))
hold off
axis([0.9 2.1 -.2 1.2])
figure(1)
************************END PROGRAM 4 **************************************************




PROGRAM 5 - THIS PROGRAM PRODUCES FIG A1
************************BEGIN PROGRAM 5*************************************************
% scale13.m  to overplot Petersen/Barney Chart
clf
clear
FC=[1:1400];
FD=(9/8)*FC;
FE=(5/4)*FC;
FF=(4/3)*FC;
FG=(3/2)*FC;
FA=(5/3)*FC;
FB=(15/8)*FC;
FCC=2*FC;
FDD=2*FD;
FEE=2*FE;
FFF=2*FF;
FGG=2*FG;
FAA=2*FA;
FBB=2*FB;
FCCC=2*FCC;
FDDD=2*FDD;
FEEE=2*FEE;
FFFF=2*FFF;
FGGG=2*FGG;
FAAA=2*FAA;
FBBB=2*FBB;
FCCCC=2*FCCC;

figure(1)
plot(FC,log10(FC),'*c')
hold on
plot(FC,log10(FD),'*c')
plot(FC,log10(FE),'*c')
plot(FC,log10(FF),'*c')
plot(FC,log10(FG),'*c')
plot(FC,log10(FA),'*c')
plot(FC,log10(FB),'*c')
plot(FC,log10(FCC),'*c')
plot(FC,log10(FDD),'*c')
plot(FC,log10(FEE),'*c')
plot(FC,log10(FFF),'*c')
plot(FC,log10(FGG),'*c')
plot(FC,log10(FAA),'*c')
plot(FC,log10(FBB),'*c')
plot(FC,log10(FCCC),'*c')
plot(FC,log10(FDDD),'*c')
plot(FC,log10(FEEE),'*c')
plot(FC,log10(FFFF),'*c')
plot(FC,log10(FGGG),'*c')
plot(FC,log10(FAAA),'*c')
plot(FC,log10(FBBB),'*c')
plot(FC,log10(FCCCC),'*c')

axis([1 1400 log10(500) log10(4000)])
hold off
figure(1)
*****************************END PROGRAM 5******************************