Simulation digital pll

[Durak]

Hello. It is impossible to calculate PLL loop filter. I suffer more than one day. Calculated using the formula:


result simulation in matlab for sampling rate Fs = 4E+3:

For Fs = 10E+3:


It works only when Fs = 4E + 3, with which it can be connected? I can not understand where the error is.
As the synchronization (time frequency lock) does not match a calculated. How to calculate the loop filter?

Code in matlab:

Kd=5.5;
Ko=1;

Fs = 4E+3;
Ts = 1/Fs;

demp = 1/sqrt(2);
dF0max = 1000;
Tf=100e-6;

B = (4*dF0max^2 / Tf)^(1/3);
Wn = 2*B/(demp+1/(4*demp));


a0 = 4*(Wn*Ts)^2 +demp*Wn*Ts+4;
a1 = 2*(Wn*Ts)^2 - 8;
a2 = (Wn*Ts)^2 - 4*demp*Wn*Ts+4;

g0 = 1/(Kd*Ko)*(a1/a0+2);
g1 = 1/(Kd*Ko)*(a2/a0-1);

N = 100;
F0 = 1E+3;
pho = 90/180*pi;
s = sin(2 * pi * F0 * Ts * (1:N) + pho);
Fnco = F0 - 800;

e(1) = 0;
phi_output(1) = 0;
nco(1) = 0;
pnco(1) = 0;

for k = 2:N

nco(k) = sin(2*pi*Fnco*Ts*(k-1) + pnco(k-1));

e(k) = Kd * s(k-1) * nco(k);

phi_output(k) = g0*e(k) + g1*e(k-1) + phi_output(k-1);

pnco(k) = Ko * phi_output(k) + pnco(k-1);

end;

n = N;

plot((1:n)*Ts, s(1:n),(1:n)*Ts, nco(1:n));
figure(2);
plot((1:n)*Ts, phi_output(1:n));