%% 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           = 1E5;
OptionsStruct.DipolarPolarAngle       = 0;
OptionsStruct.DipolarAzimuthAngle     = 0;
OptionsStruct.ScatteringLength        = 86;

OptionsStruct.TrapFrequencies         = [10, 10, 72.4];
OptionsStruct.TrapDepth               = 5;
OptionsStruct.BoxSize                 = 15;
OptionsStruct.TrapPotentialType       = 'Harmonic';

OptionsStruct.NumberOfGridPoints      = [256, 512, 256];
OptionsStruct.Dimensions              = [50, 120, 150];
OptionsStruct.CutoffType              = 'Cylindrical';
OptionsStruct.SimulationMode          = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution'
OptionsStruct.TimeStepSize            = 500E-6;                   % in s
OptionsStruct.NumberOfTimeSteps       = 2E6;                      % in s
OptionsStruct.EnergyTolerance         = 5E-10;
OptionsStruct.ResidualTolerance       = 1E-05;

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

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

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

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

%% - Create Variational2D and Calculator object with specified options

OptionsStruct = struct;

OptionsStruct.NumberOfAtoms           = 1E5;     
OptionsStruct.DipolarPolarAngle       = 0;
OptionsStruct.DipolarAzimuthAngle     = 0;
OptionsStruct.ScatteringLength        = 98.0676;     % Critical point: 102.515; Triangular phase: 98.0676; Stripe phase: 102.2518; Honeycomb phase: 102.6441

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

OptionsStruct.NumberOfGridPoints      = [128, 128];
OptionsStruct.Dimensions              = [7.5, 7.5];  % Critical point: 6.996; Triangular phase: 7.5; Stripe phase: 6.972; Honeycomb phase: 6.239 for both for Atom Number fixed to 1E5
OptionsStruct.TimeStepSize            = 5E-6;        % in s
OptionsStruct.TimeCutOff              = 1E5;         % in s
OptionsStruct.EnergyTolerance         = 5E-10;
OptionsStruct.ResidualTolerance       = 1E-05;

OptionsStruct.MaxIterations           = 20;      
OptionsStruct.VariationalWidth        = 4;       
OptionsStruct.WidthLowerBound         = 0.2;
OptionsStruct.WidthUpperBound         = 12;
OptionsStruct.WidthCutoff             = 1e-4;

OptionsStruct.PlotLive                = true;
OptionsStruct.JobNumber               = 1;
OptionsStruct.RunOnGPU                = false;
OptionsStruct.SaveData                = true;
OptionsStruct.SaveDirectory           = './Data';
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 initial wavefunction
Plotter.visualizeWavefunction2D(psi,Params,Transf)
%% - Plot GS wavefunction
Plotter.visualizeGSWavefunction2D(Params.njob)