function [flux,beta]=nma_getFlux1(A,k,dx,dy)
%flux function for PDE numerical course, UC Davis

[row,col]=size(A);
flux  = zeros(row,col);
beta  = zeros(row,col);
nx=row-2;
ny=col-2;

for i=2+nx-1:-1:2
    for j=2:1:2+ny-1
        qx=-k*(A(i,j+1)-A(i,j-1))/(2*dx);
        qy=-k*(A(i-1,j)-A(i+1,j))/(2*dy);
        flux(i,j)=sqrt(qx^2+qy^2);
        beta(i,j)=atan2(qy,qx)*180/pi;
    end
end   

end