function [xx,yy, field_out,spec] = spharm_trans_filt(lon, lat,field, trunc, wn);

% field=mSLP;
% wn=1;
% trunc=20
%routine to put the data into a gaussian quadrature grid and to preform the
%transform to the spherical harmonics and display the total wavenumbers
%specified by wn

[NN,N]=size(field);


%if nargin<4; error('Error - Not enough input parameters specified'), end;
%if nargin>4; error('Error - to0 many input parameters specified'), 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; lat=lat*pi/180; display('converting input latitudes to radians'); end;
if max(lon)>2*pi; lon=lon*pi/180; display('converting longitude input to 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);

%now interpolate the field onto the gaussian grid
ifield=interp2(LON, LAT, field, XX, YY);


N=length(xx); %new dimensions of the guassian quadrature field
NN=length(yy);

a=6.37*10^6; %radius of the Earth

num=trunc+1:-1:1; 
NNN=sum(num); %total number of basis functions
nor=zeros(NNN,1); %pre-allocate for the normalization
yy=asin(yy); %redefine the y coordinate
le=zeros(NN, trunc+1, trunc+1); %pre-allocate for legendre polynomials 
dle=zeros(NN, trunc+1, trunc+1); % pre - allocate for the anayltically defined y derivative of legendre polynomials

%put together the legendre polynomials and derivatives  as a function of n
%and m
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)=divleg(m,n+1-m,yy);
    end
end


%now put them in a column list- just pre sort some stuff here
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)/mean(abs(le(:,n, abs(m)+1)));%%%%%%%%%%%%%%%%%%%%%added 7/28/08 AD - to prevent huge numbers in basis
        %dLE(:,M+n)=dle(:,n,abs(m)+1);
    end
    M=M+n;
end

NNN=sum(num)+sum(num(2:length(num))); %would be 2*sum(num) but some are unneccesary
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);
    tn(n,1)=nn(n); %totalwavenumber
    
    basis(:,:,n+offset)=LE(:,n)*sin(mm(n)*xx);
    norm(n+offset,1)=sum(sum((basis(:,:,n).^2).*WW));
    nnn(n+offset,1)=-nn(n)*(nn(n)+1)/(a^2);
    tn(n+offset,1)=nn(n); %total wavenumber
end


for n=1:1:NNN
    spec(n,1)=sum(sum(basis(:,:,n).*WW.*ifield))/(norm(n));
end

if max(wn==0) == 1; %if the global mean is specified
    field_out=ones(NN,N)*mean(mean(WW.*ifield))/mean(mean(WW)); %start with the global mean
else
    field_out=zeros(NN,N);
end

for n=2:1:NNN
    if max(wn==tn(n)) == 1; %if the basis has a total wave# which is specfied 
        n;
        field_out=field_out+spec(n,1)*basis(:,:,n);
    end
end

figure(11); clf;
subplot(2,1,1)
load coast;
axesm('robinson','MapLatLimit',[-90 90],'MapLonLimit',[0 360])
load coast
plotm(lat,long, 'k');hold on;
hold on;
pcolorm(LAT*180/pi, LON*180/pi, field-mean(mean(cos(LAT).*field))/mean(mean(cos(LAT))));
color_trans; colormap(cmap);
MAX=max(max(abs(field-mean(mean(cos(LAT).*field))/mean(mean(cos(LAT))))));
caxis([-MAX MAX]);
colorbar;

figure(11); hold on;
subplot(2,1,2)
load coast;
axesm('robinson','MapLatLimit',[-90 90],'MapLonLimit',[0 360])
load coast
plotm(lat,long, 'k');hold on;
hold on;
color_trans; colormap(cmap);

%MAX=.4*max(max(abs(field_out)));
caxis([-MAX MAX]);
pcolorm(YY*180/pi,XX*180/pi, field_out);
colorbar;



%colorbar
%title('output (contours) and input (colors)')
