function out=sphdivx(varargin)

%ddx operator in spherical input
%input format I sphdivx(A,lat,lon) will set default order to 2
%input format II  sphydivx(A,lat,lon, 4) will make the differencing fourth order accurate  

%automatically detects periodic (full zonal domain)

if nargin<3; error('Error - Not enough input parameters specified'), end;
if nargin>4; error('Error - too many input parameters specified'), end;

A = varargin{1};
lat = varargin{2};
lon = varargin{3};

if nargin==3;
    ord=2; %default order   
else ord = varargin{4};
end

[NN N]=size(A);


ord;


if length(lon)>N;error('Error - vector of longitude values has too large of dimensions'),end;
if length(lon)<N;error('Error - vector of longitude values has too small of dimensions'),end;
if max(lon)>4; lon=lon*pi/180; end %convert to radians if input is in degree

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 max(abs(lat))>(pi/2); lat=lat*pi/180; end %convert to radians if input is in degree

out=zeros(NN,N);
a=6.37*10^6;


if ord==2
    if ((min(lon)-(max(lon)-2*pi)))<(1.1*(lon(2)-lon(1)))  %detects periodic domain
    %disp('DOMAIN IS PERIODIC- order of scheme is 2')
        
    for m=1:1:NN
        out(m,1)=(A(m,2)-A(m,N))/(a*cos(lat(m))*(lon(2)-(lon(N)-2*pi))); %second order centered differencing at periodic boundary
        out(m,N)=(A(m,1)-A(m,N-1))/(a*cos(lat(m))*(lon(1)-(lon(N-1)-2*pi)));
    end
    
    
    
    
    else
    %disp('DOMAIN IS NOT PERIODIC- order of scheme is 2')
    
    for m=1:1:NN
        out(m,1)= (-1*A(m,3)+4*A(m,2)-3*A(m,1))/(2*a*cos(lat(m))*(lon(2)-lon(1))); %second order accurate forward scheme
        out(m,N)= (3*A(m,N)-4*A(m,N-1)+1*A(m,N-2))/(2*a*cos(lat(m))*(lon(N)-lon(N-1)));
    end
    
    end %periodic if loop

    for n=2:1:N-1
        for m=1:1:NN
            out(m,n)=(A(m,n+1)-A(m,n-1))/(a*cos(lat(m))*(lon(n+1)-lon(n-1))); %second order centered differencing
        end
    end

elseif ord==4
    
    if ((min(lon)-(max(lon)-2*pi)))<(1.1*(lon(2)-lon(1)))  %detects periodic domain
    disp('DOMAIN IS PERIODIC- order of scheme is 4')
        
    for m=1:1:NN
        out(m,1)=(-A(m,3)+8*A(m,2)-8*A(m,N)+A(m,N-1))/(6*a*cos(lat(m))*(lon(2)-(lon(N)-2*pi))); %fourth order centered differencing at periodic boundaries
        out(m,2)=(-A(m,4)+8*A(m,3)-8*A(m,1)+A(m,N))/(6*a*cos(lat(m))*(lon(3)-(lon(1))));
        out(m,N)=(-A(m,2)+8*A(m,1)-8*A(m,N-1)+A(m,N-2))/(6*a*cos(lat(m))*(lon(1)-(lon(N-1)-2*pi)));
        out(m,N-1)=(-A(m,1)+8*A(m,N)-8*A(m,N-2)+A(m,N-3))/(6*a*cos(lat(m))*(lon(N)-(lon(N-2))));
    end
    
    else
        disp('DOMAIN IS NOT PERIODIC- order of scheme is 4 except at boundaries')
        
    for m=1:1:NN
        out(m,1)= (-1*A(m,3)+4*A(m,2)-3*A(m,1))/(2*a*cos(lat(m))*(lon(2)-lon(1))); %second order accurate forward scheme
        out(m,N)= (3*A(m,N)-4*A(m,N-1)+1*A(m,N-2))/(2*a*cos(lat(m))*(lon(N)-lon(N-1)));
    end
    
   
    for m=1:1:NN
        n=2;
        out(m,n)=(A(m,n+1)-A(m,n-1))/(a*cos(lat(m))*(lon(n+1)-lon(n-1))); %second order centered differencing
        n=N-1;
        out(m,n)=(A(m,n+1)-A(m,n-1))/(a*cos(lat(m))*(lon(n+1)-lon(n-1)));
    end
    
    end %periodic if loop

    for n=3:1:N-2
        for m=1:1:NN
            out(m,n)=(-A(m,n+2)+8*A(m,n+1)-8*A(m,n-1)+A(m,n-2))/(6*a*cos(lat(m))*(lon(n+1)-lon(n-1))); %fourth order centered differencing
        end
    end
else; error('input order must be 2 or 4');
end

       
    
%looks good 1/4/07