close all; clear; clc;

%%
c = @(x) cos(2*pi*x); s = @(x) sin(2*pi*x);
f = @(x,y) c(x).*s(y); % Function

% Analytical Jacobian and Hessian
J = @(x,y) -2*pi*[s(x).*s(y) -c(x).*c(y)];
H = @(x,y) -4*pi^2*... 
    [c(x).*s(y), c(y).*s(x);
     c(y).*s(x), c(x).*s(y)];

%% Plot
cm = viridis;
N = 50; x = linspace(-1/2,1/2,N);
[xx,yy] = meshgrid(x);

% Path
r = 1/4;
ps = exp(1i*linspace(0,2*pi,50+1))*r; ps(end) = [];
ps = [[1:-0.2:-1 -1:0.2:1]*r ps];

for j = 1:numel(ps) % For each point of the path
    x0 = [real(ps(j)) imag(ps(j))].';
    
    figure(1); clf;
    surf(xx,yy,f(xx,yy)); hold on
    colormap(cm);
    Q = zeros(size(xx));
    for xi = 1:N
        for yi = 1:N
            xp = [x(xi); x(yi)];
            % Taylor expansion of f
            Q(yi,xi) = f(x0(1),x0(2)) + J(x0(1),x0(2))*(xp-x0) + ...
                1/2*(xp-x0)'*H(x0(1),x0(2))*(xp-x0);
        end
    end
    surf(xx,yy,Q,'edgecolor','none','facealpha',1/2,'facecolor',cm(end,:))
    plot3(x0(1),x0(2),f(x0(1),x0(2)),'ro');
    zlim([-1 1]); caxis([-1 1]);
    view(-20,64);
    drawnow;
end
