function out=sphdivy(varargin)

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

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 %display order of finite differencing


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
    for m=1:1:N
        out(1,m)= (-1*A(3,m)+4*A(2,m)-3*A(1,m))/(2*a*(lat(2)-lat(1))); %second order accurate forward scheme
        out(NN,m)= (3*A(NN,m)-4*A(NN-1,m)+1*A(NN-2,m))/(2*a*(lat(NN)-lat(NN-1)));        
    end

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

elseif ord==4
    for m=1:1:N
        out(1,m)= (-1*A(3,m)+4*A(2,m)-3*A(1,m))/(2*a*(lat(2)-lat(1))); %second order accurate forward scheme
        out(NN,m)= (3*A(NN,m)-4*A(NN-1,m)+1*A(NN-2,m))/(2*a*(lat(NN)-lat(NN-1)));
    end
    
    for m=1:1:N
        out(2,m)=(A(3,m)-A(1,m))/(a*(lat(3)-lat(1))); %second order center differencing for the point one in
        out(NN-1,m)=(A(NN,m)-A(NN-2))/(a*(lat(NN)-lat(NN-2)));
    end 

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

       
    
%looks good 1/4/07