Latest working version
This commit is contained in:
parent
cfe15aa13d
commit
570995f3f5
@ -17,8 +17,9 @@ OptionsStruct.CutoffType = 'Cylindrical';
|
|||||||
OptionsStruct.TrapPotentialType = 'Harmonic';
|
OptionsStruct.TrapPotentialType = 'Harmonic';
|
||||||
|
|
||||||
OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution'
|
OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution'
|
||||||
OptionsStruct.TimeStep = 50e-06; % in s
|
OptionsStruct.TimeStep = 50E-6; % in s
|
||||||
OptionsStruct.SimulationTime = 4e-03; % in s
|
OptionsStruct.SimulationTime = 2E6; % in s
|
||||||
|
OptionsStruct.EnergyTolerance = 5E-10;
|
||||||
|
|
||||||
OptionsStruct.SaveData = true;
|
OptionsStruct.SaveData = true;
|
||||||
OptionsStruct.SaveDirectory = './Data';
|
OptionsStruct.SaveDirectory = './Data';
|
||||||
|
@ -8,10 +8,12 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
|
|||||||
TrapFrequencies;
|
TrapFrequencies;
|
||||||
NumberOfGridPoints;
|
NumberOfGridPoints;
|
||||||
Dimensions;
|
Dimensions;
|
||||||
|
|
||||||
SimulationMode;
|
SimulationMode;
|
||||||
TimeStep;
|
TimeStep;
|
||||||
SimulationTime;
|
SimulationTime;
|
||||||
|
EnergyTolerance;
|
||||||
|
MinimumTimeStep;
|
||||||
|
|
||||||
%Flags
|
%Flags
|
||||||
|
|
||||||
@ -46,8 +48,12 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
|
|||||||
@(x) any(strcmpi(x,{'ImaginaryTimeEvolution','RealTimeEvolution'})));
|
@(x) any(strcmpi(x,{'ImaginaryTimeEvolution','RealTimeEvolution'})));
|
||||||
addParameter(p, 'TimeStep', 5E-4,...
|
addParameter(p, 'TimeStep', 5E-4,...
|
||||||
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
|
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
|
||||||
addParameter(p, 'SimulationTime', 3,...
|
addParameter(p, 'SimulationTime', 2e6,...
|
||||||
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
|
@(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,...
|
||||||
|
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
|
||||||
addParameter(p, 'DebugMode', false,...
|
addParameter(p, 'DebugMode', false,...
|
||||||
@islogical);
|
@islogical);
|
||||||
addParameter(p, 'SaveData', false,...
|
addParameter(p, 'SaveData', false,...
|
||||||
@ -68,6 +74,8 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
|
|||||||
this.SimulationMode = p.Results.SimulationMode;
|
this.SimulationMode = p.Results.SimulationMode;
|
||||||
this.TimeStep = p.Results.TimeStep;
|
this.TimeStep = p.Results.TimeStep;
|
||||||
this.SimulationTime = p.Results.SimulationTime;
|
this.SimulationTime = p.Results.SimulationTime;
|
||||||
|
this.EnergyTolerance = p.Results.EnergyTolerance;
|
||||||
|
this.MinimumTimeStep = p.Results.MinimumTimeStep;
|
||||||
|
|
||||||
this.DebugMode = p.Results.DebugMode;
|
this.DebugMode = p.Results.DebugMode;
|
||||||
this.DoSave = p.Results.SaveData;
|
this.DoSave = p.Results.SaveData;
|
||||||
|
@ -1,34 +1,34 @@
|
|||||||
function [Params] = setupParameters(this)
|
function [Params] = setupParameters(this)
|
||||||
|
|
||||||
CONSTANTS = Helper.PhysicsConstants;
|
CONSTANTS = Helper.PhysicsConstants;
|
||||||
hbar = CONSTANTS.PlanckConstantReduced; % [J.s]
|
hbar = CONSTANTS.PlanckConstantReduced; % [J.s]
|
||||||
kbol = CONSTANTS.BoltzmannConstant; % [J/K]
|
kbol = CONSTANTS.BoltzmannConstant; % [J/K]
|
||||||
mu0 = CONSTANTS.VacuumPermeability; % [N/A^2]
|
mu0 = CONSTANTS.VacuumPermeability; % [N/A^2]
|
||||||
muB = CONSTANTS.BohrMagneton; % [J/T]
|
muB = CONSTANTS.BohrMagneton; % [J/T]
|
||||||
a0 = CONSTANTS.BohrRadius; % [m]
|
a0 = CONSTANTS.BohrRadius; % [m]
|
||||||
m0 = CONSTANTS.AtomicMassUnit; % [kg]
|
m0 = CONSTANTS.AtomicMassUnit; % [kg]
|
||||||
w0 = 2*pi*100; % Angular frequency unit [s^-1]
|
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.
|
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
|
% Number of points in each direction
|
||||||
Params.Nx = this.NumberOfGridPoints(1);
|
Params.Nx = this.NumberOfGridPoints(1);
|
||||||
Params.Ny = this.NumberOfGridPoints(2);
|
Params.Ny = this.NumberOfGridPoints(2);
|
||||||
Params.Nz = this.NumberOfGridPoints(3);
|
Params.Nz = this.NumberOfGridPoints(3);
|
||||||
|
|
||||||
% Dimensions (in units of l0)
|
% Dimensions (in units of l0)
|
||||||
Params.Lx = this.Dimensions(1);
|
Params.Lx = this.Dimensions(1);
|
||||||
Params.Ly = this.Dimensions(2);
|
Params.Ly = this.Dimensions(2);
|
||||||
Params.Lz = this.Dimensions(3);
|
Params.Lz = this.Dimensions(3);
|
||||||
|
|
||||||
% Masses
|
% Masses
|
||||||
Params.m = CONSTANTS.Dy164Mass;
|
Params.m = CONSTANTS.Dy164Mass;
|
||||||
l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length
|
l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length
|
||||||
|
|
||||||
% Atom numbers
|
% Atom numbers
|
||||||
Params.N = this.NumberOfAtoms;
|
Params.N = this.NumberOfAtoms;
|
||||||
|
|
||||||
% Dipole angle
|
% Dipole angle
|
||||||
Params.theta = this.DipolarPolarAngle; % pi/2 dipoles along x, theta=0 dipoles along z
|
Params.theta = this.DipolarPolarAngle; % pi/2 dipoles along x, theta=0 dipoles along z
|
||||||
Params.phi = this.DipolarAzimuthAngle;
|
Params.phi = this.DipolarAzimuthAngle;
|
||||||
|
|
||||||
% Dipole lengths (units of muB)
|
% Dipole lengths (units of muB)
|
||||||
@ -38,35 +38,34 @@ Params.mu = CONSTANTS.DyMagneticMoment;
|
|||||||
Params.as = this.ScatteringLength*a0;
|
Params.as = this.ScatteringLength*a0;
|
||||||
|
|
||||||
% Trapping frequencies
|
% Trapping frequencies
|
||||||
Params.wx = 2*pi*this.TrapFrequencies(1);
|
Params.wx = 2*pi*this.TrapFrequencies(1);
|
||||||
Params.wy = 2*pi*this.TrapFrequencies(2);
|
Params.wy = 2*pi*this.TrapFrequencies(2);
|
||||||
Params.wz = 2*pi*this.TrapFrequencies(3);
|
Params.wz = 2*pi*this.TrapFrequencies(3);
|
||||||
|
|
||||||
% Stochastic GPE
|
% Stochastic GPE
|
||||||
Params.gamma_S = 7.5*10^(-3); % gamma for the stochastic GPE
|
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.muchem = 12.64*Params.wz/w0; % fixing the chemical potential for the stochastic GPE
|
||||||
|
|
||||||
Params.Etol = 5e-10; % Tolerances
|
|
||||||
Params.rtol = 1e-5;
|
|
||||||
Params.cut_off = 2e6; % sometimes the imaginary time gets a little stuck
|
|
||||||
% even though the solution is good, this just stops it going on forever
|
|
||||||
Params.mindt = 1e-6; % Minimum size for a time step using adaptive dt
|
|
||||||
|
|
||||||
|
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
|
||||||
|
|
||||||
% ================ Parameters defined by those above ================ %
|
% ================ Parameters defined by those above ================ %
|
||||||
|
|
||||||
% Contact interaction strength (units of l0/m)
|
% Contact interaction strength (units of l0/m)
|
||||||
Params.gs = 4*pi*Params.as/l0;
|
Params.gs = 4*pi*Params.as/l0;
|
||||||
|
|
||||||
% Dipole lengths
|
% Dipole lengths
|
||||||
Params.add = mu0*Params.mu^2*Params.m/(12*pi*hbar^2);
|
Params.add = mu0*Params.mu^2*Params.m/(12*pi*hbar^2);
|
||||||
|
|
||||||
% DDI strength
|
% DDI strength
|
||||||
Params.gdd = 12*pi*Params.add/l0; %sometimes the 12 is a 4 --> depends on how Vdk (DDI) is defined
|
Params.gdd = 12*pi*Params.add/l0; %sometimes the 12 is a 4 --> depends on how Vdk (DDI) is defined
|
||||||
|
|
||||||
% Trap gamma
|
% Trap gamma
|
||||||
Params.gx=(Params.wx/w0)^2;
|
Params.gx = (Params.wx/w0)^2;
|
||||||
Params.gy=(Params.wy/w0)^2;
|
Params.gy = (Params.wy/w0)^2;
|
||||||
Params.gz=(Params.wz/w0)^2;
|
Params.gz = (Params.wz/w0)^2;
|
||||||
|
|
||||||
% == Calculate LHY correction == %
|
% == Calculate LHY correction == %
|
||||||
eps_dd = Params.add/Params.as;
|
eps_dd = Params.add/Params.as;
|
||||||
@ -83,6 +82,11 @@ end
|
|||||||
Params.gammaQF = 128/3*sqrt(pi*(Params.as/l0)^5)*Q5;
|
Params.gammaQF = 128/3*sqrt(pi*(Params.as/l0)^5)*Q5;
|
||||||
|
|
||||||
% Loading the rest into Params
|
% Loading the rest into Params
|
||||||
Params.hbar = hbar; Params.kbol = kbol; Params.mu0 = mu0; Params.muB = muB; Params.a0 = a0;
|
Params.hbar = hbar;
|
||||||
Params.w0 = w0; Params.l0 = l0;
|
Params.kbol = kbol;
|
||||||
|
Params.mu0 = mu0;
|
||||||
|
Params.muB = muB;
|
||||||
|
Params.a0 = a0;
|
||||||
|
Params.w0 = w0;
|
||||||
|
Params.l0 = l0;
|
||||||
end
|
end
|
Loading…
Reference in New Issue
Block a user