clear all
close all

set(0, 'DefaultLineLineWidth', 1, ...
    'DefaultLineMarkerSize', 6,...
    'DefaultAxesFontSize', 18, ...
    'DefaultAxesTickLabelInterpreter', 'latex',...
    'DefaultFigurePaperOrientation', 'portrait',...
    'DefaultFigurePaperUnits', 'inches',...
    'DefaultFigurePaperSize', [6, 4.5],...
    'DefaultFigurePaperPosition', [0,0,6,4.5],...
    'DefaultLegendInterpreter', 'latex',...
    'DefaultTextInterpreter', 'latex',...
    'DefaultFigureRenderer', 'painters')

% Set up the grid
stepsize = 0.1
[x, y] = meshgrid(-3:stepsize:3, -3:stepsize:3);  % Tighter domain to match radius ~2.24
x = x';
y = y';

% Compute the objective function
z = (x.^2) .* exp(y) + 1.0 .* x .* (y.^2);

% Define custom contour levels
custom_levels = [-32 -16 -8 -4 -2 -1 1 4 8 16 25 48];

% Plot contours
close all;
figure;
set(gcf, 'Position', [10, 10, 800, 800]);

% ---------------- TOP VIEW ----------------
subplot(1,2,1)
colormap turbo
[C,h] =contour(x, y, z, custom_levels, 'ShowText', 'on', 'LineWidth', 3);
hold on;

% Plot constraint: x^2 + y^2 = 5
theta = linspace(0, 2*pi, 300);
r = sqrt(5);
xc = r * cos(theta);
yc = r * sin(theta);
zc = (xc.^2).* exp(yc) + 1.0*xc.*(yc.^2);
plot(xc, yc, 'k', 'LineWidth', 2, 'LineStyle','--');

xline(0, 'k');
yline(0, 'k');
xlabel('$x$')
ylabel('$y$')
grid
axis equal;
title('Top View', 'Interpreter','latex', 'FontSize', 20);

% Plot optimal points as stars
[~, maxIdx] = max(zc);
x_opt = xc(maxIdx);
y_opt = yc(maxIdx);
z_opt = zc(maxIdx);
plot(x_opt, y_opt, 'kp', 'MarkerSize', 12, 'MarkerFaceColor', 'k');  % Star

% ---------------- SIDE VIEW ----------------
subplot(1,2,2)
[x, y] = meshgrid(-2.5:0.25:2.5, -2.5:0.25:2.5);
F = (x.^2).* exp(y) + 1.0*x.*(y.^2);

% Plot the surface
surf(x, y, F, 'EdgeColor', 'interp', 'FaceAlpha', 0.9);
colormap turbo;
shading interp;
hold on;

contour3(x,y,F,custom_levels,'k', 'ShowText', 'off') % contour lines added

% Projected contour lines at z = -30
% contour(x, y, F, custom_levels);  % Get contour matrix

% Parse and plot the contour segments manually
z_proj = -50;  % Projection height
i = 1;
while i < size(C,2)
    level = C(1,i);
    numPoints = C(2,i);
    x_vals = C(1,i+1:i+numPoints);
    y_vals = C(2,i+1:i+numPoints);
    z_vals = z_proj * ones(1, numPoints);

    % Only plot if level is in custom_levels (for safety)
    if any(abs(level - custom_levels) < 1e-6)
        plot3(x_vals, y_vals, z_vals, 'r', 'LineWidth', 1);
    end
    i = i + numPoints + 1;
end


% Constraint curve
theta = linspace(0, 2*pi, 400);
r = sqrt(5);
xc = r * cos(theta);
yc = r * sin(theta);
zc = (xc.^2).* exp(yc) + 1.0*xc.*(yc.^2);

plot3(xc, yc, zc+.1, 'k', 'LineWidth', 2);       % Constraint on surface
plot3(xc, yc, z_proj*ones(1,400), 'k--', 'LineWidth', 1.5);  % Projection at z = -30

% Labels and formatting
xlabel('$x$'); ylabel('$y$'); zlabel('$F(x,y)$');
title('Side View', 'Interpreter','latex', 'FontSize', 20);
view(-15, 15);  % Adjust for optimal viewing angle

% Optimal point as star
plot3(x_opt,  y_opt, z_opt, 'kp', 'MarkerSize', 12, 'LineWidth', 2, 'MarkerFaceColor', 'k');
plot3(x_opt,  y_opt, z_proj,  'kp', 'MarkerSize', 12, 'LineWidth', 2, 'MarkerFaceColor', 'k');  % Projection



% Final figure formatting
x0 = 10; y0 = 10; width = 1000; height = 500;
set(gcf,'position',[x0,y0,width,height])
tightfig
saveas(gcf,'lect2_lagrange2c.pdf')