diff --git a/Dipolar Gas Simulator/+Scripts/test.m b/Dipolar Gas Simulator/+Scripts/test.m index 1b6d93f..c2679d7 100644 --- a/Dipolar Gas Simulator/+Scripts/test.m +++ b/Dipolar Gas Simulator/+Scripts/test.m @@ -17,8 +17,9 @@ OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.TrapPotentialType = 'Harmonic'; OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' -OptionsStruct.TimeStep = 50e-06; % in s -OptionsStruct.SimulationTime = 4e-03; % in s +OptionsStruct.TimeStep = 50E-6; % in s +OptionsStruct.SimulationTime = 2E6; % in s +OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.SaveData = true; OptionsStruct.SaveDirectory = './Data'; diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/DipolarGas.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/DipolarGas.m index 166ecc2..2005af8 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/DipolarGas.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/DipolarGas.m @@ -8,10 +8,12 @@ classdef DipolarGas < handle & matlab.mixin.Copyable TrapFrequencies; NumberOfGridPoints; Dimensions; - + SimulationMode; TimeStep; SimulationTime; + EnergyTolerance; + MinimumTimeStep; %Flags @@ -46,8 +48,12 @@ classdef DipolarGas < handle & matlab.mixin.Copyable @(x) any(strcmpi(x,{'ImaginaryTimeEvolution','RealTimeEvolution'}))); addParameter(p, 'TimeStep', 5E-4,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'SimulationTime', 3,... + addParameter(p, 'SimulationTime', 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,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); addParameter(p, 'DebugMode', false,... @islogical); addParameter(p, 'SaveData', false,... @@ -68,6 +74,8 @@ classdef DipolarGas < handle & matlab.mixin.Copyable this.SimulationMode = p.Results.SimulationMode; this.TimeStep = p.Results.TimeStep; this.SimulationTime = p.Results.SimulationTime; + this.EnergyTolerance = p.Results.EnergyTolerance; + this.MinimumTimeStep = p.Results.MinimumTimeStep; this.DebugMode = p.Results.DebugMode; this.DoSave = p.Results.SaveData; diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m index c23e0c4..bf6114e 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m @@ -1,34 +1,34 @@ function [Params] = setupParameters(this) -CONSTANTS = Helper.PhysicsConstants; -hbar = CONSTANTS.PlanckConstantReduced; % [J.s] -kbol = CONSTANTS.BoltzmannConstant; % [J/K] -mu0 = CONSTANTS.VacuumPermeability; % [N/A^2] -muB = CONSTANTS.BohrMagneton; % [J/T] -a0 = CONSTANTS.BohrRadius; % [m] -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. +CONSTANTS = Helper.PhysicsConstants; +hbar = CONSTANTS.PlanckConstantReduced; % [J.s] +kbol = CONSTANTS.BoltzmannConstant; % [J/K] +mu0 = CONSTANTS.VacuumPermeability; % [N/A^2] +muB = CONSTANTS.BohrMagneton; % [J/T] +a0 = CONSTANTS.BohrRadius; % [m] +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) % Number of points in each direction -Params.Nx = this.NumberOfGridPoints(1); -Params.Ny = this.NumberOfGridPoints(2); -Params.Nz = this.NumberOfGridPoints(3); +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); +Params.Lx = this.Dimensions(1); +Params.Ly = this.Dimensions(2); +Params.Lz = this.Dimensions(3); % Masses -Params.m = CONSTANTS.Dy164Mass; -l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length +Params.m = CONSTANTS.Dy164Mass; +l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length % Atom numbers Params.N = this.NumberOfAtoms; % 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; % Dipole lengths (units of muB) @@ -38,35 +38,34 @@ Params.mu = CONSTANTS.DyMagneticMoment; Params.as = this.ScatteringLength*a0; % Trapping frequencies -Params.wx = 2*pi*this.TrapFrequencies(1); -Params.wy = 2*pi*this.TrapFrequencies(2); -Params.wz = 2*pi*this.TrapFrequencies(3); +Params.wx = 2*pi*this.TrapFrequencies(1); +Params.wy = 2*pi*this.TrapFrequencies(2); +Params.wz = 2*pi*this.TrapFrequencies(3); % 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.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.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 + % ================ Parameters defined by those above ================ % % Contact interaction strength (units of l0/m) -Params.gs = 4*pi*Params.as/l0; +Params.gs = 4*pi*Params.as/l0; % 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 -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 -Params.gx=(Params.wx/w0)^2; -Params.gy=(Params.wy/w0)^2; -Params.gz=(Params.wz/w0)^2; +Params.gx = (Params.wx/w0)^2; +Params.gy = (Params.wy/w0)^2; +Params.gz = (Params.wz/w0)^2; % == Calculate LHY correction == % eps_dd = Params.add/Params.as; @@ -83,6 +82,11 @@ end Params.gammaQF = 128/3*sqrt(pi*(Params.as/l0)^5)*Q5; % Loading the rest into Params -Params.hbar = hbar; Params.kbol = kbol; Params.mu0 = mu0; Params.muB = muB; Params.a0 = a0; -Params.w0 = w0; Params.l0 = l0; +Params.hbar = hbar; +Params.kbol = kbol; +Params.mu0 = mu0; +Params.muB = muB; +Params.a0 = a0; +Params.w0 = w0; +Params.l0 = l0; end \ No newline at end of file