clc;
clear;
close all;

disp('-----------------------------------------');
disp('NE565 PLL characteristic                 ');
disp('-----------------------------------------');

xk  = [3.83 4.01 4.33 4.62 4.81 5.05];
yk  = [4000 3500 3000 2500 2000 1500];
[p, s] = polyfit(xk, yk, 1);
x = [3.5:0.1:5.5];
y = polyval(p, x);

x2k = [40.3 57.9 79.8 99 149.8 138.4]*(pi/180);
y2k = [3.83 4.01 4.33 4.62 4.81 5.05];
[p2, s2] = polyfit(x2k, y2k, 1);
x2 = [0:1:180]*(pi/180);
y2 = polyval(p2, x2);

x3k = [1500 2000 2500 3000 3500 4000];
y3k = [138.4 149.8 99 79.8 57.9 40.3];
[p3, s3] = polyfit(x3k, y3k, 1);
x3 = [1000:1:4500];
y3 = polyval(p3, x3);

i90 = find( abs(y3 - 90) ==  min(abs(y3 - 90)) );

fo = 2900;
Vc = 12;

lgth = 33.6*fo/Vc;

Ko = 2*pi*abs(p(1));
KD = p2(1);

lgms = Ko*KD;

pe = 100*abs(lgth-lgms)/abs(lgth);

s = ['phase = 90 degrees at ',num2str(x3(i90)),' Hz'];
disp(s)
s = ['oscillator sensitivity = ',num2str(Ko),' rd/s per volt'];
disp(s)
s = ['phase detector sensitivity = ',num2str(KD),' volts per radian'];
disp(s)
s = ['loop gain theory = ',num2str(lgth)];
disp(s)
s = ['loop gain measured = ',num2str(lgms)];
disp(s)
s = ['percent error = ',num2str(pe),' %'];
disp(s)
disp(' ');

figure(1)
plot(x, y)
hold on
plot(xk, yk, 'rx')
s = ['VCO: fosc = ', num2str(p(1)),' Hz/volt x Vin + ',num2str(p(2)),' Hz'];
title(s)
xlabel('volts')
ylabel('Hz')
axis([3.5 5.5 1000 4500])

figure(2)
plot(x2, y2)
hold on
plot(x2k, y2k, 'rx')
s = ['phase detector: Vout = ', num2str(p2(1)),' volts/rad x Oin + ',num2str(p2(2)),' volts'];
title(s)
xlabel('radians')
ylabel('volts')

figure(3)
plot(x3, y3)
hold on
plot(x3k, y3k, 'rx')
s = ['y = ', num2str(p3(1)),' x + ',num2str(p3(2))];
title(s)
xlabel('Hz')
ylabel('degrees')