%% Optimize local absorbed power at the bottom (z = thickness)
% Append after the parameters in Cure Depth - Reflector.
% Uses existing I0 [mW/mm^2], thickness [mm], R, A_uvvis, L_cuvette [mm].
% Assumes fixed I0, thickness and R while alpha varies.
assert(isscalar(I0) && isfinite(I0) && I0 >= 0, 'I0 must be finite and nonnegative.');
assert(isscalar(thickness) && isfinite(thickness) && thickness > 0, 'Thickness must be finite and positive.');
assert(isscalar(R) && isfinite(R) && R >= 0 && R <= 1, 'R must lie between 0 and 1.');
assert(isscalar(A_uvvis) && isfinite(A_uvvis) && A_uvvis >= 0, 'Absorbance must be finite and nonnegative.');
assert(isscalar(L_cuvette) && isfinite(L_cuvette) && L_cuvette > 0, 'Cuvette length must be finite and positive.');

% P_bottom(alpha) = alpha * I0 * (1 + R) * exp(-alpha * thickness).
% Its derivative is proportional to exp(-alpha*d) * (1 - alpha*d).
% For I0 > 0, the maximum is at alpha = 1/d.
alpha_current = log(10) * A_uvvis / L_cuvette;
alpha_peak_location = 1 / thickness;
alpha_sweep = linspace(0, 6 * alpha_peak_location, 601);
P_bottom_no_mirror = I0 .* alpha_sweep .* exp(-alpha_sweep * thickness);
P_bottom_mirror = (1 + R) .* P_bottom_no_mirror;
P_bottom_current = I0 * (1 + R) * alpha_current * exp(-alpha_current * thickness);
P_bottom_max = I0 * (1 + R) / (exp(1) * thickness);

if I0 > 0
    alpha_opt = alpha_peak_location;
    A_uvvis_opt = alpha_opt * L_cuvette / log(10);
    fprintf('Optimal alpha: %.6g mm^-1\n', alpha_opt);
    fprintf('Equivalent UV-Vis absorbance: %.6g\n', A_uvvis_opt);
    fprintf('Maximum bottom power: %.6g mW/mm^3\n', P_bottom_max);
else
    alpha_opt = NaN;
    A_uvvis_opt = NaN;
    fprintf('I0 is zero: all alpha values give zero power; no unique optimum.\n');
end
fprintf('Current bottom power: %.6g mW/mm^3\n', P_bottom_current);

figure('Color', 'w', 'Position', [150, 150, 760, 480]);
plot(alpha_sweep, P_bottom_no_mirror, 'k:', 'LineWidth', 1.8, ...
    'DisplayName', 'No mirror'); hold on;
plot(alpha_sweep, P_bottom_mirror, 'Color', [0.47 0.25 0.29], ...
    'LineWidth', 2.5, 'DisplayName', sprintf('With mirror (R = %.2f)', R));
if I0 > 0
    plot(alpha_opt, P_bottom_max, 'd', 'Color', [0.33 0.45 0.39], ...
        'MarkerFaceColor', [0.33 0.45 0.39], 'MarkerSize', 9, ...
        'DisplayName', sprintf('Optimum: alpha = %.4g mm^{-1}', alpha_opt));
    xline(alpha_opt, '--', 'Color', [0.33 0.45 0.39], 'HandleVisibility', 'off');
end
if alpha_current <= alpha_sweep(end)
    plot(alpha_current, P_bottom_current, 'o', 'Color', [0.47 0.25 0.29], ...
        'MarkerFaceColor', 'w', 'MarkerSize', 8, 'DisplayName', 'Current setting');
else
    fprintf('Current alpha (%.6g mm^-1) is outside the plotted range [0, %.6g].\n', ...
        alpha_current, alpha_sweep(end));
end
grid on;
xlabel('Absorption coefficient \alpha (mm^{-1})');
ylabel('Local absorbed power at z = d (mW/mm^3)');
title(sprintf('Bottom absorbed power vs. alpha (d = %.4g mm)', thickness));
legend('Location', 'best');
xlim([0, alpha_sweep(end)]);
ylim([0, max(1.1 * P_bottom_max, eps)]);