%take derivatives of a legendre polynomial on the surface semi anayltically
%output is derivative with respect to y
function out=divleg(m,n,lat)  


%latitude input in radians
%n is the nth polynomial for a given n
%the traditional n is then n+m-1
a=6.37*10^6;
NN=length(lat);
z=sin(lat);
y=a*lat;
N=n+m-1;
P=legendre(N,z);
Pn=P(m+1, :);
if n>1; P=legendre(N-1, z);
    Pnminus=P(m+1,:);
else Pnminus=zeros(1,NN);
end
out=((1./(z.^2-1))').*(N*z'.*Pn-(m+N)*Pnminus).*cos(lat')/a;
%co plot the anayltic result with the numerical
%der=zeros(1,NN);
%for n=2:1:NN-1
%    der(1,n)=(Pn(1,n+1)-Pn(1,n-1))/(y(1,n+1)-y(1,n-1));
%end
%der(1,1)=(Pn(2)-Pn(1))/(y(2)-y(1));
%der(1,NN)=(Pn(NN)-Pn(NN-1))/(y(NN)-y(NN-1));
%figure(1); clf;
%plot(z, 6000000*der, 'g', z, 6000000*out, 'b', z, Pn, 'r')
%legend('numerical', 'anayltic', 'function')
%title('input legendre polynomial and its derivatives')