%% This script is testing the functionalities of the Dipolar Gas Simulator
%
%  Important:  Run only sectionwise!! 

%% - Create Simulator, Potential and Calculator object with specified options

OptionsStruct = struct;

OptionsStruct.NumberOfAtoms           = 40000;
OptionsStruct.DipolarPolarAngle       = deg2rad(0);
OptionsStruct.DipolarAzimuthAngle     = 0;
OptionsStruct.ScatteringLength        = 95;

OptionsStruct.TrapFrequencies         = [30, 60, 90];
OptionsStruct.TrapPotentialType       = 'Harmonic';

OptionsStruct.NumberOfGridPoints      = [128, 64, 64];
OptionsStruct.Dimensions              = [30, 20, 20];
OptionsStruct.UseApproximationForLHY  = true;
OptionsStruct.IncludeDDICutOff        = true;
OptionsStruct.CutoffType              = 'Cylindrical';
OptionsStruct.SimulationMode          = 'EnergyMinimization';     % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization'
OptionsStruct.GradientDescentMethod   = 'NonLinearCGD';           % 'HeavyBall' | 'NonLinearCGD' 
OptionsStruct.MaxIterationsForGD      = 1000;
OptionsStruct.TimeStepSize            = 1E-3;            % in s
OptionsStruct.MinimumTimeStepSize     = 1E-6;            % in s
OptionsStruct.TimeCutOff              = 2E6;             % in s
OptionsStruct.EnergyTolerance         = 5E-10;
OptionsStruct.ResidualTolerance       = 1E-08;
OptionsStruct.NoiseScaleFactor        = 0.010;

OptionsStruct.PlotLive                = true;
OptionsStruct.JobNumber               = 0;
OptionsStruct.RunOnGPU                = false;
OptionsStruct.SaveData                = true;
OptionsStruct.SaveDirectory           = './Results/Data_Full3D/GradientDescent'; % './Results/Data_Full3D/AnisotropicTrap/Tilted0'
options                               = Helper.convertstruct2cell(OptionsStruct);

sim                                   = Simulator.DipolarGas(options{:});
pot                                   = Simulator.Potentials(options{:});
sim.Potential                         = pot.trap(); % + pot.repulsive_chopstick();

%-% Run Simulation %-%
NumberOfOutputs                       = 5;
[Params, Transf, psi, V, VDk, stats]  = Helper.runWithProfiling(@() sim.run(), NumberOfOutputs, OptionsStruct.SaveDirectory);
fprintf('Runtime: %.3f seconds\n', stats.runtime);
fprintf('Memory used: %.2f MB\n', stats.workspaceMemoryMB);

%% - Plot numerical grid
% Plotter.visualizeSpace(Transf)
%% - Plot trap potential
Plotter.visualizeTrapPotential(sim.Potential,Params,Transf)
%% - Plot initial wavefunction
Plotter.visualizeWavefunction(psi,Params,Transf)

%% Imaginary-Time followed by Real-time

% - Imaginary-Time
OptionsStruct = struct;

OptionsStruct.NumberOfAtoms                 = 1E5;
OptionsStruct.DipolarPolarAngle             = deg2rad(0);
OptionsStruct.DipolarAzimuthAngle           = deg2rad(0);
OptionsStruct.ScatteringLength              = 110;

OptionsStruct.TrapFrequencies               = [50, 20, 150];
OptionsStruct.TrapPotentialType             = 'Harmonic';

OptionsStruct.NumberOfGridPoints            = [64, 128, 32];
OptionsStruct.Dimensions                    = [50, 50, 10];
OptionsStruct.UseApproximationForLHY        = true;
OptionsStruct.IncludeDDICutOff              = true;
OptionsStruct.CutoffType                    = 'CustomCylindrical';
OptionsStruct.CustomCylindricalCutOffRadius = 20.0;
OptionsStruct.CustomCylindricalCutOffHeight = 4.0;

OptionsStruct.SimulationMode                = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization'
OptionsStruct.TimeStepSize                  = 5E-3;           % in s
OptionsStruct.MinimumTimeStepSize           = 1E-6;           % in s
OptionsStruct.TimeCutOff                    = 1E2;            % in s
OptionsStruct.EnergyTolerance               = 5E-10;
OptionsStruct.ResidualTolerance             = 1E-08;
OptionsStruct.NoiseScaleFactor              = 0.05;
OptionsStruct.PlotLive                      = true;
OptionsStruct.JobNumber                     = 0;
OptionsStruct.RunOnGPU                      = false;
OptionsStruct.SaveData                      = false;
OptionsStruct.SaveDirectory                 = './Results/Data_Full3D/RealTimeDynamics';
options                                     = Helper.convertstruct2cell(OptionsStruct);

sim                                         = Simulator.DipolarGas(options{:});
pot                                         = Simulator.Potentials(options{:});
sim.Potential                               = pot.trap();

%-% Run Simulation %-%
[Params, Transf, psi, V, VDk]               = sim.run();

% Save only final state as initial wavefunction for real-time propagation
mkdir(sprintf(OptionsStruct.SaveDirectory))
save(sprintf(strcat(OptionsStruct.SaveDirectory, '/psi_init.mat'),Params.njob),'psi','Transf','Params','VDk','V');

OptionsStruct.SimulationMode                = 'RealTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization'
OptionsStruct.EquilibrationTime             = 10E-3; 
OptionsStruct.QuenchTime                    = 10E-3;
OptionsStruct.HoldTime                      = 10E-3;
OptionsStruct.QuenchScatteringLength        = false;
OptionsStruct.RotateDipoles                 = true;
OptionsStruct.FinalScatteringLength         = 88;
OptionsStruct.FinalDipolarPolarAngle        = deg2rad(35);
OptionsStruct.FinalDipolarAzimuthAngle      = deg2rad(0);
OptionsStruct.TimeStepSize                  = 1E-3;            % in s
OptionsStruct.NoiseScaleFactor              = 0.010;

OptionsStruct.PlotLive                      = true;
OptionsStruct.JobNumber                     = 0;
OptionsStruct.RunOnGPU                      = false;
OptionsStruct.SaveData                      = true;
options                                     = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct

sim                                         = Simulator.DipolarGas(options{:});
pot                                         = Simulator.Potentials(options{:});
sim.Potential                               = pot.trap(); 

%-% Run Simulation %-%
[Params, Transf, psi, V, VDk]               = sim.run();

%%
% To reproduce results from the Blair Blakie paper:
% (n*add^2, as/add)
% Critical point: (0.0978, 0.784); Triangular phase: (0.0959, 0.750); Stripe phase: (0.144, 0.765); Honeycomb phase: (0.192, 0.780)

% N  = ((n*add^2)/Params.add^2) * (Params.Lx *1E-6)^2
% Critical point: N = 2.0427e+07; Triangular phase: N = 2.0030e+07; Stripe phase: N = 3.0077e+07; Honeycomb phase: N = 4.0102e+07 for dimensions fixed to 100

% as = ((as/add)*Params.add)/Params.a0
% Critical point: 102.5133; Triangular phase: 98.0676; Stripe phase: 100.0289; Honeycomb phase: 101.9903

%% - Create Variational2D and Calculator object with specified options

OptionsStruct = struct;

OptionsStruct.NumberOfAtoms           = 2.0030e+07;     
OptionsStruct.DipolarPolarAngle       = 0;
OptionsStruct.DipolarAzimuthAngle     = 0;
OptionsStruct.ScatteringLength        = 98.0676;        

OptionsStruct.TrapFrequencies         = [0, 0, 72.4];
OptionsStruct.TrapPotentialType       = 'None';

OptionsStruct.NumberOfGridPoints      = [128, 128];
OptionsStruct.Dimensions              = [100, 100];      
OptionsStruct.TimeStepSize            = 0.005;           % in s
OptionsStruct.MinimumTimeStepSize     = 1E-5;            % in s
OptionsStruct.TimeCutOff              = 1E6;             % in s
OptionsStruct.EnergyTolerance         = 5E-10;
OptionsStruct.ResidualTolerance       = 1E-05;
OptionsStruct.NoiseScaleFactor        = 0.05;

OptionsStruct.MaxIterations           = 10;      
OptionsStruct.VariationalWidth        = 5.0;       
OptionsStruct.WidthLowerBound         = 1;
OptionsStruct.WidthUpperBound         = 12;
OptionsStruct.WidthCutoff             = 5e-3;

OptionsStruct.PlotLive                = true;
OptionsStruct.JobNumber               = 1;
OptionsStruct.RunOnGPU                = false;
OptionsStruct.SaveData                = true;
OptionsStruct.SaveDirectory           = './Results/Data_TriangularPhase';
options                               = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct

solver                                = VariationalSolver2D.DipolarGas(options{:});
pot                                   = VariationalSolver2D.Potentials(options{:});
solver.Potential                      = pot.trap();

%-% Run Solver %-%
[Params, Transf, psi, V, VDk]         = solver.run();
%%
% Solve BdG equations
% Load data
Data = load(sprintf(horzcat(solver.SaveDirectory, '/Run_%03i/psi_gs.mat'),solver.JobNumber),'psi','Transf','Params','VParams','V');

Transf  = Data.Transf; 
Params  = Data.Params;
VParams = Data.VParams;
V       = Data.V;

if isgpuarray(Data.psi)
    psi = gather(Data.psi);
else
    psi = Data.psi;
end

% == DDI potential == %
VDk = solver.Calculator.calculateVDkWithCutoff(Transf, Params, VParams.ell);

% == Chemical potential == %
muchem = solver.Calculator.calculateChemicalPotential(psi,Params,VParams,Transf,VDk,V);

[evals, modes] = BdGSolver2D.solveBogoliubovdeGennesIn2D(psi, Params, VDk, VParams, Transf, muchem);

% Save the eigenvalues and eigenvectors to a .mat file
save(sprintf(strcat(solver.SaveDirectory, '/Run_%03i/bdg_eigen_data.mat'),solver.JobNumber), 'evals', 'modes', '-v7.3');

%% - Create Variational2D and Calculator object with specified options

OptionsStruct = struct;

OptionsStruct.NumberOfAtoms           = 3.0077e+07;     
OptionsStruct.DipolarPolarAngle       = 0;
OptionsStruct.DipolarAzimuthAngle     = 0;
OptionsStruct.ScatteringLength        = 100.0289;        

OptionsStruct.TrapFrequencies         = [0, 0, 72.4];
OptionsStruct.TrapPotentialType       = 'None';

OptionsStruct.NumberOfGridPoints      = [128, 128];
OptionsStruct.Dimensions              = [100, 100];      
OptionsStruct.TimeStepSize            = 0.005;           % in s
OptionsStruct.MinimumTimeStepSize     = 1E-5;            % in s
OptionsStruct.TimeCutOff              = 1E6;             % in s
OptionsStruct.EnergyTolerance         = 5E-10;
OptionsStruct.ResidualTolerance       = 1E-05;
OptionsStruct.NoiseScaleFactor        = 0.05;

OptionsStruct.MaxIterations           = 10;      
OptionsStruct.VariationalWidth        = 5;       
OptionsStruct.WidthLowerBound         = 1;
OptionsStruct.WidthUpperBound         = 12;
OptionsStruct.WidthCutoff             = 5e-3;

OptionsStruct.PlotLive                = true;
OptionsStruct.JobNumber               = 1;
OptionsStruct.RunOnGPU                = false;
OptionsStruct.SaveData                = true;
OptionsStruct.SaveDirectory           = './Results/Data_StripePhase';
options                               = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct

solver                                = VariationalSolver2D.DipolarGas(options{:});
pot                                   = VariationalSolver2D.Potentials(options{:});
solver.Potential                      = pot.trap();

%-% Run Solver %-%
[Params, Transf, psi, V, VDk]         = solver.run();

%% - Create Variational2D and Calculator object with specified options

OptionsStruct = struct;

OptionsStruct.NumberOfAtoms           = 4.2e+07;     
OptionsStruct.DipolarPolarAngle       = 0;
OptionsStruct.DipolarAzimuthAngle     = 0;
OptionsStruct.ScatteringLength        = 101.35;          

OptionsStruct.TrapFrequencies         = [0, 0, 72.4];
OptionsStruct.TrapPotentialType       = 'None';

OptionsStruct.NumberOfGridPoints      = [128, 128];
OptionsStruct.Dimensions              = [100, 100];      
OptionsStruct.TimeStepSize            = 0.005;           % in s
OptionsStruct.MinimumTimeStepSize     = 1E-5;            % in s
OptionsStruct.TimeCutOff              = 2E6;             % in s
OptionsStruct.EnergyTolerance         = 5E-10;
OptionsStruct.ResidualTolerance       = 1E-05;
OptionsStruct.NoiseScaleFactor        = 0.05;

OptionsStruct.MaxIterations           = 10;      
OptionsStruct.VariationalWidth        = 6;       
OptionsStruct.WidthLowerBound         = 1;
OptionsStruct.WidthUpperBound         = 12;
OptionsStruct.WidthCutoff             = 5e-3;

OptionsStruct.PlotLive                = true;
OptionsStruct.JobNumber               = 1;
OptionsStruct.RunOnGPU                = false;
OptionsStruct.SaveData                = true;
OptionsStruct.SaveDirectory           = './Results/Data_HoneycombPhase';
options                               = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct

solver                                = VariationalSolver2D.DipolarGas(options{:});
pot                                   = VariationalSolver2D.Potentials(options{:});
solver.Potential                      = pot.trap();

%-% Run Solver %-%
[Params, Transf, psi, V, VDk]         = solver.run();

%% - Plot numerical grid
% Plotter.visualizeSpace2D(Transf)
%% - Plot trap potential
% Plotter.visualizeTrapPotential2D(solver.Potential,Params,Transf)
%% - Plot initial wavefunction
Plotter.visualizeWavefunction2D(psi,Params,Transf)
%% - Plot GS wavefunction
Plotter.visualizeGSWavefunction2D(solver.SaveDirectory, solver.JobNumber)

%% - Analysis
SaveDirectory = './Results/Data_TriangularPhase';
JobNumber     = 1;
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)

%% - Analysis
SaveDirectory = './Results/Data_StripePhase';
JobNumber     = 2;
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)

%% - Analysis
SaveDirectory = './Results/Data_HoneycombPhase';
JobNumber     = 3;
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)

%% Tilting of the dipoles
% Atom Number      = 1.00e+05
% System size      = [10, 10]

% - Analysis
SaveDirectory = './Results/Data_TiltingOfDipoles';
JobNumber     = 5;
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)

%% Tilting of the dipoles
% Atom Number      = 1.00e+05
% System size      = [5 * l_rot, 5 * l_rot]

%% - Analysis
SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz500';
JobNumber     = 2;
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction(SaveDirectory, JobNumber);

%% - Analysis
SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz750';
JobNumber     = 1;
% Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction(SaveDirectory, JobNumber);

%% - Analysis
SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz1000';
JobNumber     = 1;
% Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction(SaveDirectory, JobNumber);

%% - Analysis
SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz2000';
JobNumber     = 1;
% Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction(SaveDirectory, JobNumber);

%% - Analysis
SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/Hz500';
NumberOfJobs  = 16;
SaveOption    = 'video';
[contrast_array, periodX_array, periodY_array] = Scripts.analyzeGSWavefunction_multipleruns(SaveDirectory, NumberOfJobs, SaveOption);

%% - Analysis
SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/Hz750';
NumberOfJobs  = 16;
SaveOption    = 'video';
[contrast_array, periodX_array, periodY_array] = Scripts.analyzeGSWavefunction_multipleruns(SaveDirectory, NumberOfJobs, SaveOption);

%% - Analysis
SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/Hz1000';
NumberOfJobs  = 16;
SaveOption    = 'video';
[contrast_array, periodX_array, periodY_array] = Scripts.analyzeGSWavefunction_multipleruns(SaveDirectory, NumberOfJobs, SaveOption);

%% - Analysis
SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/Hz2000';
NumberOfJobs  = 16;
SaveOption    = 'video';
[contrast_array, periodX_array, periodY_array] = Scripts.analyzeGSWavefunction_multipleruns(SaveDirectory, NumberOfJobs, SaveOption);
%% - Analysis
SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle_newruns/Hz1000';
JobNumber     = 0;
% Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction(SaveDirectory, JobNumber);
%% - Analysis
SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/Hz500';
JobNumber     = 0;
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber);

%% - Analysis
SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/Hz750';
JobNumber     = 1;
% Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber);

%% - Analysis
SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/Hz1000';
JobNumber     = 1;
% Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber);

%% - Analysis
SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/Hz2000';
JobNumber     = 1;
% Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber);

%%
SaveDirectory = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSize/Hz500/Degree0';
% Define the desired range for Lx and Ly
Lx_Range      = 4:0.2:6; % Extend Lx from 5 to 10
Ly_Range      = 4:0.2:6; % Extend Ly from 5 to 10
MinEnergyDataArray = Scripts.analyzeGSWavefunction_optimal_system_size(SaveDirectory, Lx_Range, Ly_Range);

%%
% SaveDirectory      = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSize/Hz500/Degree0';
SaveDirectory      = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSizeWithBiasAnsatz/Hz500/Degree0';
% Define the desired range
% LatticeSpacing     = 1.0:0.05:4.0; 
LatticeSpacing     = linspace(1.6,2.25,61); 
MinEnergyDataArray = Scripts.analyzeGSWavefunction_constrained_optimal_system_size(SaveDirectory, LatticeSpacing);
[val, idx]         = min(MinEnergyDataArray(:,3));
OptimalSize        = MinEnergyDataArray(idx,1:2);
Plotter.visualizeGSWavefunction2D(SaveDirectory, idx)
[contrast, periodX, periodY] = Scripts.analyzeGSWavefunction(SaveDirectory, idx, false);
%%
% SaveDirectory      = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSize/Hz500/Degree5';
SaveDirectory      = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSizeWithBiasAnsatz/Hz500/Degree5';
% Define the desired range
% LatticeSpacing     = 1.0:0.05:4.0; 
LatticeSpacing     = linspace(1.6,2.25,61); 
MinEnergyDataArray = Scripts.analyzeGSWavefunction_constrained_optimal_system_size(SaveDirectory, LatticeSpacing);
[val, idx]         = min(MinEnergyDataArray(:,3));
OptimalSize        = MinEnergyDataArray(idx,1:2);
Plotter.visualizeGSWavefunction2D(SaveDirectory, idx)
[contrast, periodX, periodY] = Scripts.analyzeGSWavefunction(SaveDirectory, idx, false);

%%
% SaveDirectory      = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSize/Hz500/Degree7_5';
SaveDirectory      = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSizeWithBiasAnsatz/Hz500/Degree7_5';
% Define the desired range
% LatticeSpacing     = 1.0:0.05:4.0; 
LatticeSpacing     = linspace(1.6,2.25,61); 
MinEnergyDataArray = Scripts.analyzeGSWavefunction_constrained_optimal_system_size(SaveDirectory, LatticeSpacing);
[val, idx]         = min(MinEnergyDataArray(:,3));
OptimalSize        = MinEnergyDataArray(idx,1:2);
Plotter.visualizeGSWavefunction2D(SaveDirectory, idx)
[contrast, periodX, periodY] = Scripts.analyzeGSWavefunction(SaveDirectory, idx, false);

%%
% SaveDirectory      = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSize/Hz500/Degree10';
SaveDirectory      = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSizeWithBiasAnsatz/Hz500/Degree10';
% Define the desired range
% LatticeSpacing     = 1.0:0.05:4.0; 
LatticeSpacing     = linspace(1.6,2.25,61);
MinEnergyDataArray = Scripts.analyzeGSWavefunction_constrained_optimal_system_size(SaveDirectory, LatticeSpacing);
[val, idx]         = min(MinEnergyDataArray(:,3));
OptimalSize        = MinEnergyDataArray(idx,1:2);
Plotter.visualizeGSWavefunction2D(SaveDirectory, idx)
[contrast, periodX, periodY] = Scripts.analyzeGSWavefunction(SaveDirectory, idx, false);

%%
% SaveDirectory      = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSize/Hz500/Degree15';
SaveDirectory      = './Results/Data_TiltingOfDipoles/TransitionAngle/OptimalSystemSizeWithBiasAnsatz/Hz500/Degree15';
% Define the desired range
% LatticeSpacing     = 1.0:0.05:4.0; 
LatticeSpacing     = linspace(1.6,2.25,61);
MinEnergyDataArray = Scripts.analyzeGSWavefunction_constrained_optimal_system_size(SaveDirectory, LatticeSpacing);
[val, idx]         = min(MinEnergyDataArray(:,3));
OptimalSize        = MinEnergyDataArray(idx,1:2);
Plotter.visualizeGSWavefunction2D(SaveDirectory, idx)
[contrast, periodX, periodY] = Scripts.analyzeGSWavefunction(SaveDirectory, idx, false);

%% - Analysis
SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/AspectRatio/AR2_8';
JobNumber     = 0; % 79
% JobNumber     = 1; % 80
% JobNumber     = 2; % 81
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber);

%% - Analysis
SaveDirectory = './Results/Data_TiltingOfDipoles/HarmonicTrap/AspectRatio/AR3_7';
% JobNumber     = 0; % 80
% JobNumber     = 1; % 81
JobNumber     = 2; % 82
% JobNumber     = 3; % 83
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber);

%% - Plot GS wavefunction
% SaveDirectory = './Results/Data_Full3D/CompleteLHY/AspectRatio2_8';
% SaveDirectory = './Results/Data_Full3D/CompleteLHY/AspectRatio3_7';
SaveDirectory = './Results/Data_Full3D/CompleteLHY/AspectRatio3_8';
JobNumber     = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)

%%
% SaveDirectory = './Results/Data_Full3D/CompleteLHY/BeyondSSD_SSD';
% SaveDirectory = './Results/Data_Full3D/CompleteLHY/BeyondSSD_Labyrinth';
SaveDirectory = './Results/Data_Full3D/CompleteLHY/BeyondSSD_Honeycomb';
JobNumber     = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)

%% - Plot GS wavefunction
SaveDirectory = './Results/Data_Full3D/ApproximateLHY/AspectRatio2_8';
% SaveDirectory = './Results/Data_Full3D/ApproximateLHY/AspectRatio3_7';
% SaveDirectory = './Results/Data_Full3D/ApproximateLHY/AspectRatio3_8';
% SaveDirectory = './Results/Data_Full3D/ApproximateLHY/AspectRatio3_9';
JobNumber     = 3;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)

%%
% SaveDirectory = './Results/Data_Full3D/ApproximateLHY/BeyondSSD_SSD';
% SaveDirectory = './Results/Data_Full3D/ApproximateLHY/BeyondSSD_Labyrinth';
SaveDirectory = './Results/Data_Full3D/ApproximateLHY/BeyondSSD_Honeycomb';
JobNumber     = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)

%%
SaveDirectory = './Results/Data_Full3D/TiltedDipoles0';
JobNumber     = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)

%%
SaveDirectory = './Results/Data_Full3D/TiltedDipoles15';
JobNumber     = 2;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)

%%
SaveDirectory = './Results/Data_Full3D/TiltedDipoles30';
JobNumber     = 2;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)

%%
SaveDirectory = './Results/Data_Full3D/TiltedDipoles45';
JobNumber     = 2;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)

%%
SaveDirectory = './Results/Data_Full3D/RealTimeDynamics';
JobNumber     = 0;
Plotter.makeMovie(SaveDirectory, JobNumber)

%%
SaveDirectory = './Results/Data_Full3D/AnisotropicTrap/TiltedDipoles15';
JobNumber     = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)

%%
SaveDirectory = './Results/Data_Full3D/AnisotropicTrap/TiltedDipoles30';
JobNumber     = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
SaveDirectory = './Results/Data_Full3D/AnisotropicTrap/TiltedDipoles35';
JobNumber     = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)

%%
SaveDirectory = './Results/Data_Full3D/AnisotropicTrap/TiltedDipoles40';
JobNumber     = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)

%%
SaveDirectory = './Results/Data_Full3D/AnisotropicTrap/TiltedDipoles45';
JobNumber     = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)

%%
SaveDirectory = './Results/Data_Full3D/PhaseDiagram/GradientDescent/';
JobNumber     = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)

%%
SaveDirectory = './Results/Data_Full3D/PhaseDiagram/ImagTimePropagation/aS_8.312000e+01_theta_020_phi_000_N_100000';
JobNumber     = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)

%%
SaveDirectory = 'D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_712500';
JobNumber     = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)

%% Identify and count droplets
Radius                  = 2;  % The radius within which peaks will be considered duplicates
PeakThreshold           = 3E3;
SaveDirectory           = 'D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_712500';
JobNumber               = 0;
SuppressPlotFlag        = false;
DropletNumber           = Scripts.identifyAndCountDroplets(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag);

%% Extract average peak column density
Radius                  = 2;  % The radius within which peaks will be considered duplicates
PeakThreshold           = 3E3;
SaveDirectory           = 'D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_712500';
JobNumber               = 0;
SuppressPlotFlag        = false;
AveragePCD               = Scripts.extractAveragePeakColumnDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag);

%% Extract average unit cell density
Radius                  = 2;  % The radius within which peaks will be considered duplicates
PeakThreshold           = 3E3;
SaveDirectory           = 'D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_9.562000e+01_theta_000_phi_000_N_712500';
JobNumber               = 0;
SuppressPlotFlag        = false;
UCD                     = Scripts.extractAverageUnitCellDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag);

%% Plot number of droplets
% Parameters
Radius                  = 2;
PeakThreshold           = 6E3;
JobNumber               = 0;
TitleString             = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 0^\circ";
SuppressPlotFlag        = true;

SCATTERING_LENGTH_RANGE = [85.21 86.25 87.29 88.33 89.38];

NUM_ATOMS_LIST          = [100000 304167 508333 712500 916667 1120833 1325000 ...
                           1529167 1733333 1937500 2141667 2345833 2550000 ...
                           2754167 2958333 3162500 3366667 3570833];

% Prepare figure
figure(1);
clf;
set(gcf,'Position', [100, 100, 1000, 700]);
hold on;

% Color order for better visibility
colors = lines(length(SCATTERING_LENGTH_RANGE));

% Store legend labels
legendEntries = cell(1, length(SCATTERING_LENGTH_RANGE));

% Loop over scattering lengths
for j = 1:length(SCATTERING_LENGTH_RANGE)
    aS = SCATTERING_LENGTH_RANGE(j);
    
    % Format scattering length in scientific notation with 6 decimal places
    aS_string = sprintf('%.6e', aS);
    
    % Construct base directory for this aS
    baseDir = ['D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_' ...
               aS_string '_theta_000_phi_000_N_'];
    
    % Preallocate results for this curve
    DropletNumbers = zeros(size(NUM_ATOMS_LIST));
    
    % Loop over all atom numbers
    for i = 1:length(NUM_ATOMS_LIST)
        N = NUM_ATOMS_LIST(i);
        SaveDirectory = [baseDir num2str(N)];
        
        % Call your droplet counting function
        DropletNumbers(i) = Scripts.identifyAndCountDroplets(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag);
    end
    
    % Plot curve
    plot(NUM_ATOMS_LIST, DropletNumbers, 'o-', ...
         'Color', colors(j, :), 'LineWidth', 1.5);
    
    % Store legend entry
    legendEntries{j} = ['a_s = ' num2str(aS) 'a_o'];
end

% Finalize plot
xlabel('Number of Atoms', 'Interpreter', 'tex', 'FontSize', 16);
ylabel('Number of Droplets', 'Interpreter', 'tex', 'FontSize', 16);
title(TitleString, 'Interpreter', 'tex', 'FontSize', 18);
legend(legendEntries, 'Interpreter', 'tex', 'FontSize', 12, 'Location', 'bestoutside');
grid on;
hold off;

%% Plot average peak column density
% Parameters
Radius                  = 2;
PeakThreshold           = 400;
JobNumber               = 0;
TitleString             = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 0^\circ";
SuppressPlotFlag        = true;

SCATTERING_LENGTH_RANGE = [95.62];

NUM_ATOMS_LIST          = [50000 54545 59091 63636 68182 72727 77273 81818 86364 90909 95455 100000 304167 508333 712500 916667 1120833 1325000 ...
                           1529167 1733333 1937500 2141667 2345833 2550000 2754167 2958333 3162500 3366667 3570833 3775000 3979167 4183333 4387500 ...
                           4591667 4795833 5000000];

% Prepare figure
figure(1);
clf;
set(gcf,'Position', [100, 100, 1000, 700]);
hold on;

% Color order for better visibility
colors = lines(length(SCATTERING_LENGTH_RANGE));

% Store legend labels
legendEntries = cell(1, length(SCATTERING_LENGTH_RANGE));

% Loop over scattering lengths
for j = 1:length(SCATTERING_LENGTH_RANGE)
    aS = SCATTERING_LENGTH_RANGE(j);
    
    % Format scattering length in scientific notation with 6 decimal places
    aS_string = sprintf('%.6e', aS);
    
    % Construct base directory for this aS
    baseDir = ['D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/aS_' ...
               aS_string '_theta_000_phi_000_N_'];
    
    % Preallocate results for this curve
    AverageCDs = zeros(size(NUM_ATOMS_LIST));
    
    % Loop over all atom numbers
    for i = 1:length(NUM_ATOMS_LIST)
        N = NUM_ATOMS_LIST(i);
        SaveDirectory = [baseDir num2str(N)];
        
        % Call your droplet counting function
        AverageCDs(i) = Scripts.extractAveragePeakColumnDensity(SaveDirectory, JobNumber, Radius, PeakThreshold, SuppressPlotFlag);
    end
    
    % Plot curve
    plot(NUM_ATOMS_LIST, AverageCDs, 'o-', ...
         'Color', colors(j, :), 'LineWidth', 1.5);
    
    % Store legend entry
    legendEntries{j} = ['a_s = ' num2str(aS) 'a_o'];
end

% Finalize plot
xlabel('Number of Atoms', 'Interpreter', 'tex', 'FontSize', 16);
ylabel('Average Peak Column Density', 'Interpreter', 'tex', 'FontSize', 16);
title(TitleString, 'Interpreter', 'tex', 'FontSize', 18);
legend(legendEntries, 'Interpreter', 'tex', 'FontSize', 12, 'Location', 'bestoutside');
grid on;
hold off;

%% Plot TF radii of unmodulated states
% Parameters
JobNumber               = 0;
TitleString             = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 40^\circ";
SuppressPlotFlag        = true;

SCATTERING_LENGTH_RANGE = [97.71];

NUM_ATOMS_LIST          = [50000 54545 59091 63636 68182 72727 77273 81818 86364 90909 95455];
% NUM_ATOMS_LIST          = [100000 304167 508333 712500 916667 1120833 1325000 1529167 1733333 1937500 2141667 2345833 2550000 2754167 2958333 3162500 3366667 3570833 3775000 3979167 4183333 4387500 4591667 4795833 5000000];

% Loop over scattering lengths
for j = 1:length(SCATTERING_LENGTH_RANGE)
    aS = SCATTERING_LENGTH_RANGE(j);
    
    % Format scattering length in scientific notation with 6 decimal places
    aS_string = sprintf('%.6e', aS);
    
    % Construct base directory for this aS
    baseDir = ['D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta40/LowN/aS_' ...
               aS_string '_theta_040_phi_000_N_'];
    
   % Preallocate results for this curve: rows = N, cols = [Rx, Ry]
   TF_Radii = zeros(length(NUM_ATOMS_LIST), 2);
    
    % Loop over all atom numbers
    for i = 1:length(NUM_ATOMS_LIST)
        N = NUM_ATOMS_LIST(i);
        SaveDirectory = [baseDir num2str(N)];
    
        % Extract both Rx and Ry
        try
            [Rx, Ry] = Scripts.extractThomasFermiRadii(SaveDirectory, JobNumber, SuppressPlotFlag);
        catch ME
            warning(ME.message)
            Rx       = NaN;
            Ry       = NaN;
        end

        % Store as row [Rx, Ry]
        TF_Radii(i, :) = [Rx, Ry];
    end
end

% Plot curve
% Prepare figure
fig = figure(1); 
clf
set(gcf,'Position', [100, 100, 1200, 500])
t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact');
% Color order for better visibility
colors = lines(length(SCATTERING_LENGTH_RANGE));
% Original density plot
nexttile;
plot(NUM_ATOMS_LIST, TF_Radii(:, 1), 'o-', ...
     'Color', colors(j, :), 'LineWidth', 1.5);
% Store legend entry
legendEntries{j} = ['a_s = ' num2str(aS) 'a_o'];

% Finalize plot
xlabel('Number of Atoms', 'Interpreter', 'tex', 'FontSize', 16);
ylabel('TF Radius - X ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16);
legend(legendEntries,'Interpreter', 'tex', 'FontSize', 12, 'Location', 'bestoutside');
axis square
grid on;

nexttile
plot(NUM_ATOMS_LIST, TF_Radii(:, 2), 'o-', ...
     'Color', colors(j, :), 'LineWidth', 1.5);
% Store legend entry
legendEntries{j} = ['a_s = ' num2str(aS) 'a_o'];

% Finalize plot
xlabel('Number of Atoms', 'Interpreter', 'tex', 'FontSize', 16);
ylabel('TF Radius - Y ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16);
legend(legendEntries,'Interpreter', 'tex', 'FontSize', 12, 'Location', 'bestoutside');
axis square
grid on;

sgtitle(TitleString, 'Interpreter', 'tex', 'FontSize', 18);

%% Overlay curves

% File list and associated theta labels
fileList = {'TFFermi_Theta0.mat', 'TFFermi_Theta20.mat', 'TFFermi_Theta40.mat'};
thetaLabels = {'\theta = 0^\circ', '\theta = 20^\circ', '\theta = 40^\circ'};

% Set up figure
fig = figure(1); 
clf
set(gcf,'Position', [100, 100, 1200, 500])
t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact');

% Color order for visibility
colors = lines(length(fileList));

legendEntries = cell(1, length(fileList));

% Loop over files
for j = 1:length(fileList)
    data = load(fileList{j});
    
    aS              = data.SCATTERING_LENGTH_RANGE;  % scalar
    NUM_ATOMS_LIST  = data.NUM_ATOMS_LIST;
    TF_Radii        = data.TF_Radii;  % [N x 2]

    % Store legend entry (can choose to include theta or a_s)
    legendEntries{j} = sprintf('%s, a_s = %.2f a_0', thetaLabels{j}, aS);

    % Plot Rx (TF_Radii(:,1))
    nexttile(1);
    plot(NUM_ATOMS_LIST, TF_Radii(:, 1), 'o-', ...
         'Color', colors(j, :), 'LineWidth', 1.5); hold on;

    % Plot Ry (TF_Radii(:,2))
    nexttile(2);
    plot(NUM_ATOMS_LIST, TF_Radii(:, 2), 'o-', ...
         'Color', colors(j, :), 'LineWidth', 1.5); hold on;
end

% Finalize Rx plot
nexttile(1);
xlabel('Number of Atoms', 'Interpreter', 'tex', 'FontSize', 16);
ylabel('TF Radius - X ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16);
legend(legendEntries, 'Interpreter', 'tex', 'FontSize', 12, 'Location', 'bestoutside');
axis square
grid on;

% Finalize Ry plot
nexttile(2);
xlabel('Number of Atoms', 'Interpreter', 'tex', 'FontSize', 16);
ylabel('TF Radius - Y ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16);
legend(legendEntries, 'Interpreter', 'tex', 'FontSize', 12, 'Location', 'bestoutside');
axis square
grid on;

% Add super title
sgtitle('[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz', ...
        'Interpreter', 'tex', 'FontSize', 18);

%% Generate lists

% Set display format
format longG

% Generate the lists
list1 = linspace(1E5, 5E6, 25);
list2 = linspace(80, 105, 25);

% Convert to strings with no scientific notation
str_list1 = compose('%.0f', list1);
str_list2 = compose('%.2f', list2);

% Join as space-separated strings
row1 = strjoin(str_list1', ' ');
row2 = strjoin(str_list2', ' ');

% Display results
disp(row1)
disp(row2)

%% Use space-separated floating-point/integer values
SCATTERING_LENGTH_RANGE = "[80.00 81.04 82.08 83.12 84.17 85.21 86.25 87.29 88.33 89.38 90.42 91.46 92.50 93.54 94.58 95.62 96.67 97.71 98.75 99.79 100.83 101.88 102.92 103.96 105.00]";
NUM_ATOMS_LIST          = "[100000 304167 508333 712500 916667 1120833 1325000 1529167 1733333 1937500 2141667 2345833 2550000 2754167 2958333 3162500 3366667 3570833 3775000 3979167 4183333 4387500 4591667 4795833 5000000]";

SCATTERING_LENGTH_RANGE ="[75.00 76.09 77.18 78.27 79.36 80.00 81.04 82.08 83.12 84.17 85.21 86.25 87.29]";
NUM_ATOMS_LIST          ="[50000 54545 59091 63636 68182 72727 77273 81818 86364 90909 95455]";

%% Explore the states
SCATTERING_LENGTH_RANGE = [95.62];
NUM_ATOMS_LIST          = [100000 304167 508333 712500 916667 1120833 1325000 1529167 1733333 1937500 2141667 2345833 2550000 2754167 2958333 3162500 3366667 3570833 3775000 3979167 4183333 4387500 4591667 4795833 5000000];
% NUM_ATOMS_LIST          = [4183333 4387500 4591667 4795833 5000000];

JobNumber = 0;

for i = 1:length(SCATTERING_LENGTH_RANGE)
    aS = SCATTERING_LENGTH_RANGE(i);
    for j = 1:length(NUM_ATOMS_LIST)
        N = NUM_ATOMS_LIST(j);
        
        SaveDirectory = sprintf('D:/Results - Numerics/Data_Full3D/PhaseDiagram/ImagTimePropagation/Theta0/HighN/aS_%.6e_theta_000_phi_000_N_%d', aS, N);
        fprintf('Processing JobNumber %d: %s\n', JobNumber, SaveDirectory);

        % Call the plotting function
        try
            Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
        catch ME
            warning(ME.message)
            continue;
        end
    
        % Pause to inspect plot before continuing
        pause(1.5)
    end
end
%% Missing files

data = {
    ''
};

% Preallocate arrays
numEntries = numel(data);
aS_values = zeros(numEntries, 1);
N_values = zeros(numEntries, 1);

% Regular expression pattern
pattern = 'aS_([\d.eE+-]+)_.*?_N_(\d+)';

% Extract values using regular expressions
for i = 1:numEntries
    tokens = regexp(data{i}, pattern, 'tokens');
    if ~isempty(tokens)
        aS_values(i) = str2double(tokens{1}{1});
        N_values(i) = str2double(tokens{1}{2});
    end
end

% Create and display table
T = table(aS_values, N_values, 'VariableNames', {'aS', 'N'});
disp(T);


%% Visualize phase diagram
load('./Results/Data_Full3D/PhaseDiagramTilted_Theta_20.mat')
PhaseDiagramMatrix   = M;
ScatteringLengths    = SCATTERING_LENGTH_RANGE;
NumberOfAtoms        = round(NUM_ATOMS_LIST * 1E-5,2);
xlen                 = length(NumberOfAtoms);
ylen                 = length(ScatteringLengths);

figure(1);
set(gcf,'Position',[100 100 950 750])
fig.WindowState = 'maximized';
clf
set(gcf,'Position', [100, 100, 1600, 900])
set(gca,'FontSize',16,'Box','On','Linewidth',2);
hImg = imagesc(PhaseDiagramMatrix);
set(gca, 'YDir', 'normal');
colormap(parula);
colorbar;
axis equal tight;
xticks(1:xlen);
yticks(1:ylen);
xticklabels(string(NumberOfAtoms));
yticklabels(string(ScatteringLengths));
xlabel('Number of Atoms', 'Interpreter', 'tex', 'FontSize', 16);
ylabel('Scattering Length (\times a_o)', 'Interpreter', 'tex', 'FontSize', 16);
grid on;

%% Edit phase diagram
load('./Results/Data_Full3D/PhaseDiagramTilted_Theta_20.mat')
PhaseDiagramMatrix   = M;
ScatteringLengths    = SCATTERING_LENGTH_RANGE;
NumberOfAtoms        = NUM_ATOMS_LIST;

Scripts.editPhaseDiagram(PhaseDiagramMatrix, ScatteringLengths, NumberOfAtoms)

%% Smoothen phase diagram
load('./Results/Data_Full3D/PhaseDiagramTilted_Theta_20.mat');  % Load M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST
TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 40^\circ";

Scripts.plotSmoothedPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, TitleString);

%% Select boundary points from original phase diagram
load('./Results/Data_Full3D/PhaseDiagramTilted_Theta_20.mat')
M = M;
N_atoms = round(NUM_ATOMS_LIST * 1E-5, 2);
scatt_lengths = SCATTERING_LENGTH_RANGE;
savePath = './Results/Data_Full3D//BoundaryPoints/SelectedPoints_Theta20_1.mat';

Scripts.selectPhaseDiagramPoints(M, N_atoms, scatt_lengths, savePath);

%% Interactively modify selected boundary points
load('./Results/Data_Full3D/PhaseDiagramTilted_Theta_20.mat')
TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 20^\circ";

Scripts.plotSmoothedContourOnPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, ...
    TitleString, ...
    './Results/Data_Full3D/BoundaryPoints/SelectedPoints_Theta20_1', './Results/Data_Full3D/BoundaryPoints/SelectedPoints_Theta20_2', ...
    'interpMethod', 'pchip', 'showPoints', true);

%% Edit modified points
load('./Results/Data_Full3D/PhaseDiagramTilted_Theta_20.mat')
TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 20^\circ";

Scripts.plotSmoothedContourOnPhaseDiagram(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, ...
    TitleString, ...
    './Results/Data_Full3D/BoundaryPoints/ModifiedPoints_Theta20_1', './Results/Data_Full3D/BoundaryPoints/ModifiedPoints_Theta20_2',  ...
    'interpMethod', 'pchip', 'showPoints', true);

%% Plot with interpolated phase boundary
load('./Results/Data_Full3D/PhaseDiagramTilted_Theta_40.mat')
PhaseDiagramMatrix   = M;
ScatteringLengths    = SCATTERING_LENGTH_RANGE;
NumberOfAtoms        = NUM_ATOMS_LIST;
PhaseBoundary        = load("./Results/Data_Full3D/BoundaryPoints/PhaseBoundary_Tilted_Theta40.mat");
TitleString = "[ \omega_x, \omega_y, \omega_z ] = 2 \pi \times [ 50, 20, 150 ] Hz; \theta = 40^\circ";

Scripts.plotPhaseDiagramWithBoundaries(M, SCATTERING_LENGTH_RANGE, NUM_ATOMS_LIST, PhaseBoundary, TitleString);

%% Density modulation determination

% SaveDirectory = './Results/Data_Full3D/GradientDescent/Phi020/aS_090_theta_020_phi_000_N_100000';
% JobNumber     = 0;

SaveDirectory = './Results/Data_Full3D/PhaseDiagram';
JobNumber     = 1;

% Load data
Data          = load(sprintf(horzcat(SaveDirectory, '/Run_%03i/psi_gs.mat'),JobNumber),'psi','Params','Transf','Observ');

Params        = Data.Params;
Transf        = Data.Transf; 
Observ        = Data.Observ; 

if isgpuarray(Data.psi)
    psi = gather(Data.psi);
else
    psi = Data.psi;
end
if isgpuarray(Data.Observ.residual)
    Observ.residual = gather(Data.Observ.residual); 
else
    Observ.residual = Data.Observ.residual; 
end

ModulationFlag = determineDensityModulation(psi, Params, Transf);

if ModulationFlag
    disp('The state is modulated');
else
    disp('The state is not modulated');
end