From 88f2554b8eb2d14b82302ae6b3dd6cfd433ce748 Mon Sep 17 00:00:00 2001
From: Karthik Chandrashekara <karthik@physi.uni-heidelberg.de>
Date: Wed, 26 Mar 2025 19:50:59 +0100
Subject: [PATCH] Added custom cylindrical cut off, yet to be fully tested,
 minor mods to plotting routines.

---
 Dipolar-Gas-Simulator/+Plotter/plotLive.m     |  4 +--
 .../+Plotter/visualizeGSWavefunction.m        |  6 ++---
 Dipolar-Gas-Simulator/+Scripts/run_locally.m  | 17 +++++++-----
 .../+Scripts/run_on_cluster.m                 |  6 ++---
 .../+Simulator/@Calculator/Calculator.m       | 26 ++++++++++++-------
 .../+Simulator/@Calculator/calculateVDk.m     | 14 ++++++++++
 .../+Simulator/@DipolarGas/DipolarGas.m       |  2 +-
 7 files changed, 51 insertions(+), 24 deletions(-)

diff --git a/Dipolar-Gas-Simulator/+Plotter/plotLive.m b/Dipolar-Gas-Simulator/+Plotter/plotLive.m
index cd033be..bec2d08 100644
--- a/Dipolar-Gas-Simulator/+Plotter/plotLive.m
+++ b/Dipolar-Gas-Simulator/+Plotter/plotLive.m
@@ -32,7 +32,7 @@ function plotLive(psi,Params,Transf,Observ)
     colormap(gca, Helper.Colormaps.plasma())
     xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
     ylabel('$z$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
-    title('$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14) 
+    title('$|\Psi(x,z)|^2$', 'Interpreter', 'latex', 'FontSize', 14) 
 
     nexttile;
     plotyz = pcolor(y,z,nyz');
@@ -43,7 +43,7 @@ function plotLive(psi,Params,Transf,Observ)
     colormap(gca, Helper.Colormaps.plasma())
     xlabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
     ylabel('$z$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
-    title('$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14) 
+    title('$|\Psi(y,z)|^2$', 'Interpreter', 'latex', 'FontSize', 14) 
 
     nexttile;
     plotxy = pcolor(x,y,nxy');
diff --git a/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m b/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m
index 8071c95..dbe589a 100644
--- a/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m
+++ b/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m
@@ -1,9 +1,9 @@
-function visualizeGSWavefunction(run_index)
+function visualizeGSWavefunction(folder_path, run_index)
+    
     set(0,'defaulttextInterpreter','latex')
     set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
     
-    folder_path = './Data';
-
+    % Load data
     Data   = load(sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'),run_index),'psi','Params','Transf','Observ');
     
     Params = Data.Params;
diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m
index 7736548..8c90eee 100644
--- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m
+++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m
@@ -11,14 +11,17 @@ OptionsStruct.DipolarPolarAngle       = deg2rad(0);
 OptionsStruct.DipolarAzimuthAngle     = 0;
 OptionsStruct.ScatteringLength        = 85;
 
-OptionsStruct.TrapFrequencies         = [125, 125, 350];
+AspectRatio                           = 2.8;
+HorizontalTrapFrequency               = 125;
+VerticalTrapFrequency                 = AspectRatio * HorizontalTrapFrequency;
+OptionsStruct.TrapFrequencies         = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency];
 OptionsStruct.TrapPotentialType       = 'Harmonic';
 
-OptionsStruct.NumberOfGridPoints      = [128, 128, 64];
+OptionsStruct.NumberOfGridPoints      = [256, 256, 128];
 OptionsStruct.Dimensions              = [30, 30, 15];
 OptionsStruct.CutoffType              = 'Cylindrical';
 OptionsStruct.SimulationMode          = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution'
-OptionsStruct.TimeStepSize            = 0.002;           % in s
+OptionsStruct.TimeStepSize            = 0.001;           % in s
 OptionsStruct.MinimumTimeStepSize     = 1E-6;            % in s
 OptionsStruct.TimeCutOff              = 1E6;             % in s
 OptionsStruct.EnergyTolerance         = 5E-10;
@@ -26,10 +29,10 @@ OptionsStruct.ResidualTolerance       = 1E-05;
 OptionsStruct.NoiseScaleFactor        = 0.05;
 
 OptionsStruct.PlotLive                = true;
-OptionsStruct.JobNumber               = 1;
+OptionsStruct.JobNumber               = 0;
 OptionsStruct.RunOnGPU                = false;
 OptionsStruct.SaveData                = true;
-OptionsStruct.SaveDirectory           = './Results/Data_3D';
+OptionsStruct.SaveDirectory           = sprintf('./Results/Data_3D/AspectRatio%s', strrep(num2str(AspectRatio), '.', '_'));
 options                               = Helper.convertstruct2cell(OptionsStruct);
 clear OptionsStruct
 
@@ -47,7 +50,9 @@ Plotter.visualizeTrapPotential(sim.Potential,Params,Transf)
 %% - Plot initial wavefunction
 Plotter.visualizeWavefunction(psi,Params,Transf)
 %% - Plot GS wavefunction
-Plotter.visualizeGSWavefunction(Params.njob)
+SaveDirectory = './Results/Data_3D/AspectRatio2_8';
+JobNumber     = 0;
+Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
 
 %%
 
diff --git a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m
index 5f558fd..fa7956d 100644
--- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m
+++ b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m
@@ -17,7 +17,7 @@ OptionsStruct.NumberOfGridPoints      = [256, 256, 128];
 OptionsStruct.Dimensions              = [30, 30, 15];
 OptionsStruct.CutoffType              = 'Cylindrical';
 OptionsStruct.SimulationMode          = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution'
-OptionsStruct.TimeStepSize            = 0.002;           % in s
+OptionsStruct.TimeStepSize            = 0.001;           % in s
 OptionsStruct.MinimumTimeStepSize     = 1E-6;            % in s
 OptionsStruct.TimeCutOff              = 1E6;             % in s
 OptionsStruct.EnergyTolerance         = 5E-10;
@@ -58,7 +58,7 @@ OptionsStruct.NumberOfGridPoints      = [256, 256, 128];
 OptionsStruct.Dimensions              = [30, 30, 15];
 OptionsStruct.CutoffType              = 'Cylindrical';
 OptionsStruct.SimulationMode          = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution'
-OptionsStruct.TimeStepSize            = 0.002;           % in s
+OptionsStruct.TimeStepSize            = 0.001;           % in s
 OptionsStruct.MinimumTimeStepSize     = 1E-6;            % in s
 OptionsStruct.TimeCutOff              = 1E6;             % in s
 OptionsStruct.EnergyTolerance         = 5E-10;
@@ -66,7 +66,7 @@ OptionsStruct.ResidualTolerance       = 1E-05;
 OptionsStruct.NoiseScaleFactor        = 0.05;
 
 OptionsStruct.PlotLive                = false;
-OptionsStruct.JobNumber               = 1;
+OptionsStruct.JobNumber               = 0;
 OptionsStruct.RunOnGPU                = true;
 OptionsStruct.SaveData                = true;
 OptionsStruct.SaveDirectory           = sprintf('./Results/Data_3D/AspectRatio%s', strrep(num2str(AspectRatio), '.', '_'));
diff --git a/Dipolar-Gas-Simulator/+Simulator/@Calculator/Calculator.m b/Dipolar-Gas-Simulator/+Simulator/@Calculator/Calculator.m
index a80c05c..638d6a8 100644
--- a/Dipolar-Gas-Simulator/+Simulator/@Calculator/Calculator.m
+++ b/Dipolar-Gas-Simulator/+Simulator/@Calculator/Calculator.m
@@ -8,7 +8,9 @@ classdef Calculator < handle & matlab.mixin.Copyable
                                         'OrderParameter',                                          1, ...
                                         'PhaseCoherence',                                          1, ...
                                         'TotalEnergy',                                             1, ...
-                                        'CutoffType',                                  'Cylindrical');
+                                        'CutoffType',                                  'Cylindrical', ...
+                                        'CustomCylindricalCutOffRadius',                           1, ...
+                                        'CustomCylindricalCutOffHeight',                           1);
 
     end
 
@@ -20,6 +22,8 @@ classdef Calculator < handle & matlab.mixin.Copyable
         PhaseCoherence;
         TotalEnergy;
         CutoffType;
+        CustomCylindricalCutOffRadius
+        CustomCylindricalCutOffHeight
     end
     
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -31,7 +35,7 @@ classdef Calculator < handle & matlab.mixin.Copyable
             p = inputParser;
             p.KeepUnmatched = true;
             addParameter(p, 'CutoffType', this.CalculatorDefaults.CutoffType,... 
-                @(x) any(strcmpi(x,{'Cylindrical','CylindricalInfiniteZ', 'Spherical'})));
+                @(x) any(strcmpi(x,{'Cylindrical','CylindricalInfiniteZ', 'Spherical', 'CustomCylindrical'})));
             
             p.parse(varargin{:});
             
@@ -42,16 +46,20 @@ classdef Calculator < handle & matlab.mixin.Copyable
             this.PhaseCoherence                = this.CalculatorDefaults.PhaseCoherence;
             this.TotalEnergy                   = this.CalculatorDefaults.TotalEnergy;
             this.CutoffType                    = p.Results.CutoffType;
+            this.CustomCylindricalCutOffRadius = this.CalculatorDefaults.CustomCylindricalCutOffRadius;
+            this.CustomCylindricalCutOffHeight = this.CalculatorDefaults.CustomCylindricalCutOffHeight;
         end
 
         function restoreDefaults(this)
-            this.ChemicalPotential          = this.CalculatorDefaults.ChemicalPotential;
-            this.EnergyComponents           = this.CalculatorDefaults.EnergyComponents;
-            this.NormalizedResiduals        = this.CalculatorDefaults.NormalizedResiduals;
-            this.OrderParameter             = this.CalculatorDefaults.OrderParameter;
-            this.PhaseCoherence             = this.CalculatorDefaults.PhaseCoherence;
-            this.TotalEnergy                = this.CalculatorDefaults.TotalEnergy;
-            this.CutoffType                 = this.CalculatorDefaults.CutoffType;
+            this.ChemicalPotential             = this.CalculatorDefaults.ChemicalPotential;
+            this.EnergyComponents              = this.CalculatorDefaults.EnergyComponents;
+            this.NormalizedResiduals           = this.CalculatorDefaults.NormalizedResiduals;
+            this.OrderParameter                = this.CalculatorDefaults.OrderParameter;
+            this.PhaseCoherence                = this.CalculatorDefaults.PhaseCoherence;
+            this.TotalEnergy                   = this.CalculatorDefaults.TotalEnergy;
+            this.CutoffType                    = this.CalculatorDefaults.CutoffType;
+            this.CustomCylindricalCutOffRadius = this.CalculatorDefaults.CustomCylindricalCutOffRadius;
+            this.CustomCylindricalCutOffHeight = this.CalculatorDefaults.CustomCylindricalCutOffHeight;
         end
         
     end %
diff --git a/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateVDk.m b/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateVDk.m
index 1c9fbc7..84bff55 100644
--- a/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateVDk.m
+++ b/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateVDk.m
@@ -57,6 +57,20 @@ function VDk = calculateVDk(this,Params,Transf,TransfRad,IncludeDDICutOff)
         
                 K       = sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
                 VDk     = (cos(alph).^2-1/3).*(1 + 3*cos(Rcut*K)./(Rcut^2.*K.^2) - 3*sin(Rcut*K)./(Rcut^3.*K.^3));
+                
+            case 'CustomCylindrical'
+                alph    = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
+                alph(1) = pi/2;
+                VDk     = cos(alph).^2-1/3;
+                % Implementing a custom cylindrical cutoff:
+                VDr     = ifftn(VDk);
+                VDr     = fftshift(VDr);
+                % Define the cylindrical mask
+                cylinder_mask = (Transf.X.^2 + Transf.Y.^2 <= this.CustomCylindricalCutOffRadius^2) & (abs(Transf.Z) <= this.CustomCylindricalCutOffHeight / 2);
+                % Multiply the mask to apply the cut-off
+                VDr     = VDr.*double(cylinder_mask);
+                VDr     = ifftshift(VDr);
+                VDk     = fftn(VDr);
             
             otherwise
                 disp('Choose a valid DDI cutoff type!')
diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m
index f3b8693..292a56c 100644
--- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m
+++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m
@@ -55,7 +55,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
             addParameter(p, 'SimulationMode', 'ImaginaryTimeEvolution',... 
                 @(x) assert(any(strcmpi(x,{'ImaginaryTimeEvolution','RealTimeEvolution'}))));
             addParameter(p, 'CutoffType', 'Cylindrical',... 
-                @(x) assert(any(strcmpi(x,{'Cylindrical','CylindricalInfiniteZ', 'Spherical'}))));
+                @(x) assert(any(strcmpi(x,{'Cylindrical','CylindricalInfiniteZ', 'Spherical', 'CustomCylindrical'}))));
             addParameter(p, 'TimeStepSize',         5E-4,...
                 @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));  
             addParameter(p, 'MinimumTimeStepSize',   1e-6,...