function [Dxf,Dyf,Dzf] = DF_2D(xx,yy,ff)

% Derivada com a Transformada de Fourier via FFT
% Versao 2D, retorna as derivadas x e y
% Assume que (xx,yy) sao dadas por um grid de Ny x Nx pontos
%
% considerando o par (f(x),F(k)) definido (em cada direcao) por 
%
% F(k) = \int f(x)*exp(-i*2*pi*k*x) dx
% f(x) = \int F(k)*exp( i*2*pi*k*x) dk
%

I = sqrt(-1); % unidade imaginaria

[Ny,Nx] = size(xx);  % Extrai dimensoes do grid
x = xx(1,:);         % {x,y}: vetores que definem o grid
y = yy(:,1)';

% Adiciona um ponto artificial no final dos vetores {x,y} do grid

a  = x(1);
dx = x(2)-x(1);
b  = x(Nx)+dx;

c  = y(1);
dy = y(2)-y(1);
d  = y(Ny)+dy;

% Extensao do grid

x_esq = x - (b-a);
x_dir = x + (b-a);

y_esq = y - (d-c);
y_dir = y + (d-c);

[xx_ext,yy_ext]=meshgrid([x_esq,x,x_dir],[y_esq,y,y_dir]);

% Extensao da funcao f

ff_ext = zeros(size(xx_ext));

for i = 1:Ny % estensao em x, para cada valor de y

  f = ff(i,:);  % fatia correspondente a y=y(i)

  D = 10/(b-a);
  X = x_esq - x(1);
  A = f(1);
  B = (f(2)-f(1))/dx - A*D;
  C = (f(3)-2*f(2)+f(1))/dx^2 - 2*B*D - A*D^2;
  
  f_esq = (A + B*X + 0.5*C*X.^2).*exp(D*X);
  
  D = -D;
  X = x_dir - x(Nx);
  A = f(Nx);
  B = (f(Nx)-f(Nx-1))/dx - A*D;
  C = (f(Nx)-2*f(Nx-1)+f(Nx-2))/dx^2 - 2*B*D - A*D^2;

  f_dir = (A + B*X + 0.5*C*X.^2).*exp(D*X);

  ff_ext(Ny+i,:) = [f_esq,f,f_dir];
  
end

for j = 1:3*Nx % estensao em y, para cada valor de x (x ja' estendido)
  
  f = ff_ext(Ny+1:2*Ny,j);  % fatia correspondente a x=x(j)

  D = 10/(d-c);
  Y = y_esq - y(1);
  A = f(1);
  B = (f(2)-f(1))/dy - A*D;
  C = (f(3)-2*f(2)+f(1))/dy^2 - 2*B*D - A*D^2;
  
  f_esq = (A + B*Y + 0.5*C*Y.^2).*exp(D*Y);
  
  D = -D;
  Y = y_dir - y(Ny);
  A = f(Ny);
  B = (f(Ny)-f(Ny-1))/dy - A*D;
  C = (f(Ny)-2*f(Ny-1)+f(Ny-2))/dy^2 - 2*B*D - A*D^2;

  f_dir = (A + B*Y + 0.5*C*Y.^2).*exp(D*Y);

  ff_ext(:,j) = [f_esq';f;f_dir'];
  
end

% Grid do dominio da frequencia

A = 2*a-b;             % nova extremidade 'a esquerda
B = 2*b - a;           % nova extremidade 'a direita

k = (-3*Nx/2+1):(3*Nx/2); 
kappax = (k - 1)/(B-A);

C = 2*c-d;             % nova extremidade abaixo
D = 2*d - c;           % nova extremidade acima

k = (-3*Ny/2+1):(3*Ny/2); 
kappay = (k - 1)/(D-C);

[kpx,kpy] = meshgrid(kappax,kappay);

% Transformada via FFT

F = dx*dy*fft2(ff_ext);   % resultado no dominio da frequecia [0,N]^2
F = fftshift(F);          % redistribui o dominio da freq. para [-N/2,N/2]^2
F = exp(-I*2*pi*(A*kpx+C*kpy)).*F; % correcao das extremidades

% Calcula as derivadas parciais no dominio da frequencia

dxF = I*(2*pi)*kpx.*F;
dyF = I*(2*pi)*kpy.*F;
dzF = (2*pi)*sqrt(kpx.^2+kpy.^2).*F;

% Transformada inversa

dxF = exp(I*2*pi*(A*kpx+C*kpy)).*dxF; % desfaz a correcao das extremidades
dxF = ifftshift(dxF);                 % desfaz o shift
Dxf = (1/(dx*dy))*ifft2(dxF);         % aplica a FFT inversa
Dxf = real(Dxf);                      % remove a parte imaginaria

dyF = exp(I*2*pi*(A*kpx+C*kpy)).*dyF; % idem para dyF 
dyF = ifftshift(dyF);                 
Dyf = (1/(dx*dy))*ifft2(dyF);
Dyf = real(Dyf);

dzF = exp(I*2*pi*(A*kpx+C*kpy)).*dzF; % idem para dzF 
dzF = ifftshift(dzF);                 
Dzf = (1/(dx*dy))*ifft2(dzF);
Dzf = real(Dzf);

% Corta a resposta no intervalo [a,b]x[c,d]

indx = Nx+1:2*Nx; 
indy = Ny+1:2*Ny; 

xxc = xx_ext(indy,indx);
yyc = yy_ext(indy,indx);
Dxf = Dxf(indy,indx);
Dyf = Dyf(indy,indx);
Dzf = Dzf(indy,indx);

% Graficos (se nao houver argumentos de saida)

if(nargout==0)

  close all
  figure

  dxf = (ff(:,3:Nx)-ff(:,1:Nx-2))/(2*dx);  % derivadas calculadas 
  dxf(:,Nx) = dxf(:,Nx-2);                 % por difs finitas 
  dxf(:,2:Nx-1) = dxf(:,1:Nx-2);           % (para comparar)
  
  dyf = (ff(3:Ny,:)-ff(1:Ny-2,:))/(2*dy); 
  dyf(Ny,:) = dyf(Ny-2,:);
  dyf(2:Ny-1,:) = dyf(1:Ny-2,:);
 
  subplot(2,2,1)
     mesh(xx,yy,dxf)
     title('Derivada espacial x')
     view(30,30)
  subplot(2,2,2)
     mesh(xx,yy,dyf)
     title('Derivada espacial y')
     view(30,30)
  subplot(2,2,3)
     mesh(xxc,yyc,Dxf)
     title('Derivada via FFT x')
     view(30,30)
  subplot(2,2,4)
     mesh(xxc,yyc,Dyf)
     title('Derivada via FFT y')
     view(30,30)

  Dxf = []; Dyf = []; Dzf = []; 

end
