% beginning of preprocessor
clear all 					% be safe 
fdel = input('Enter the size of the frequency step in Hertz > ');
fn = input('Enter the loop natural frequency in Hertz > ');
zeta = input('Enter zeta, the loop damping factor > ');
npts = 2000;           	% default number of simulation points
fs = 2000;					% default sampling frequency
T = 1/fs;
t = (0:(npts-1))/fs;    % time vector
nsettle = fix(npts/10); % set nsettle time as 0.1*npts
Kt=4*pi*zeta*fn;			% loop gain
a=pi*fn/zeta;				% loop filter parameter
filt_in_last = 0; filt_out_last = 0;
vco_in_last = 0; vco_out = 0; vco_out_last = 0;
% end of preprocessor
%
% beginning of simulation loop
for i=1:npts	
   if i<nsettle
      fin(i) = 0;
      phin = 0;
   else
      fin(i) = fdel;
      phin = 2*pi*fdel*T*(i-nsettle);
   end
        s1=phin - vco_out;
        s2=sin(s1);      % sinusoidal phase detector
        s3=Kt*s2; 
        filt_in = a*s3;
        filt_out = filt_out_last + (T/2)*(filt_in + filt_in_last);
        filt_in_last = filt_in;
        filt_out_last = filt_out;
        vco_in = s3 + filt_out;
        vco_out = vco_out_last + (T/2)*(vco_in + vco_in_last);
        vco_in_last = vco_in;
        vco_out_last = vco_out;
        phierror(i)=s1;
        fvco(i)=vco_in/(2*pi);
        freqerror(i) = fin(i)-fvco(i);
end
% end of simulation loop
%
% beginning of postprcessor
kk = 0;
while kk == 0
k = menu('Phase Lock Loop Postprocessor',...
        'Input Frequency and VCO Frequency',...
        'Phase Plane Plot',...
        'Exit Program');
     if k == 1
        plot(t,fin,t,fvco)
        title('Input Frequency and VCO Freqeuncy')
        xlabel('Time - Seconds')
        ylabel('Frequency - Hertz')
        pause
     elseif k == 2
        plot(phierror/2/pi,freqerror)
        title('Phase Plane')
        xlabel('Phase Error / pi')
        ylabel('Frequency Error - Hz')
        pause
     elseif k == 3
        kk = 1;
     end
end
% end of postprcssor