function [xx,yy, WW,uout, vout, basis, vbasis, ubasis, norm,nnn,gdiv] = prefilter(lon, lat, u, v, trunc);

%prefilter routine to put the data into a gaussian quadrature grid and to
%serperate the velocity field into stream function and velocity potential

c=size(u);
N=c(1,2);
NN=c(1,1);


if nargin<5; error('Error - Not enough input parameters specified'), end;
if nargin>5; error('Error - to0 many input parameters specified'), end;
if size(u)~=size(v); error('Error- vector component matrices must be the same size'), end;
if length(lat)>NN;error('Error - vector of latitude values has too large of dimensions'),end;
if length(lat)<NN;error('Error - vector of latitude values has too small of dimensions'),end;
if length(lon)>N;error('Error - vector of longitude values has too large of dimensions'),end;
if length(lon)<NN;error('Error - vector of longitude values has too small of dimensions'),end;
if max(lat)>pi; error('Error- latitude input should be in radians'), end;
if max(lon)>2*pi; error('Error- longitude input should be in radians'), end; 


[LON LAT]=meshgrid(lon,lat);
GPy=round((3*trunc+1)/2+.999);   %number of Gaussian latitudes
[yy, wy]=gaussquad(GPy);          %sin(lat) and weights of Gaussian quadrature points 

%now do the same in x
GPx=2*GPy;   
xx=0:360/GPx:360-360/GPx;
xx=xx*pi/180;
xx=xx';

[XX YY]=meshgrid(xx,asin(yy));
[XX WW]=meshgrid(xx,wy);

vort=cur(u,v, lon,lat);
DIV=div(u,v,lon,lat);


%now interpolate the velocity field onto the gaussian grid
gvort=interp2(LON, LAT, vort, XX, YY);
gdiv=interp2(LON, LAT, DIV, XX, YY);

N=length(xx);
NN=length(yy);

a=6.37*10^6;

num=trunc+1:-1:1;
NNN=sum(num); 
M=0;
nor=zeros(NNN,1);
yy=asin(yy);
le=zeros(NN, trunc+1, trunc+1);
dle=zeros(NN, trunc+1, trunc+1);
for n=0:1:(trunc)
    P=legendre(n,sin(yy));
    for m=0:1:n;
        le(:,n+1-m, m+1)=(P(m+1,:))';
        nnn(n+1-m, m+1)=n;
        dle(:,n+1-m,m+1)=divy(le(:, n+1-m, m+1), yy);
    end
end

M=0;

xx=xx';
LE=zeros(NN, NNN);
dLE=zeros(NN, NNN);
for m=0:trunc
    for n=1:1:(1+trunc-abs(m))
        nor(M+n,1)=sum(sum(((le(:,n,abs(m)+1)*cos(m*xx)).^2).*WW));
        nn(M+n,1)=nnn(n,abs(m)+1);
        mm(M+n,1)=m;
        LE(:,M+n)=le(:,n, abs(m)+1);
        dLE(:,M+n)=dle(:,n,abs(m)+1);
    end
    M=M+n;
end

NNN=sum(num)+sum(num(2:length(num)));
basis=zeros(NN,N,NNN);
vbasis=zeros(NN,N,NNN);
ubasis=zeros(NN,N,NNN);
vvpbasis=zeros(NN,N,NNN);
uvpbasis=zeros(NN,N,NNN);

norm=zeros(NNN,1);
nnn=zeros(NNN,1);
offset=sum(num(2:length(num)));

for n=1:sum(num)
    basis(:,:,n)=LE(:,n)*cos(mm(n)*xx);
    norm(n,1)=sum(sum((basis(:,:,n).^2).*WW));
    nnn(n,1)=-nn(n)*(nn(n)+1)/(a^2);
    ubasis(:,:,n)=-dLE(:,n)*cos(mm(n)*xx);
    vvpbasis(:,:,n)=dLE(:,n)*cos(mm(n)*xx);
    for m=1:1:NN
        vbasis(m,:,n)=(-mm(n)/(a*cos(yy(m))))*LE(m,n)*sin(mm(n)*xx);
        uvpbasis(m,:,n)=(-mm(n)/(a*cos(yy(m))))*LE(m,n)*sin(mm(n)*xx);
    end
    
        
    basis(:,:,n+offset)=LE(:,n)*sin(mm(n)*xx);
    ubasis(:,:,n+offset)=-dLE(:,n)*sin(mm(n)*xx);
    vvpbasis(:,:,n+offset)=dLE(:,n)*sin(mm(n)*xx);
    for m=1:1:NN
        vbasis(m,:,n+offset)=(mm(n)/(a*cos(yy(m))))*LE(m,n)*cos(mm(n)*xx);
        uvpbasis(m,:,n+offset)=(mm(n)/(a*cos(yy(m))))*LE(m,n)*cos(mm(n)*xx);
    end
    
    norm(n+offset,1)=sum(sum((basis(:,:,n).^2).*WW));
    nnn(n+offset,1)=-nn(n)*(nn(n)+1)/(a^2);
end


%now invert the divergence and vorticity
psic=zeros(NNN,1);
chic=zeros(NNN,1);

for n=1:1:NNN
    psic(n,1)=sum(sum(basis(:,:,n).*WW.*gvort))/(norm(n)*nnn(n));
    chic(n,1)=sum(sum(basis(:,:,n).*WW.*gdiv))/(norm(n)*nnn(n));
end

psi=zeros(NN,N);
u=zeros(NN,N);
v=zeros(NN,N);
chi=zeros(NN,N);
uvp=zeros(NN,N);
vvp=zeros(NN,N);


for n=2:1:NNN
    psi=psi+psic(n,1)*basis(:,:,n);
    u=u+psic(n,1)*ubasis(:,:,n);
    v=v+psic(n,1)*vbasis(:,:,n);
    chi=chi+chic(n,1)*basis(:,:,n);
    uvp=uvp+chic(n,1)*uvpbasis(:,:,n);
    vvp=vvp+chic(n,1)*vvpbasis(:,:,n);
end

figure(1); clf;
axesm('robinson','MapLatLimit',[-90 90],'MapLonLimit',[0 360])
load coast
framem; patchm(lat,long, 'g')
contourm(YY*180/pi, XX*180/pi, chi, 15); 
hold on
quiverm(YY*180/pi,XX*180/pi, vvp,uvp)
title('MODERN velocity potential and associated velocities')


figure(7); clf;
axesm('robinson','MapLatLimit',[-90 90],'MapLonLimit',[0 360])
load coast
framem; patchm(lat,long, 'g')
contourm(YY*180/pi, XX*180/pi, psi, 15); 
hold on
quiverm(YY*180/pi,XX*180/pi, v,u)
title('LGM STREAMFUNCTION AND ASSOCIATED VELOCITIES')


uout=uvp+u;
vout=vvp+v;
