function[ddx, ddy, ddx2, ddy2, lap,av]=lapoprectdomain(X,Y) ;
%define linear operator and matrix for the grided system
%where x and y are the input domain defined from meshgrid

delx=X(1,2)-X(1,1); dely=Y(2,1)-Y(1,1); %spacing in x and y
c=size(X); NN=c(1,1); N=c(1,2);

%ddx2=trid(1/(delx^2), -2/(delx^2), 1/(delx^2), N);

e=ones(N,1);
ddx2=spdiags((1/(delx^2))*[1*e -2*e 1*e], [-1 0 1], N, N);
ddx2(N,1)=1/(delx^2); ddx2(1,N)=1/(delx^2);
%operator such that second x  derivative = psi*ddx2 w where psi is a NN by N grid  

%and congruently the first derivative
%ddx=trid(-1/(2*delx), 0, 1/(2*delx), N);ddx(N,1)=1/(2*delx); ddx(1,N)=-1/(2*delx);

ddx=spdiags((1/(2*delx))*[-1*e 1*e], [-1  1], N, N);ddx(N,1)=1/(2*delx); ddx(1,N)=-1/(2*delx);
avx=spdiags([1*e 2*e 1*e], [-1 0 1], N, N);

%avx=trid(1,2,1, N); %averaging algorithm
%operator such that the x derivative = psi*ddx w where psi is a NN by N grid  

ddx2=DIAGo(ddx2, NN, N);ddx=DIAGo(ddx, NN, N); avx=DIAGo(avx, NN,N);

%ddx2=full(ddx2);

%operator such that second x derivative is ddx2 * psi where psi is column
%vector labeled across rows (and same for ddx

%PSI0=reshape(psi0', NN*N,1);

%ddx2(1,1)=0; ddx2(NN*N, NN*N)=0; 
%to make a non singular operator matrix
%converts from a 2d gridded matrix to a column where entries adjacent in x
%are adjacent in the colum vecotr

XX=reshape(X',N*NN,1); YY=reshape(Y', N*NN,1);
%puts mesh into column vector form

%now repeat with y derivative

%ddy2=trid(1/(dely^2), -2/(dely^2), 1/(dely^2),NN); ddy2(NN,1)=1/(dely^2); ddy2(1,NN)=1/(dely^2);
%%% will have to put in the right numbers into this once we get those

e=ones(NN,1);
ddy2=spdiags((1/(dely^2))*[1*e -2*e 1*e], [-1 0 1], NN, NN);
ddy2(1,1:4)=[2/(dely^2) -5/(dely^2) 4/(dely^2) -1/(dely^2)]; %implement the forward and backwards finite differencing
ddy2(NN,NN-3:1:NN)=[-1/(dely^2) 4/(dely^2) -5/(dely^2) 2/(dely^2)]; %at the boundaries


ddy22=spdiags((1/(dely^2))*[1*e -2*e 1*e], [-1 0 1], NN, NN);
ddy22(1,1:4)=[0/(dely^2) -5/(dely^2) 4/(dely^2) -1/(dely^2)]; %implement the forward and backwards finite differencing
ddy22(NN,NN-3:1:NN)=[-1/(dely^2) 4/(dely^2) -5/(dely^2) 0/(dely^2)]; %at the boundaries
%ddy2=full(ddy2)

%ddy2(1, 1:2)=[-1/(dely^2) 1/(dely^2)]; ddy2(NN, NN-1:NN)=[-1/(dely^2) 1/(dely^2)];
%for the single periodic domain

%operator such that second y derivative is ddy2 * psi

%and congruently for the first derivative
%ddy=trid(-1/(2*dely), 0, 1/(dely*2),NN);

ddy=spdiags((1/(2*dely))*[-1*e 1*e], [-1  1], NN, NN);
ddy(1,1:3)=[-3/(dely*2) 4/(dely*2) -1/(dely*2)]; %implement the forward finite difference at boundary
ddy(NN,NN-2:1:NN)=[1/(2*dely) -4/(2*dely) 3/(2*dely)]; %and the backward
avy=trid(1,2,1,NN);

%operator such that second y derivative is ddy2 * psig

%diagnolize as if we reversed x and y in our numbering system
DDY2=DIAGo(ddy2, N, NN); DDY=DIAGo(ddy, N, NN); AVY=DIAGo(avy, N,NN); DDY22=DIAGo(ddy22, N, NN);
%operator such that second y derivative would be DDY2 * psi IF PSI WAS
%COLUMN VECTOR OF REVERSED(REVERSED) numbering system (number down columns)

for n=1:N
    for nn=1:NN
        ddy2(n+N*(nn-1),1:N*NN)=reshape(reshape((DDY2(nn+NN*(n-1),:)),NN,N)',1,N*NN);
        ddy(n+N*(nn-1),1:N*NN)=reshape(reshape((DDY(nn+NN*(n-1),:)),NN,N)',1,N*NN);
        avy(n+N*(nn-1),1:N*NN)=reshape(reshape((AVY(nn+NN*(n-1),:)),NN,N)',1,N*NN);
        ddy22(n+N*(nn-1),1:N*NN)=reshape(reshape((DDY22(nn+NN*(n-1),:)),NN,N)',1,N*NN);
    end
end
%double reshape reverses the numbering system

%VORT00=reshape(vort00,N,NN)' %useful for converting back to grid for
%countour

ddx=sparse(ddx); ddy=sparse(ddy); ddx2=sparse(ddx2); 
X=reshape(XX,N,NN)'; Y=reshape(YY,N,NN)';

lap=ddx2+ddy22; lap=full(lap);
% for n=1:1:1
%      lap(n,n)=0;
% end
% lap(NN*N-N+1,NN*N-N+1)=0; 
av=(avy+avx)./8;
av=sparse(av);
