diff --git a/Dipolar Gas Simulator/+Helper/PhysicsConstants.m b/Dipolar Gas Simulator/+Helper/PhysicsConstants.m index 9a37257..14fbc82 100644 --- a/Dipolar Gas Simulator/+Helper/PhysicsConstants.m +++ b/Dipolar Gas Simulator/+Helper/PhysicsConstants.m @@ -25,6 +25,7 @@ classdef PhysicsConstants < handle % Dy specific constants Dy164Mass = 163.929174751*1.660539066E-27; Dy164IsotopicAbundance = 0.2826; + DyMagneticMoment = 9.93*9.274009994E-24; end methods diff --git a/Dipolar Gas Simulator/+Plotter/visualizeDDIPotential.m b/Dipolar Gas Simulator/+Plotter/visualizeDDIPotential.m new file mode 100644 index 0000000..eae0eab --- /dev/null +++ b/Dipolar Gas Simulator/+Plotter/visualizeDDIPotential.m @@ -0,0 +1,94 @@ +function visualizeDDIPotential(Transf) + + x = Transf.x; + y = Transf.y; + z = Transf.z; + + [X,Y] = meshgrid(x,y); + + height = 20; + width = 45; + figure(1) + clf + set(gcf, 'Units', 'centimeters') + set(gcf, 'Position', [2 4 width height]) + set(gcf, 'PaperPositionMode', 'auto') + + subplot(1,2,1) + zlin = ones(size(X, 1)) * z(1); % Generate z data + mesh(x, y, zlin, 'EdgeColor','b') % Plot the surface + hold on + + for idx = 2:length(z) + zlin = ones(size(X, 1))* z(idx); % Generate z data + mesh(x, y, zlin, 'EdgeColor','b') % Plot the surface + end + + % set the axes labels' properties + xlabel(gca, {'$x / l_o ~ (m)$'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + ylabel(gca, {'$y / l_o ~ (m)$'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + zlabel(gca, {'$z / l_o ~ (m)$'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + title(gca, 'Real Space', ... + 'Interpreter', 'tex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + + x = Transf.kx; + y = Transf.ky; + z = Transf.kz; + + [X,Y] = meshgrid(x,y); + + subplot(1,2,2) + zlin = ones(size(X, 1)) * z(1); % Generate z data + mesh(x, y, zlin, 'EdgeColor','b') % Plot the surface + hold on + + for idx = 2:length(z) + zlin = ones(size(X, 1))* z(idx); % Generate z data + mesh(x, y, zlin, 'EdgeColor','b') % Plot the surface + end + + % set the axes labels' properties + xlabel(gca, {'$k_x / l_o ~ (m^{-1})$'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + ylabel(gca, {'$k_y / l_o ~ (m^{-1})$'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + zlabel(gca, {'$k_z / l_o ~ (m^{-1})$'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + title(gca, 'Fourier Space', ... + 'Interpreter', 'tex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') +end \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Plotter/visualizeSpace.m b/Dipolar Gas Simulator/+Plotter/visualizeSpace.m index ac5646c..5785318 100644 --- a/Dipolar Gas Simulator/+Plotter/visualizeSpace.m +++ b/Dipolar Gas Simulator/+Plotter/visualizeSpace.m @@ -50,12 +50,6 @@ function visualizeSpace(Transf) 'FontWeight', 'normal', ... 'FontAngle', 'normal') - x = Transf.kx; - y = Transf.ky; - z = Transf.kz; - - [X,Y] = meshgrid(x,y); - subplot(1,2,2) zlin = ones(size(X, 1)) * z(1); % Generate z data mesh(x, y, zlin, 'EdgeColor','b') % Plot the surface diff --git a/Dipolar Gas Simulator/+Plotter/visualizeTrapPotential.m b/Dipolar Gas Simulator/+Plotter/visualizeTrapPotential.m new file mode 100644 index 0000000..5cc278a --- /dev/null +++ b/Dipolar Gas Simulator/+Plotter/visualizeTrapPotential.m @@ -0,0 +1,95 @@ +function visualizeTrapPotential(V,Params,Transf) + + set(0,'defaulttextInterpreter','latex') + set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex'); + + format long + x = Transf.x*Params.l0*1e6; + y = Transf.y*Params.l0*1e6; + z = Transf.z*Params.l0*1e6; + + dx = x(2)-x(1); dy = y(2)-y(1); dz = z(2)-z(1); + + %Plotting + height = 10; + width = 35; + figure(1) + clf + set(gcf, 'Units', 'centimeters') + set(gcf, 'Position', [2 4 width height]) + set(gcf, 'PaperPositionMode', 'auto') + + subplot(1,3,1) + n = V; + nxz = squeeze(trapz(n*dy,2)); + nyz = squeeze(trapz(n*dx,1)); + nxy = squeeze(trapz(n*dz,3)); + + plotxz = pcolor(x,z,nxz'); + set(plotxz, 'EdgeColor', 'none'); + xlabel(gca, {'$x$ ($\mu$m)'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + ylabel(gca, {'$z$ ($\mu$m)'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + title(gca, {'$V_{xz}$'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + colorbar + + subplot(1,3,2) + plotyz = pcolor(y,z,nyz'); + set(plotyz, 'EdgeColor', 'none'); + xlabel(gca, {'$y$ ($\mu$m)'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + ylabel(gca, {'$z$ ($\mu$m)'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + title(gca, {'$V_{yz}$'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + colorbar + + subplot(1,3,3) + plotxy = pcolor(x,y,nxy'); + set(plotxy, 'EdgeColor', 'none'); + xlabel(gca, {'$x$ ($\mu$m)'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + ylabel(gca, {'$y$ ($\mu$m)'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + title(gca, {'$V_{xy}$'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + colorbar +end \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Plotter/visualizeWavefunction.m b/Dipolar Gas Simulator/+Plotter/visualizeWavefunction.m new file mode 100644 index 0000000..807e51e --- /dev/null +++ b/Dipolar Gas Simulator/+Plotter/visualizeWavefunction.m @@ -0,0 +1,92 @@ +function visualizeWavefunction(psi,Params,Transf) + + set(0,'defaulttextInterpreter','latex') + set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex'); + + format long + x = Transf.x*Params.l0*1e6; + y = Transf.y*Params.l0*1e6; + z = Transf.z*Params.l0*1e6; + + dx = x(2)-x(1); dy = y(2)-y(1); dz = z(2)-z(1); + + %Plotting + height = 10; + width = 35; + figure(1) + clf + set(gcf, 'Units', 'centimeters') + set(gcf, 'Position', [2 4 width height]) + set(gcf, 'PaperPositionMode', 'auto') + + subplot(1,3,1) + n = abs(psi).^2; + nxz = squeeze(trapz(n*dy,2)); + nyz = squeeze(trapz(n*dx,1)); + nxy = squeeze(trapz(n*dz,3)); + + plotxz = pcolor(x,z,nxz'); + set(plotxz, 'EdgeColor', 'none'); + xlabel(gca, {'$x$ ($\mu$m)'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + ylabel(gca, {'$z$ ($\mu$m)'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + title(gca, {'$|\Psi_{xz}|^2$'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + + subplot(1,3,2) + plotyz = pcolor(y,z,nyz'); + set(plotyz, 'EdgeColor', 'none'); + xlabel(gca, {'$y$ ($\mu$m)'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + ylabel(gca, {'$z$ ($\mu$m)'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + title(gca, {'$|\Psi_{yz}|^2$'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + + subplot(1,3,3) + plotxy = pcolor(x,y,nxy'); + set(plotxy, 'EdgeColor', 'none'); + xlabel(gca, {'$x$ ($\mu$m)'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + ylabel(gca, {'$y$ ($\mu$m)'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') + title(gca, {'$|\Psi_{xy}|^2$'}, ... + 'Interpreter', 'latex', ... + 'FontName', 'Times New Roman', ... + 'FontSize', 14, ... + 'FontWeight', 'normal', ... + 'FontAngle', 'normal') +end \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Scripts/test.m b/Dipolar Gas Simulator/+Scripts/test.m index 45fce31..23d5ea2 100644 --- a/Dipolar Gas Simulator/+Scripts/test.m +++ b/Dipolar Gas Simulator/+Scripts/test.m @@ -5,13 +5,23 @@ %% - Create Simulator, Potential and Calculator object with specified options OptionsStruct = struct; -OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' + +OptionsStruct.NumberOfAtoms = 1E6; +OptionsStruct.DipolarAngle = pi/2; +OptionsStruct.ScatteringLength = 86; +OptionsStruct.TrapFrequencies = [125, 125, 250]; +OptionsStruct.NumberOfPoints = [64, 64, 48]; +OptionsStruct.Dimensions = [40, 40, 20]; OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.PotentialType = 'Harmonic'; + +OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' OptionsStruct.TimeStep = 50e-06; % in s OptionsStruct.SimulationTime = 4e-03; % in s + OptionsStruct.SaveData = true; OptionsStruct.SaveDirectory = './Data'; + options = Helper.convertstruct2cell(OptionsStruct); clear OptionsStruct @@ -23,4 +33,8 @@ pot = Simulator.Potentials(options{:}); [Params, Transf, psi, V, VDk] = sim.runSimulation(calc); %% - Plot results -Plotter.visualizeSpace(Transf) \ No newline at end of file +% Plotter.visualizeSpace(Transf) + +Plotter.visualizeTrapPotential(V,Params,Transf) + +% Plotter.visualizeWavefunction(psi,Params,Transf) \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m index 141de91..3128b91 100644 --- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m +++ b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m @@ -1,7 +1,7 @@ -function VDkSemi = calculateNumericalHankelTransform(kr, kz, Rmax, Zmax, Nr) +function VDkSemi = calculateNumericalHankelTransform(this,kr,kz,Rmax,Zmax,Nr) % accuracy inputs for numerical integration -if(nargin==4) +if(nargin==5) Nr = 5e4; end diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateVDCutoff.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateVDCutoff.m index aac3e4b..ef41bf3 100644 --- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateVDCutoff.m +++ b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateVDCutoff.m @@ -1,4 +1,4 @@ -function VDk = calculateVDCutoff(this, Params, Transf) +function VDk = calculateVDCutoff(this,Params,Transf,TransfRad) % makes the dipolar interaction matrix, size numel(Params.kr) * numel(Params.kz) % Rmax and Zmax are the interaction cutoffs % VDk needs to be multiplied by Cdd @@ -16,10 +16,7 @@ VDk = VDk + exp(-Zcutoff*sqrt(Transf.KX.^2+Transf.KY.^2)).*( sinsq .* cos(Zcut % Nonanalytic part % For a cylindrical cutoff, we need to construct a kr grid based on the 3D parameters using Bessel quadrature -Params.Lr = 0.5*min(Params.Lx,Params.Ly); -Params.Nr = max(Params.Nx,Params.Ny); -[TransfRad] = SetupSpaceRadial(Params); -VDkNon = this.calculateNumericalHankelTransform(TransfRad.kr, TransfRad.kz, TransfRad.Rmax, Zcutoff); +VDkNon = this.calculateNumericalHankelTransform(TransfRad.kr, TransfRad.kz, TransfRad.Rmax, Zcutoff); % Interpolating the nonanalytic part onto 3D grid fullkr = [-flip(TransfRad.kr)',TransfRad.kr']; diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/DipolarGas.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/DipolarGas.m index d9cef99..db39c89 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/DipolarGas.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/DipolarGas.m @@ -1,20 +1,21 @@ classdef DipolarGas < handle & matlab.mixin.Copyable properties (Access = public) + NumberOfAtoms; + DipolarAngle; + ScatteringLength; + TrapFrequencies; + NumberOfPoints; + Dimensions; + CutoffType; + PotentialType; + SimulationMode; TimeStep; SimulationTime; - NumberOfAtoms; - % Etol = 5e-10; % Tolerances - % rtol = 1e-5; - % cut_off = 2e6; % sometimes the imaginary time gets a little stuck - % even though the solution is good, this just stops it going on forever - % mindt = 1e-6; %Minimum size for a time step using adaptive dt - %Flags - DebugMode; DoSave; SaveDirectory; @@ -28,27 +29,41 @@ classdef DipolarGas < handle & matlab.mixin.Copyable p = inputParser; p.KeepUnmatched = true; + addParameter(p, 'NumberOfAtoms', 1E6,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'DipolarAngle', pi/2,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > -2*pi) && (x < 2*pi))); + addParameter(p, 'ScatteringLength', 120,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > -150) && (x < 150))); + addParameter(p, 'TrapFrequencies', 100 * ones(1,3),... + @(x) assert(isnumeric(x) && isvector(x) && all(x > 0))); + addParameter(p, 'NumberOfPoints', 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))); addParameter(p, 'SimulationMode', 'ImaginaryTimeEvolution',... @(x) any(strcmpi(x,{'ImaginaryTimeEvolution','RealTimeEvolution'}))); - addParameter(p, 'NumberOfAtoms', 5000,... + addParameter(p, 'TimeStep', 5E-4,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'TimeStep', 0.0005,... + addParameter(p, 'SimulationTime', 3,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'SimulationTime', 3e-03,... - @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'DebugMode', false,... + addParameter(p, 'DebugMode', false,... @islogical); - addParameter(p, 'SaveData', false,... + addParameter(p, 'SaveData', false,... @islogical); - addParameter(p, 'SaveDirectory', './Data',... + addParameter(p, 'SaveDirectory', './Data',... @ischar); p.parse(varargin{:}); - this.SimulationMode = p.Results.SimulationMode; - this.ErrorEstimationMethod= p.Results.ErrorEstimationMethod; - this.NumberOfAtoms = p.Results.NumberOfAtoms; + this.DipolarAngle = p.Results.DipolarAngle; + this.ScatteringLength = p.Results.ScatteringLength; + this.TrapFrequencies = p.Results.TrapFrequencies; + this.NumberOfPoints = p.Results.NumberOfPoints; + this.Dimensions = p.Results.Dimensions; + + this.SimulationMode = p.Results.SimulationMode; this.TimeStep = p.Results.TimeStep; this.SimulationTime = p.Results.SimulationTime; @@ -81,7 +96,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable ret = this.SimulationTime; end function set.NumberOfAtoms(this, val) - assert(val <= 50000, '!!Not time efficient to compute for atom numbers larger than 50,000!!'); + assert(val <= 1E6, '!!Not time efficient to compute for atom numbers larger than 1,000,000!!'); this.NumberOfAtoms = val; end function ret = get.NumberOfAtoms(this) diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m index 7ee30c5..f7a5e32 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m @@ -1,4 +1,4 @@ -function [psi,V,VDk] = initialize(this,calcObj,Params,Transf) +function [psi,V,VDk] = initialize(this,calcObj,Params,Transf,TransfRad) format long X = Transf.X; Y = Transf.Y; Z = Transf.Z; @@ -11,12 +11,12 @@ if isfile(strcat(this.SaveDirectory, '/VDk_M.mat')) VDk = load(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat'))); VDk = VDk.VDk; else - VDk = calcObj.calculateVDCutoff(Params, Transf); + VDk = calcObj.calculateVDCutoff(Params,Transf,TransfRad); save(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat')),'VDk'); end disp('Finished DDI') % == Setting up the initial wavefunction == % -psi = this.SetupWavefunction(); +psi = this.setupWavefunction(Params,Transf); end \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/runSimulation.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/runSimulation.m index 9e23767..74ef82b 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/runSimulation.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/runSimulation.m @@ -4,11 +4,12 @@ function [Params, Transf, psi,V,VDk] = runSimulation(this,calcObj) % --- Set up spatial grids and transforms --- [Transf] = this.setupSpace(Params); - + [TransfRad] = this.setupSpaceRadial(Params); + % --- Initialize --- mkdir(sprintf(this.SaveDirectory)) - [psi,V,VDk] = this.initialize(calcObj,Params,Transf); + [psi,V,VDk] = this.initialize(calcObj,Params,Transf,TransfRad); Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = []; t_idx = 1; %Start at t = 0; diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m index fdd000a..7bb6005 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m @@ -1,6 +1,5 @@ function [Params] = setupParameters(this) -%========= Constants =========% CONSTANTS = Helper.PhysicsConstants; hbar = CONSTANTS.PlanckConstantReduced; % [J.s] kbol = CONSTANTS.BoltzmannConstant; % [J/K] @@ -11,44 +10,47 @@ 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 = 64; -Params.Ny = 64; -Params.Nz = 48; +Params.Nx = this.NumberOfPoints(1); +Params.Ny = this.NumberOfPoints(2); +Params.Nz = this.NumberOfPoints(3); % Dimensions (in units of l0) -Params.Lx = 40; -Params.Ly = 40; -Params.Lz = 20; +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 % Atom numbers -% Params.ppum = 2500; % particles per micron -% Params.N = Params.Lz*Params.ppum*l0*1e6; -Params.N = 10^6; +Params.N = this.NumberOfAtoms; % Dipole angle -Params.theta = pi/2; % pi/2 dipoles along x, theta=0 dipoles along z +Params.theta = this.DipolarAngle; % pi/2 dipoles along x, theta=0 dipoles along z +Params.phi = 0; % Dipole lengths (units of muB) -Params.mu = 9.93*muB; +Params.mu = CONSTANTS.DyMagneticMoment; % Scattering lengths -Params.as = 86*a0; +Params.as = this.ScatteringLength*a0; % Trapping frequencies -Params.wx = 2*pi*125; -Params.wy = 2*pi*125; -Params.wz = 2*pi*250; +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.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 % ================ Parameters defined by those above ================ % @@ -59,7 +61,7 @@ Params.gs = 4*pi*Params.as/l0; 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; diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m index 2000863..f9349e2 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m @@ -1,10 +1,14 @@ function [Transf] = setupSpaceRadial(this,Params,morder) + +Params.Lr = 0.5*min(Params.Lx,Params.Ly); +Params.Nr = max(Params.Nx,Params.Ny); + Zmax = 0.5*Params.Lz; Rmax = Params.Lr; Nz = Params.Nz; Nr = Params.Nr; -if(nargin==1) +if(nargin==2) morder=0; %only do Bessel J0 end @@ -42,7 +46,7 @@ Transf.dz=z(2)-z(1); Transf.dkz=kz(2)-kz(1); %b1=Transf; -function s_HT = hankelmatrix(order, rmax, Nr, eps_roots) +function s_HT = hankelmatrix(order,rmax,Nr,eps_roots) %HANKEL_MATRIX: Generates data to use for Hankel Transforms % % s_HT = hankel_matrix(order, rmax, Nr, eps_roots) @@ -132,7 +136,7 @@ s_HT.wk=2./((s_HT.rmax)^2*abs(Jp1).^2); return -function z = bessel_zeros(d, a, n, e) +function z = bessel_zeros(d,a,n,e) %BESSEL_ZEROS: Finds the first n zeros of a bessel function % % z = bessel_zeros(d, a, n, e) @@ -293,7 +297,7 @@ function FI = fi(y) end return -function Jr = bessr(d, a, x) +function Jr = bessr(d,a,x) switch (d) case (1) Jr = besselj(a, x)./besselj(a+1, x); diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupWavefunction.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupWavefunction.m index 756a553..e16cf48 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupWavefunction.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupWavefunction.m @@ -1,4 +1,6 @@ -function [psi] = setupWavefunction(this) +function [psi] = setupWavefunction(this,Params,Transf) + +X = Transf.X; Y = Transf.Y; Z = Transf.Z; ellx = sqrt(Params.hbar/(Params.m*Params.wx))/Params.l0; elly = sqrt(Params.hbar/(Params.m*Params.wy))/Params.l0; @@ -7,15 +9,15 @@ ellz = sqrt(Params.hbar/(Params.m*Params.wz))/Params.l0; Rx = 8*ellx; Ry = 8*elly; Rz = 8*ellz; X0 = 0.0*Transf.Xmax; Y0 = 0.0*Transf.Ymax; Z0 = 0*Transf.Zmax; -psi = exp(-(X-X0).^2/Rx^2-(Y-Y0).^2/Ry^2-(Z-Z0).^2/Rz^2); +psi = exp(-(X-X0).^2/Rx^2-(Y-Y0).^2/Ry^2-(Z-Z0).^2/Rz^2); cur_norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; -psi = psi/sqrt(cur_norm); +psi = psi/sqrt(cur_norm); -%add some noise -r = normrnd(0,1,size(X)); +% add some noise +r = normrnd(0,1,size(X)); theta = rand(size(X)); noise = r.*exp(2*pi*1i*theta); -psi = 0*psi + 1*noise; +psi = psi + 0.01*noise; Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; psi = sqrt(Params.N)*psi/sqrt(Norm);