function [dIydr,dIydz] = meqIyfd(Iy,dr,dz)
% [dIydr,dIydz] = meqIyfd(Iy,ry,zy)
% Calculate the spatial derivative in the r and z direction of Iy using
% second order finite difference and assuming the value of Iy to be 0
% outside the grid
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

  %grid size
  [nz,nr] = size(Iy);
  
  dIydr = zeros(nz,nr);
  dIydz = zeros(nz,nr);
  
  %derivatives wrt r
  dIydr(:,1)      =   Iy(:,2)                     ./ ( 2*dr );
  dIydr(:,nr)     =  -Iy(:,nr-1)                  ./ ( 2*dr );
  dIydr(:,2:nr-1) = ( Iy(:,3:nr) - Iy(:,1:nr-2) ) ./ ( 2*dr );
  
  %derivatives wrt z
  dIydz(1,:)      =   Iy(2,:)                     ./ ( 2*dz );
  dIydz(nz,:)     =  -Iy(nz-1,:)                  ./ ( 2*dz );
  dIydz(2:nz-1,:) = ( Iy(3:nz,:) - Iy(1:nz-2,:) ) ./ ( 2*dz );
end