karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

Expanded the Potential class.

Browse Source
This commit is contained in:
Karthik 2024-06-17 20:24:00 +02:00
parent 6bb354de7a
commit e67f82b564
11 changed files with 385 additions and 163 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
Dipolar Gas Simulator/+Scripts/test.m
View File

@ -10,28 +10,31 @@ OptionsStruct.NumberOfAtoms = 1E6;
OptionsStruct.DipolarPolarAngle = pi/2;
OptionsStruct.DipolarAzimuthAngle = 0;
OptionsStruct.ScatteringLength = 86;
OptionsStruct.TrapFrequencies = [125, 125, 250];
OptionsStruct.TrapDepth = 5;
OptionsStruct.BoxSize = 15;
OptionsStruct.TrapPotentialType = 'Harmonic';
OptionsStruct.NumberOfGridPoints = [64, 64, 48];
OptionsStruct.Dimensions = [40, 40, 20];
OptionsStruct.CutoffType = 'Cylindrical';
OptionsStruct.TrapPotentialType = 'Harmonic';
OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution'
OptionsStruct.TimeStep = 50E-6; % in s
OptionsStruct.SimulationTime = 2E6; % in s
OptionsStruct.TimeStepSize = 50E-6; % in s
OptionsStruct.NumberOfTimeSteps = 2E6; % in s
OptionsStruct.EnergyTolerance = 5E-10;
OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = './Data';
options = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct
sim = Simulator.DipolarGas(options{:});
pot = Simulator.Potentials(options{:});
sim = Simulator.DipolarGas(options{:});
pot = Simulator.Potentials(options{:});
sim.Potential = pot.trap(); % + pot.repulsive_stirrer();
%-% Run Simulation %-%
[Params, Transf, psi, V, VDk] = sim.runSimulation();
[Params, Transf, psi, V, VDk] = sim.run();
%% - Plot numerical grid
Plotter.visualizeSpace(Transf)

4
Dipolar Gas Simulator/+Simulator/@Calculator/Calculator.m
View File

@ -42,8 +42,6 @@ classdef Calculator < handle & matlab.mixin.Copyable
this.PhaseCoherence = this.CalculatorDefaults.PhaseCoherence;
this.TotalEnergy = this.CalculatorDefaults.TotalEnergy;
this.CutoffType = p.Results.CutoffType;
this.CalculatorDefaults.CutoffType = this.CutoffType;
end
function restoreDefaults(this)
@ -56,7 +54,7 @@ classdef Calculator < handle & matlab.mixin.Copyable
this.CutoffType = this.CalculatorDefaults.CutoffType;
end
end % - lifecycle
end %
% methods

57
Dipolar Gas Simulator/+Simulator/@DipolarGas/DipolarGas.m
View File

@ -8,14 +8,17 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
TrapFrequencies;
NumberOfGridPoints;
Dimensions;
Potential;
SimulationMode;
TimeStep;
SimulationTime;
TimeStepSize;
NumberOfTimeSteps;
EnergyTolerance;
MinimumTimeStep;
MinimumTimeStepSize;
Calculator;
SimulationParameters;
%Flags
@ -47,16 +50,16 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
addParameter(p, 'Dimensions', 10 * ones(1,3),...
@(x) assert(isnumeric(x) && isvector(x) && all(x > 0)));
addParameter(p, 'SimulationMode', 'ImaginaryTimeEvolution',...
@(x) any(strcmpi(x,{'ImaginaryTimeEvolution','RealTimeEvolution'})));
@(x) assert(any(strcmpi(x,{'ImaginaryTimeEvolution','RealTimeEvolution'}))));
addParameter(p, 'CutoffType', 'Cylindrical',...
@(x) any(strcmpi(x,{'Cylindrical','CylindricalInfiniteZ', 'Spherical'})));
addParameter(p, 'TimeStep', 5E-4,...
@(x) assert(any(strcmpi(x,{'Cylindrical','CylindricalInfiniteZ', 'Spherical'}))));
addParameter(p, 'TimeStepSize', 5E-4,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'SimulationTime', 2e6,...
addParameter(p, 'NumberOfTimeSteps', 2e6,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'EnergyTolerance', 1e-10,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'MinimumTimeStep', 1e-6,...
addParameter(p, 'MinimumTimeStepSize', 1e-6,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'DebugMode', false,...
@islogical);
@ -74,12 +77,12 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
this.TrapFrequencies = p.Results.TrapFrequencies;
this.NumberOfGridPoints = p.Results.NumberOfGridPoints;
this.Dimensions = p.Results.Dimensions;
this.Potential = NaN;
this.SimulationMode = p.Results.SimulationMode;
this.TimeStep = p.Results.TimeStep;
this.SimulationTime = p.Results.SimulationTime;
this.TimeStepSize = p.Results.TimeStepSize;
this.NumberOfTimeSteps = p.Results.NumberOfTimeSteps;
this.EnergyTolerance = p.Results.EnergyTolerance;
this.MinimumTimeStep = p.Results.MinimumTimeStep;
this.MinimumTimeStepSize = p.Results.MinimumTimeStepSize;
this.DebugMode = p.Results.DebugMode;
this.DoSave = p.Results.SaveData;
@ -87,29 +90,25 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
this.Calculator = Simulator.Calculator('CutoffType', p.Results.CutoffType);
switch this.SimulationMode
case "ImaginaryTimeEvolution"
% Development In progress
case "RealTimeEvolution"
% Development In progress
end
this.SimulationParameters = this.setupParameters();
end
end
methods
function set.TimeStep(this, val)
function set.TimeStepSize(this, val)
assert(val > 1e-06, 'Not time efficient to compute for time steps smaller than 1 microsecond!');
this.TimeStep = val;
this.TimeStepSize = val;
end
function ret = get.TimeStep(this)
ret = this.TimeStep;
function ret = get.TimeStepSize(this)
ret = this.TimeStepSize;
end
function set.SimulationTime(this, val)
% assert(val <= 5e-03, 'Not time efficient to compute for time spans longer than 5 milliseconds!');
this.SimulationTime = val;
function set.NumberOfTimeSteps(this, val)
assert(val <= 2E6, 'Not time efficient to compute for time spans longer than 2E6 seconds!');
this.NumberOfTimeSteps = val;
end
function ret = get.SimulationTime(this)
ret = this.SimulationTime;
function ret = get.NumberOfTimeSteps(this)
ret = this.NumberOfTimeSteps;
end
function set.NumberOfAtoms(this, val)
assert(val <= 1E6, '!!Not time efficient to compute for atom numbers larger than 1,000,000!!');
@ -118,7 +117,6 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
function ret = get.NumberOfAtoms(this)
ret = this.NumberOfAtoms;
end
function set.DebugMode(this, val)
this.DebugMode = val;
end
@ -137,7 +135,6 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
function ret = get.SaveDirectory(this)
ret = this.SaveDirectory;
end
end % - setters and getters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

9
Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m
View File

@ -1,9 +1,8 @@
function [psi,V,VDk] = initialize(this,Params,Transf,TransfRad)
format long
X = Transf.X; Y = Transf.Y; Z = Transf.Z;
% == Trap potential == %
V = 0.5*(Params.gx.*X.^2+Params.gy.*Y.^2+Params.gz*Z.^2);
% == User-defined potential == %
V = this.Potential;
assert(~anynan(V), 'Potential not defined! Specify as <SimulatorObject>.Potential = <PotentialsObject>.trap() + <AdditionalTerms>.');
% == Calculating the DDIs == %
if isfile(strcat(this.SaveDirectory, '/VDk_M.mat'))
@ -13,7 +12,7 @@ function [psi,V,VDk] = initialize(this,Params,Transf,TransfRad)
VDk = this.Calculator.calculateVDCutoff(Params,Transf,TransfRad);
save(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat')),'VDk');
end
fprintf('Computed and saved DDI potential in Fourier space with %s cutoff.', this.Calculator.CutoffType)
fprintf('Computed and saved DDI potential in Fourier space with %s cutoff.\n', this.Calculator.CutoffType)
% == Setting up the initial wavefunction == %
psi = this.setupWavefunction(Params,Transf);

170
Dipolar Gas Simulator/+Simulator/@DipolarGas/propagateWavefunction.m Normal file
View File

@ -0,0 +1,170 @@
function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,njob,t_idx,Observ)
set(0,'defaulttextInterpreter','latex')
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
switch this.SimulationMode
case 'ImaginaryTimeEvolution'
dt=-1j*abs(this.TimeStepSize);
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
Observ.residual = 1; Observ.res = 1;
muchem = this.Calculator.ChemicalPotential(psi,Params,Transf,VDk,V);
AdaptIdx = 0;
while t_idx < Params.sim_time_cut_off
%kin
psi = fftn(psi);
psi = psi.*exp(-0.5*1i*dt*KEop);
psi = ifftn(psi);
%DDI
frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk));
%Real-space
psi = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gdd*Phi + Params.gammaQF*abs(psi).^3 - muchem));
%kin
psi = fftn(psi);
psi = psi.*exp(-0.5*1i*dt*KEop);
psi = ifftn(psi);
%Renorm
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
psi = sqrt(Params.N)*psi/sqrt(Norm);
muchem = this.Calculator.ChemicalPotential(psi,Params,Transf,VDk,V);
if mod(t_idx,1000) == 0
%Change in Energy
E = this.Calculator.TotalEnergy(psi,Params,Transf,VDk,V);
E = E/Norm;
Observ.EVec = [Observ.EVec E];
%Chemical potential
Observ.mucVec = [Observ.mucVec muchem];
%Normalized residuals
res = this.Calculator.NormalizedResiduals(psi,Params,Transf,VDk,V,muchem);
Observ.residual = [Observ.residual res];
Observ.res_idx = Observ.res_idx + 1;
save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
%Adaptive time step -- Careful, this can quickly get out of control
relres = abs(Observ.residual(Observ.res_idx)-Observ.residual(Observ.res_idx-1))/Observ.residual(Observ.res_idx);
if relres <1e-5
if AdaptIdx > 4 && abs(dt) > Params.mindt
dt = dt / 2;
fprintf('Time step changed to '); disp(dt);
AdaptIdx = 0;
elseif AdaptIdx > 4 && abs(dt) < Params.mindt
break
else
AdaptIdx = AdaptIdx + 1;
end
else
AdaptIdx = 0;
end
end
if any(isnan(psi(:)))
disp('NaNs encountered!')
break
end
t_idx=t_idx+1;
end
%Change in Energy
E = this.Calculator.TotalEnergy(psi,Params,Transf,VDk,V);
E = E/Norm;
Observ.EVec = [Observ.EVec E];
% Phase coherence
[PhaseC] = this.Calculator.PhaseCoherence(psi,Transf,Params);
Observ.PCVec = [Observ.PCVec PhaseC];
Observ.res_idx = Observ.res_idx + 1;
save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
case 'RealTimeEvolution'
dt = abs(this.TimeStepSize);
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
muchem = chemicalpotential(psi,Params,Transf,VDk,V);
while t_idx < Params.sim_time_cut_off
% Parameters at time t
%kin
psi = fftn(psi);
psi = psi.*exp(-0.5*1i*dt*KEop);
psi = ifftn(psi);
%DDI
frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk));
%Real-space
psi = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gammaQF*abs(psi).^3 + Params.gdd*Phi - muchem));
%kin
psi = fftn(psi);
psi = psi.*exp(-0.5*1i*dt*KEop);
psi = ifftn(psi);
muchem = chemicalpotential(psi,Params,Transf,VDk,V);
if mod(t_idx,1000)==0
%Change in Normalization
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; % normalisation
Observ.NormVec = [Observ.NormVec Norm];
%Change in Energy
E = energytotal(psi,Params,Transf,VDk,V);
E = E/Norm;
Observ.EVec = [Observ.EVec E];
% Phase coherence
[PhaseC] = PhaseCoherence(psi,Transf);
Observ.PCVec = [Observ.PCVec PhaseC];
Observ.tVecPlot = [Observ.tVecPlot tVal];
Observ.res_idx = Observ.res_idx + 1;
save(sprintf('./Data/Run_%03i/TimeEvolution/psi_%i.mat',njob,Observ.res_idx),'psi','muchem','Observ','t_idx');
end
if any(isnan(psi(:)))
disp('NaNs encountered!')
break
end
t_idx=t_idx+1;
end
%Change in Normalization
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; % normalisation
Observ.NormVec = [Observ.NormVec Norm];
%Change in Energy
E = energytotal(psi,Params,Transf,VDk,V);
E = E/Norm;
Observ.EVec = [Observ.EVec E];
% Phase coherence
[PhaseC] = PhaseCoherence(psi,Transf);
Observ.PCVec = [Observ.PCVec PhaseC];
Observ.tVecPlot = [Observ.tVecPlot tVal];
Observ.res_idx = Observ.res_idx + 1;
save(sprintf('./Data/Run_%03i/TimeEvolution/psi_%i.mat',njob,Observ.res_idx),'psi','muchem','Observ','t_idx');
otherwise
disp('Choose a valid DDI cutoff type!')
return
end
end

12
Dipolar Gas Simulator/+Simulator/@DipolarGas/runSimulation.m → Dipolar Gas Simulator/+Simulator/@DipolarGas/run.m
View File

@ -1,7 +1,8 @@
function [Params, Transf, psi,V,VDk] = runSimulation(this)
function [Params, Transf, psi,V,VDk] = run(this)
% --- Obtain simulation parameters ---
[Params] = this.setupParameters();
this.SimulationParameters = Params;
% --- Set up spatial grids and transforms ---
[Transf] = this.setupSpace(Params);
[TransfRad] = this.setupSpaceRadial(Params);
@ -15,10 +16,7 @@ function [Params, Transf, psi,V,VDk] = runSimulation(this)
t_idx = 1; %Start at t = 0;
Observ.res_idx = 1;
% --- Job Settings ---
njob = 6;
mkdir(sprintf('./Data/Run_%03i',njob))
% --- Run Simulation ---
% [psi] = this.solver(psi,Params,Transf,VDk,V,njob,t_idx,Observ);
% mkdir(sprintf('./Data/Run_%03i',Params.njob))
% [psi] = this.propagateWavefunction(psi,Params,Transf,VDk,V,njob,t_idx,Observ);
end

14
Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m
View File

@ -8,7 +8,7 @@ function [Params] = setupParameters(this)
m0 = CONSTANTS.AtomicMassUnit; % [kg]
w0 = 2*pi*100; % Angular frequency unit [s^-1]
mu0factor = 0.3049584233607396; % =(m0/me)*pi*alpha^2 -- me=mass of electron, alpha=fine struct. const.
% mu0=mu0factor *hbar^2*a0/(m0*muB^2)
% mu0=mu0factor *hbar^2*a0/(m0*muB^2)
% Number of points in each direction
Params.Nx = this.NumberOfGridPoints(1);
Params.Ny = this.NumberOfGridPoints(2);
@ -19,7 +19,7 @@ function [Params] = setupParameters(this)
Params.Ly = this.Dimensions(2);
Params.Lz = this.Dimensions(3);
% Masses
% Mass, length scale
Params.m = CONSTANTS.Dy164Mass;
l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length
@ -45,10 +45,12 @@ function [Params] = setupParameters(this)
Params.gamma_S = 7.5*10^(-3); % gamma for the stochastic GPE
Params.muchem = 12.64*Params.wz/w0; % fixing the chemical potential for the stochastic GPE
Params.Etol = this.EnergyTolerance; % Tolerances
Params.cut_off = this.SimulationTime; % sometimes the imaginary time gets a little stuck
% even though the solution is good, this just stops it going on forever
Params.mindt = this.MinimumTimeStep; % Minimum size for a time step using adaptive dt
Params.Etol = this.EnergyTolerance; % Tolerances
Params.sim_time_cut_off = this.NumberOfTimeSteps; % sometimes the imaginary time gets a little stuck
% even though the solution is good, this just stops it going on forever
Params.mindt = this.MinimumTimeStepSize; % Minimum size for a time step using adaptive dt
Params.njob = 6;
% ================ Parameters defined by those above ================ %

2
Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpace.m
View File

@ -1,4 +1,4 @@
function [Transf] = setupSpace(this,Params)
function [Transf] = setupSpace(~,Params)
Transf.Xmax = 0.5*Params.Lx;
Transf.Ymax = 0.5*Params.Ly;
Transf.Zmax = 0.5*Params.Lz;

89
Dipolar Gas Simulator/+Simulator/@DipolarGas/solver.m
View File

@ -1,89 +0,0 @@
function [psi] = propagate(this,psi,Params,Transf,VDk,V,njob,t_idx,Observ)
set(0,'defaulttextInterpreter','latex')
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
dt=-1j*abs(this.TimeStep);
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
Observ.residual = 1; Observ.res = 1;
muchem = this.Calculator.ChemicalPotential(psi,Params,Transf,VDk,V);
AdaptIdx = 0;
while t_idx < Params.cut_off
%kin
psi = fftn(psi);
psi = psi.*exp(-0.5*1i*dt*KEop);
psi = ifftn(psi);
%DDI
frho = fftn(abs(psi).^2);
Phi = real(ifftn(frho.*VDk));
%Real-space
psi = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gdd*Phi + Params.gammaQF*abs(psi).^3 - muchem));
%kin
psi = fftn(psi);
psi = psi.*exp(-0.5*1i*dt*KEop);
psi = ifftn(psi);
%Renorm
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
psi = sqrt(Params.N)*psi/sqrt(Norm);
muchem = this.Calculator.ChemicalPotential(psi,Params,Transf,VDk,V);
if mod(t_idx,1000) == 0
%Change in Energy
E = this.Calculator.TotalEnergy(psi,Params,Transf,VDk,V);
E = E/Norm;
Observ.EVec = [Observ.EVec E];
%Chemical potential
Observ.mucVec = [Observ.mucVec muchem];
%Normalized residuals
res = this.Calculator.NormalizedResiduals(psi,Params,Transf,VDk,V,muchem);
Observ.residual = [Observ.residual res];
Observ.res_idx = Observ.res_idx + 1;
save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
%Adaptive time step -- Careful, this can quickly get out of control
relres = abs(Observ.residual(Observ.res_idx)-Observ.residual(Observ.res_idx-1))/Observ.residual(Observ.res_idx);
if relres <1e-5
if AdaptIdx > 4 && abs(dt) > Params.mindt
dt = dt / 2;
fprintf('Time step changed to '); disp(dt);
AdaptIdx = 0;
elseif AdaptIdx > 4 && abs(dt) < Params.mindt
break
else
AdaptIdx = AdaptIdx + 1;
end
else
AdaptIdx = 0;
end
end
if any(isnan(psi(:)))
disp('NaNs encountered!')
break
end
t_idx=t_idx+1;
end
%Change in Energy
E = this.Calculator.TotalEnergy(psi,Params,Transf,VDk,V);
E = E/Norm;
Observ.EVec = [Observ.EVec E];
% Phase coherence
[PhaseC] = this.Calculator.PhaseCoherence(psi,Transf,Params);
Observ.PCVec = [Observ.PCVec PhaseC];
Observ.res_idx = Observ.res_idx + 1;
save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
end

136
Dipolar Gas Simulator/+Simulator/@Potentials/Potentials.m
View File

@ -3,12 +3,19 @@ classdef Potentials < handle & matlab.mixin.Copyable
properties (Access = private)
PotentialDefaults = struct('TrapPotentialType', 'Harmonic', ...
'TrapFrequency', 100e3);
'TrapFrequencies', 100 * ones(1,3), ...
'TrapDepth', 10, ...
'BoxSize', 5);
end
properties (Access = public)
TrapPotentialType;
TrapFrequency;
TrapFrequencies;
TrapDepth;
BoxSize;
NumberOfGridPoints;
Dimensions;
SimulationParameters;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -16,18 +23,107 @@ classdef Potentials < handle & matlab.mixin.Copyable
methods
function this = Potentials(varargin)
this.TrapFrequency = this.PotentialDefaults.TrapFrequency;
end
function restoreDefaults(this)
this.TrapFrequency = this.PotentialDefaults.TrapFrequency;
end
end % - lifecycle
p = inputParser;
p.KeepUnmatched = true;
addParameter(p, 'TrapPotentialType', this.PotentialDefaults.TrapPotentialType,...
@(x) assert(any(strcmpi(x,{'None','Harmonic','SquareBox','RoundBox'}))));
addParameter(p, 'TrapFrequencies', this.PotentialDefaults.TrapFrequencies,...
@(x) assert(isnumeric(x) && isvector(x) && all(x > 0)));
addParameter(p, 'TrapDepth', this.PotentialDefaults.TrapDepth,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'BoxSize', this.PotentialDefaults.BoxSize,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'NumberOfGridPoints', 128 * ones(1,3),...
@(x) assert(isnumeric(x) && isvector(x) && all(x > 0)));
addParameter(p, 'Dimensions', 10 * ones(1,3),...
@(x) assert(isnumeric(x) && isvector(x) && all(x > 0)));
p.parse(varargin{:});
% methods
this.TrapPotentialType = p.Results.TrapPotentialType;
this.TrapFrequencies = p.Results.TrapFrequencies;
this.TrapDepth = p.Results.TrapDepth;
this.BoxSize = p.Results.BoxSize;
this.NumberOfGridPoints = p.Results.NumberOfGridPoints;
this.Dimensions = p.Results.Dimensions;
this.SimulationParameters = this.setupParameters();
end
% end % - setters and getters
function [ret] = trap(this)
format long
switch this.TrapPotentialType
case 'None'
x = linspace(-0.5*this.SimulationParameters.Lx,0.5*this.SimulationParameters.Lx-this.SimulationParameters.Lx/this.SimulationParameters.Nx,this.SimulationParameters.Nx);
y = linspace(-0.5*this.SimulationParameters.Ly,0.5*this.SimulationParameters.Ly-this.SimulationParameters.Ly/this.SimulationParameters.Ny,this.SimulationParameters.Ny);
z = linspace(-0.5*this.SimulationParameters.Lz,0.5*this.SimulationParameters.Lz-this.SimulationParameters.Lz/this.SimulationParameters.Nz,this.SimulationParameters.Nz);
[X,Y,Z] = ndgrid(x,y,z);
ret = 0.0 * (X+Y+Z);
case 'Harmonic'
x = linspace(-0.5*this.SimulationParameters.Lx,0.5*this.SimulationParameters.Lx-this.SimulationParameters.Lx/this.SimulationParameters.Nx,this.SimulationParameters.Nx);
y = linspace(-0.5*this.SimulationParameters.Ly,0.5*this.SimulationParameters.Ly-this.SimulationParameters.Ly/this.SimulationParameters.Ny,this.SimulationParameters.Ny);
z = linspace(-0.5*this.SimulationParameters.Lz,0.5*this.SimulationParameters.Lz-this.SimulationParameters.Lz/this.SimulationParameters.Nz,this.SimulationParameters.Nz);
[X,Y,Z] = ndgrid(x,y,z);
ret = 0.5*(this.SimulationParameters.gx.*X.^2+this.SimulationParameters.gy.*Y.^2+this.SimulationParameters.gz*Z.^2);
case 'SquareBox'
x = linspace(-0.5*this.SimulationParameters.Lx,0.5*this.SimulationParameters.Lx-this.SimulationParameters.Lx/this.SimulationParameters.Nx,this.SimulationParameters.Nx);
y = linspace(-0.5*this.SimulationParameters.Ly,0.5*this.SimulationParameters.Ly-this.SimulationParameters.Ly/this.SimulationParameters.Ny,this.SimulationParameters.Ny);
z = linspace(-0.5*this.SimulationParameters.Lz,0.5*this.SimulationParameters.Lz-this.SimulationParameters.Lz/this.SimulationParameters.Nz,this.SimulationParameters.Nz);
[X,Y,Z] = ndgrid(x,y,z);
heaviside_X = this.custom_heaviside(-abs(X) + this.SimulationParameters.boxsize/2);
heaviside_Y = this.custom_heaviside(-abs(Y) + this.SimulationParameters.boxsize/2);
heaviside_Z = this.custom_heaviside(-abs(Z) + this.SimulationParameters.boxsize/2);
ret = this.SimulationParameters.boxdepth * (1 - heaviside_X .* heaviside_Y .* heaviside_Z);
case 'RoundBox'
x = linspace(-0.5*this.SimulationParameters.Lx,0.5*this.SimulationParameters.Lx-this.SimulationParameters.Lx/this.SimulationParameters.Nx,this.SimulationParameters.Nx);
y = linspace(-0.5*this.SimulationParameters.Ly,0.5*this.SimulationParameters.Ly-this.SimulationParameters.Ly/this.SimulationParameters.Ny,this.SimulationParameters.Ny);
z = linspace(-0.5*this.SimulationParameters.Lz,0.5*this.SimulationParameters.Lz-this.SimulationParameters.Lz/this.SimulationParameters.Nz,this.SimulationParameters.Nz);
[X,Y,Z] = ndgrid(x,y,z);
ret = this.SimulationParameters.boxdepth * (1 - this.custom_heaviside((this.SimulationParameters.boxsize/2)^2 - X.^2 - Y.^2 - Z.^2));
end
end
function restoreDefaults(this)
this.TrapPotentialType = this.PotentialDefaults.TrapPotentialType;
this.TrapFrequencies = this.PotentialDefaults.TrapFrequencies;
end
end %
methods
function set.TrapFrequencies(this, val)
assert(isnumeric(val) && isvector(val) && all(val > 0), 'Incorrectly defined trap frequencies!');
this.TrapFrequencies = val;
end
function ret = get.TrapFrequencies(this)
ret = this.TrapFrequencies;
end
function set.TrapPotentialType(this, val)
assert(any(strcmpi(val,{'None','Harmonic','SquareBox','RoundBox'})), 'Trap potential of the specified type cannot be generated!');
this.TrapPotentialType = val;
end
function ret = get.TrapPotentialType(this)
ret = this.TrapPotentialType;
end
function set.NumberOfGridPoints(this, val)
assert(isnumeric(val) && isvector(val) && all(val > 0), 'Incorrectly defined grid!');
this.NumberOfGridPoints = val;
end
function ret = get.NumberOfGridPoints(this)
ret = this.NumberOfGridPoints;
end
function set.Dimensions(this, val)
assert(isnumeric(val) && isvector(val) && all(val > 0), 'Incorrectly defined dimensions!');
this.Dimensions = val;
end
function ret = get.Dimensions(this)
ret = this.Dimensions;
end
end % - setters and getters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- Methods
@ -50,7 +146,19 @@ classdef Potentials < handle & matlab.mixin.Copyable
end
methods (Static)
function y = custom_heaviside(x)
% This function computes the Heaviside step function with a custom value for Heaviside(0).
% x: input array
% H0: value for Heaviside(0)
% Use MATLAB's built-in heaviside function
y = heaviside(x);
% Replace the default value for Heaviside(0) with the custom value H0
y(x == 0) = 1;
end
% Creates an Instance of Class, ensures singleton behaviour (that there
% can only be one Instance of this class
function singleObj = getInstance(varargin)
@ -63,4 +171,4 @@ classdef Potentials < handle & matlab.mixin.Copyable
end
end
end
end

36
Dipolar Gas Simulator/+Simulator/@Potentials/setupParameters.m Normal file
View File

@ -0,0 +1,36 @@
function [Params] = setupParameters(this)
CONSTANTS = Helper.PhysicsConstants;
hbar = CONSTANTS.PlanckConstantReduced; % [J.s]
w0 = 2*pi*100; % Angular frequency unit [s^-1]
% Mass, length scale
Params.m = CONSTANTS.Dy164Mass;
l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length
% Number of points in each direction
Params.Nx = this.NumberOfGridPoints(1);
Params.Ny = this.NumberOfGridPoints(2);
Params.Nz = this.NumberOfGridPoints(3);
% Dimensions (in units of l0)
Params.Lx = this.Dimensions(1);
Params.Ly = this.Dimensions(2);
Params.Lz = this.Dimensions(3);
% Trapping frequencies
Params.wx = 2*pi*this.TrapFrequencies(1);
Params.wy = 2*pi*this.TrapFrequencies(2);
Params.wz = 2*pi*this.TrapFrequencies(3);
% Trap depth and box size for box potentials
Params.boxdepth = this.TrapDepth/(Params.wz/w0); % The depth of the box
Params.boxsize = this.BoxSize; % The size of the box
% ================ Parameters defined by those above ================ %
% Trap gamma
Params.gx = (Params.wx/w0)^2;
Params.gy = (Params.wy/w0)^2;
Params.gz = (Params.wz/w0)^2;
end