y1=0
x1=0
x2=0
z1=1
y=1
T=10u
pi=3.14159
f=40k
w=2*pi*f
s=sin(w*T*180/pi)
c=cos(w*T*180/pi)
gain = -2*pi*f*T
cor=tan(f*pi*T*180/pi)/(pi*T*f)
n=0
len=100
xx=vector(len)
yy=vector(len)
time=vector(len)
while n < len
time[n]=n*T
x=y1*s+2*c*x1-x2
xx[n]=x
z=cor*gain*(x+x1)/2+z1
yy[n]=z
n=n+1
x2=x1
x1=x
y1=y
z1=z
y=0
end
plot xx time
plot yy
gain
cor