function [ax,ay ]=gradient_perso(a,dx,dy,n);
%GRADIENT_PERSO - gradient computation
%
%   INPUT a,dx,dy,n where a is the scalar value, dx and dy the spacing of one pixel
%   in a and n is the width wanted in the gradient computation, in pixels.
%
%   Usage:
%      [ax,ay ]=gradient_perso(a,dx,dy,n)

s=size(a);
ax=zeros(s(1),s(2));
ay=zeros(s(1),s(2));

for k=n+1:s(1)-n,
   if mod(k,10)==0,
      disp(k/s(1)*100);
   end
   
   for j=n+1:s(2)-n,
      if isnan(a(k,j)),
         ax(k,j)=NaN;
         ay(k,j)=NaN;
      else
         temp=a(k,j);
         temp2=a(k,j);
         
         count=1;
         while ~isnan(a(k,j+count)),
            temp=[temp a(k,j+count)];
            count=count+1;
            if count>n,
               count=count-1;
               break;
            end
            
         end
         count=1;
         while ~isnan(a(k,j-count)),
            temp=[a(k,j-count) temp];
            count=count+1;
            if count>n,
               count=count-1;
               break;
            end
            
         end
         
         count=1;
         while ~isnan(a(k+count,j)),
            temp2=[temp2 a(k+count,j)];
            count=count+1;
               if count>n,
               count=count-1;
               break;
            end
         
         end
         count=1;
         while ~isnan(a(k-count,j)),
            temp2=[a(k-count,j) temp2];
            count=count+1;
               if count>n,
               count=count-1;
               break;
            end
         
         end
    	   
         if length(temp)==1,
            ax(k,j)=NaN;
         else
            ax(k,j)=(temp(length(temp))-temp(1))/(length(temp)-1)/dx;
         end
         
         if length(temp2)==1,
            ay(k,j)=NaN;
         else
            ay(k,j)=(temp2(length(temp2))-temp2(1))/(length(temp2)-1)/dy;
         end
         
      end
   end
end

         
          
          
          
          
          
      
