clear; close all

Nx = 2e2;     % amostragem em x
Ny = 1e2;     % amostragem em y

% Variavel independente no dominio espacial [a,b) x [c,d)

a = -4;
b = 6;
dx = (b-a)/Nx;       % incremento espacial em x
x = a + dx*(0:Nx-1);
c = -2;
d = 2;
dy = (d-c)/Ny;       % incremento espacial em y
y = c + dy*(0:Ny-1);

[xx,yy] = meshgrid(x,y); % grid espacial

% Funcao a ser transformada:
%
% f(x,y,z) = 1/sqrt(x^2+y^2+z^2)
%
% neste caso, f_xx + f_yy + f_zz = 0 e 
%
% f_z(x,y,z) = -1./sqrt(xx.^2+yy.^2+1).^3
% (o sinal de menos e' retirado, para considerar o crescimento 
% de cima p/ baixo na direcao z - Blakely, p. 326)

xo = (a + b)/2;
yo = (c + d)/2;
ff = 1./sqrt(xx.^2+yy.^2+1);

DF_2D(xx,yy,ff);  % compara derivadas espaciais

[dxf,dyf,dzf] = DF_2D(xx,yy,ff);

figure
subplot(1,2,1)
     mesh(xx,yy,1./sqrt(xx.^2+yy.^2+1).^3)
     title('dzf exata')
     view(30,30)
     za = zlim;
subplot(1,2,2)
     mesh(xx,yy,dzf)
     title('dzf aproximada')
     zlim(za);
     view(30,30)
