function [F,valmis] =integroi(f,z0,w0)
%Integroi kompleksisen funktion
%
%[F] =integroi(f,z0,w0)

[T,~]=size(f);
TT=T-1;
L=1/TT;

F=0*f;
FF=F;
FFF=F;
lkm=F;
valmis=F;

%W=round(TT*[real(z0) imag(z0)]);
W=z0;
a=w0;
F(W(1),W(2))=a;				%kiinnitetään arvo
valmis(W(1),W(2))=1;			%1 piste valmis

s=2/T;
j2=sqrt(2);
%for e=1:10
while any(valmis==0),
[n,m]=find(valmis==1);
N=length(n);

for k=1:N
	q=m(k);
	p=n(k)-1;
	if p~=0,
		if valmis(p,q)~=1,
			FF(p,q,k)=-f(p+1,q)*s+F(n(k),m(k));
			lkm(p,q)=lkm(p,q)+1;
		end
	end
	p=n(k)+1;
	if p~=T+1,
		if valmis(p,q)~=1
			FF(p,q,k)=f(p-1,q)*s+F(n(k),m(k));
			lkm(p,q)=lkm(p,q)+1;
		end
	end
	p=n(k);
	q=m(k)-1;
	if q~=0,
		if valmis(p,q)~=1,
			FF(p,q,k)=-1i*f(p,q+1)*s+F(n(k),m(k));
			lkm(p,q)=lkm(p,q)+1;
		end
	end
	q=m(k)+1;
	if q~=T+1,
		if valmis(p,q)~=1
			FF(p,q,k)=1i*f(p,q-1)*s+F(n(k),m(k));
			lkm(p,q)=lkm(p,q)+1;
		end
	end

	q=m(k)-1;
	p=n(k)-1;
	r=q*p;
	if r~=0,
		if valmis(p,q)~=1,
			FF(p,q,k)=(-1-1i)/j2*f(p+1,q+1)*s+F(n(k),m(k));
			lkm(p,q)=lkm(p,q)+1;
		end
	end
	p=n(k)+1;
	r=q*(p-T-1);
	if r~=0,
		if valmis(p,q)~=1
			FF(p,q,k)=(1-1i)/j2*f(p-1,q+1)*s+F(n(k),m(k));
			lkm(p,q)=lkm(p,q)+1;
		end
	end
	q=m(k)+1;
	r=(q-T-1)*(p-T-1);
	if r~=0,
		if valmis(p,q)~=1,
			FF(p,q,k)=(1+1i)/j2*f(p-1,q-1)*s+F(n(k),m(k));
			lkm(p,q)=lkm(p,q)+1;
		end
	end
	p=n(k)-1;
	r=(q-T-1)*p;
	if r~=0,
		if valmis(p,q)~=1
			FF(p,q,k)=(-1+1i)/j2*f(p+1,q-1)*s+F(n(k),m(k));
			lkm(p,q)=lkm(p,q)+1;
		end
    end
end

for j=1:T
for k=1:T
if lkm(j,k)~=0,
FFF(j,k)=sum(FF(j,k,:))/lkm(j,k);
valmis(j,k)=1;
end
end
end
F=F+FFF;
mean(mean(valmis))
end
end

