From eb0f6375c60b54fe5438687921232f8a91660171 Mon Sep 17 00:00:00 2001
From: Karthik Chandrashekara <karthik@physi.uni-heidelberg.de>
Date: Mon, 24 Mar 2025 15:08:16 +0100
Subject: [PATCH] Plot of the cost function for the mloop optimizer.

---
 Estimations/EstimateCostFunctionMLoop.m | 87 +++++++++++++++++++++++++
 1 file changed, 87 insertions(+)
 create mode 100644 Estimations/EstimateCostFunctionMLoop.m

diff --git a/Estimations/EstimateCostFunctionMLoop.m b/Estimations/EstimateCostFunctionMLoop.m
new file mode 100644
index 0000000..9d8747f
--- /dev/null
+++ b/Estimations/EstimateCostFunctionMLoop.m
@@ -0,0 +1,87 @@
+%%
+% Parameters
+alpha = 0.5;         % Example value for alpha
+N = linspace(1e4, 1e6, 100); % Atom number range
+OD = linspace(0, 3, 100);   % Optical depth range
+N1 = 100;            % Example parameter for N1 controlling transition steepness
+
+% Smoothed Heaviside function as given in the appendix
+SmoothedHeaviside = @(N) (N > 0) .* (2 ./ (1 + exp(N1 ./ N)));
+
+% Initialize the cost function
+CostFunction = zeros(length(N), length(OD));
+
+% Calculate the cost function for each N and OD
+for i = 1:length(N)
+    for j = 1:length(OD)
+        % Smoothed Heaviside step function
+        H_N = SmoothedHeaviside(N(i));
+        % Cost function: OD^3 * N^(alpha - 9/5) multiplied by smoothed Heaviside
+        CostFunction(i, j) = H_N * OD(j)^3 * N(i)^(alpha - 9/5);
+    end
+end
+
+% Create a meshgrid for plotting
+[OD_grid, N_grid] = meshgrid(OD, N);
+
+% Plot the cost function as a surface plot
+figure(1);
+clf
+set(gcf,'Position',[50 50 950 750])
+set(gca,'FontSize',16,'Box','On','Linewidth',2);
+surf(OD_grid, N_grid, CostFunction);
+xlabel('Optical Depth (OD)','FontSize',16);
+ylabel('Number of Atoms (N)','FontSize',16);
+zlabel('Cost Function','FontSize',16);
+title(['|C(X)| = OD^3 * N^{\alpha - 9/5} * Smoothed Heaviside Function with Alpha = ', num2str(alpha)],'Interpreter', 'tex','FontSize',16);
+% title('Cost Function: OD^3 * N^{\alpha - 9/5} with Smoothed Heaviside Function');
+%%
+% Parameters
+N = linspace(1e3, 1e6, 100);  % Atom number range
+OD = linspace(0, 3, 100);    % Optical depth range
+N1 = 1e5;                     % Example parameter for N1 controlling transition steepness
+
+% Smoothed Heaviside function as given in the appendix
+
+
+% Create figure and initial plot
+fig = figure(2);
+clf
+set(gcf,'Position',[50 50 950 750])
+set(gca,'FontSize',16,'Box','On','Linewidth',2);
+alpha = 0.5;  % Initial value for alpha
+
+% Nested function to calculate and update the cost function plot
+function update_plot(alpha_value, N, N1, OD)
+    SmoothedHeaviside = @(x) (x > 0) .* (2 ./ (1 + exp(N1 ./ x)));
+    % Calculate the cost function for the given alpha
+    CostFunction = zeros(length(N), length(OD));
+    for i = 1:length(N)
+        for j = 1:length(OD)
+            H_N = SmoothedHeaviside(N(i));  % Smoothed Heaviside function
+            CostFunction(i, j) = H_N * OD(j)^3 * N(i)^(alpha_value - 9/5);
+        end
+    end
+    % Create a meshgrid for plotting
+    [OD_grid, N_grid] = meshgrid(OD, N);
+    % Plot the cost function as a surface plot
+    surf(OD_grid, N_grid, CostFunction);
+    xlabel('Optical Depth (OD)','FontSize',16);
+    ylabel('Number of Atoms (N)','FontSize',16);
+    zlabel('Cost Function','FontSize',16);
+    title(['|C(X)| = OD^3 * N^{\alpha - 9/5} * Smoothed Heaviside Function with Alpha = ', num2str(alpha_value)],'Interpreter', 'tex','FontSize',16);
+end
+
+% Initial plot
+update_plot(alpha, N, N1, OD);
+
+% Add a slider to control alpha
+uicontrol('Style', 'slider', ...
+    'Min', 0.1, 'Max', 2, 'Value', alpha, ...  % Set the range for alpha
+    'Position', [100 20 200 20], ...
+    'Callback', @(src, event) update_plot(src.Value, N, N1, OD));
+
+% Add a label for the slider
+uicontrol('Style', 'text', ...
+    'Position', [150 45 120 20], ...
+    'String', 'Adjust Alpha');