function [xx,yy,nnn, basis, WW, norm] = basis_setup2(trunc);

%routine to put the data into a gaussian quadrature grid and to preform the
%transform to the spherical harmonics
 


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);

WW=WW/mean(g2c(WW));


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;       
    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);

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);        
    
    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);
end

% basis(1:GPy,1:GPx,1)=zeros(GPy, GPx);
% 
% 
% for nn=2:NNN
%     basis(1:GPy, 1:GPx,nn)=basis(1:GPy,1:GPx,nn)/((mean((g2c(squeeze(basis(1:GPy, 1:GPx,nn))).^2).*g2c(WW)))^.5);
%     norm(nn,1)=mean(mean((basis(:,:,nn).^2).*WW));
% end
% norm(1,1)=1;


