function [VV,V] =domainbasis(N,K)
%Gives orthonormal polynomials for a polygon K
%T=1000;
[T,~]=size(K);
T=T-1;
L=1/T;
[X,Y] = meshgrid(0:L:1, 0:L:1);

V=fliplr(eye(N));
VV=zeros(N);
for j=1:N
    for k=1:j
        u=V(j,:);
        v=V(k,:);
        w=polyval(u,X+1i*Y).*conj(polyval(v,X+1i*Y)).*K;
        a=real(w);b=imag(w);
        A=sum(sum(a))/(T+1).^2;B=sum(sum(b))/(T+1).^2;
        VV(j,k)=A+1i*B;
    end
    j/N
end
V(1,:)=V(1,:)/sqrt(VV(1,1));
for j=2:N
    D=zeros(1,N);
    for k=1:j-1
        D=D+VV(j,k)/VV(k,k)*V(k,:);
    end
    V(j,:)=V(j,:)-D;
end
for j=1:N
    u=V(j,:);
    w=polyval(u,X+1i*Y).*conj(polyval(u,X+1i*Y)).*K;
    a=real(w);b=imag(w);
    A=sum(sum(a))/(T+1).^2;B=sum(sum(b))/(T+1).^2;
    VV(j,j)=A+1i*B;
    V(j,:)=V(j,:)/sqrt(VV(j,j));
end
end