clc
close all
clear all
k = 45.7;
kp=0.5;
%T1 = 0.0001;
%T2 = 0.01;
tau = 6e-4;
%T22 = 1/tau.^2;

t = 0:tau:500;
 f = [0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 2 3 4 5 6 7 8 9 10 20 30 40 50 60 70 80 90 100 200 300 400 500 600 700 800 900 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 20000 30000 40000 50000];
x = zeros(length(t),length(f));
y = zeros(length(t),length(f));
 Ak = zeros(1,length(f));
for i = 1:length(f)
    x(1,i) = 0;
    y(1,i) = 0;
    x(2,i) = 0;
    y(2,i) = 0;
    for n = 3:(length(t)-1)
        x(n,i) = sin(f(i)*t(n));
        y(n,i) = kp*(x(n,i)-x(n-1,i))+y(n-1,i)+k*x(n,i)*tau;
    end
    Ak(i) = pi*sum(abs(y(:,i)-mean(y(:,i))))*0.5/length(y(:,i));
end





w = 0.1:0.1:50000;
A = sqrt(kp.^2.*w.^2+k.^2)./w;

    loglog(w,((A)));
grid on
hold on
   loglog(f,((Ak)),'*r');
% figure
%     plot(log10(w),fi);