function out=sphlap(A,lat,lon)

%for non periodic domain

[NN N]=size(A);


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;

%do the y part first

[LON,LAT] = meshgrid(lon,lat);
Asinlat=sin(LAT).*sphdivy(A,lat,lon,2);
outy=(sphdivy(Asinlat,lat,lon,2))./sin(LAT);

%now do the x
if ((min(lon)-(max(lon)-2*pi)))<(1.1*(lon(2)-lon(1)))  %detects periodic domain

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

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

    
out=outx+outy;