%%

% To reproduce results from the Blair Blakie paper:

% 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  = ((nadd^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         = [10, 10, 72.4];
OptionsStruct.TrapPotentialType       = 'None';

OptionsStruct.NumberOfGridPoints      = [256, 256];
OptionsStruct.Dimensions              = [100, 100];      
OptionsStruct.TimeStepSize            = 100E-6;          % in s
OptionsStruct.MinimumTimeStepSize     = 1E-5;            % in s
OptionsStruct.TimeCutOff              = 2E6;             % in s
OptionsStruct.EnergyTolerance         = 5E-10;
OptionsStruct.ResidualTolerance       = 1E-04;
OptionsStruct.NoiseScaleFactor        = 4;

OptionsStruct.MaxIterations           = 20;      
OptionsStruct.VariationalWidth        = 5.7;       
OptionsStruct.WidthLowerBound         = 0.2;
OptionsStruct.WidthUpperBound         = 20;
OptionsStruct.WidthCutoff             = 1e-2;

OptionsStruct.PlotLive                = false;
OptionsStruct.JobNumber               = 1;
OptionsStruct.RunOnGPU                = true;
OptionsStruct.SaveData                = true;
OptionsStruct.SaveDirectory           = './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();

%% - 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         = [10, 10, 72.4];
OptionsStruct.TrapPotentialType       = 'None';

OptionsStruct.NumberOfGridPoints      = [256, 256];
OptionsStruct.Dimensions              = [100, 100];      
OptionsStruct.TimeStepSize            = 100E-6;          % in s
OptionsStruct.MinimumTimeStepSize     = 1E-5;            % in s
OptionsStruct.TimeCutOff              = 2E6;             % in s
OptionsStruct.EnergyTolerance         = 5E-10;
OptionsStruct.ResidualTolerance       = 1E-04;
OptionsStruct.NoiseScaleFactor        = 4;

OptionsStruct.MaxIterations           = 20;      
OptionsStruct.VariationalWidth        = 5.7;       
OptionsStruct.WidthLowerBound         = 0.2;
OptionsStruct.WidthUpperBound         = 20;
OptionsStruct.WidthCutoff             = 1e-2;

OptionsStruct.PlotLive                = false;
OptionsStruct.JobNumber               = 1;
OptionsStruct.RunOnGPU                = true;
OptionsStruct.SaveData                = true;
OptionsStruct.SaveDirectory           = './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.0102e+07;     
OptionsStruct.DipolarPolarAngle       = 0;
OptionsStruct.DipolarAzimuthAngle     = 0;
OptionsStruct.ScatteringLength        = 101.9903;        

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

OptionsStruct.NumberOfGridPoints      = [256, 256];
OptionsStruct.Dimensions              = [100, 100];      
OptionsStruct.TimeStepSize            = 100E-6;          % in s
OptionsStruct.MinimumTimeStepSize     = 1E-5;            % in s
OptionsStruct.TimeCutOff              = 2E6;             % in s
OptionsStruct.EnergyTolerance         = 5E-10;
OptionsStruct.ResidualTolerance       = 1E-04;
OptionsStruct.NoiseScaleFactor        = 4;

OptionsStruct.MaxIterations           = 20;      
OptionsStruct.VariationalWidth        = 5.7;       
OptionsStruct.WidthLowerBound         = 0.2;
OptionsStruct.WidthUpperBound         = 20;
OptionsStruct.WidthCutoff             = 1e-2;

OptionsStruct.PlotLive                = false;
OptionsStruct.JobNumber               = 1;
OptionsStruct.RunOnGPU                = true;
OptionsStruct.SaveData                = true;
OptionsStruct.SaveDirectory           = './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();