From cfe15aa13d414cf654ebfc69d55c21afb64e6c3c Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Thu, 13 Jun 2024 11:30:12 +0200 Subject: [PATCH] Latest working version. --- .../+Plotter/visualizeSpace.m | 6 ++ .../+Plotter/visualizeTrapPotential.m | 17 ++-- .../+Plotter/visualizeWavefunction.m | 4 +- Dipolar Gas Simulator/+Scripts/test.m | 17 ++-- .../+Simulator/@Calculator/Calculator.m | 20 ++++- .../@Calculator/calculateVDCutoff.m | 80 +++++++++++++------ .../+Simulator/@DipolarGas/DipolarGas.m | 18 +++-- .../+Simulator/@DipolarGas/initialize.m | 2 +- .../+Simulator/@DipolarGas/setupParameters.m | 10 +-- .../+Simulator/@Potentials/Potentials.m | 5 +- 10 files changed, 121 insertions(+), 58 deletions(-) diff --git a/Dipolar Gas Simulator/+Plotter/visualizeSpace.m b/Dipolar Gas Simulator/+Plotter/visualizeSpace.m index 5785318..ac5646c 100644 --- a/Dipolar Gas Simulator/+Plotter/visualizeSpace.m +++ b/Dipolar Gas Simulator/+Plotter/visualizeSpace.m @@ -50,6 +50,12 @@ 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 index 5cc278a..b561e06 100644 --- a/Dipolar Gas Simulator/+Plotter/visualizeTrapPotential.m +++ b/Dipolar Gas Simulator/+Plotter/visualizeTrapPotential.m @@ -1,5 +1,4 @@ function visualizeTrapPotential(V,Params,Transf) - set(0,'defaulttextInterpreter','latex') set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex'); @@ -9,22 +8,26 @@ function visualizeTrapPotential(V,Params,Transf) 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; + width = 45; figure(1) clf set(gcf, 'Units', 'centimeters') - set(gcf, 'Position', [2 4 width height]) + set(gcf, 'Position', [2 8 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)); + nxz = nxz./max(nxz(:)); + nyz = nyz./max(nyz(:)); + nxy = nxy./max(nxy(:)); + plotxz = pcolor(x,z,nxz'); set(plotxz, 'EdgeColor', 'none'); xlabel(gca, {'$x$ ($\mu$m)'}, ... @@ -46,7 +49,7 @@ function visualizeTrapPotential(V,Params,Transf) 'FontWeight', 'normal', ... 'FontAngle', 'normal') colorbar - + subplot(1,3,2) plotyz = pcolor(y,z,nyz'); set(plotyz, 'EdgeColor', 'none'); @@ -69,7 +72,7 @@ function visualizeTrapPotential(V,Params,Transf) 'FontWeight', 'normal', ... 'FontAngle', 'normal') colorbar - + subplot(1,3,3) plotxy = pcolor(x,y,nxy'); set(plotxy, 'EdgeColor', 'none'); diff --git a/Dipolar Gas Simulator/+Plotter/visualizeWavefunction.m b/Dipolar Gas Simulator/+Plotter/visualizeWavefunction.m index 807e51e..db29f0f 100644 --- a/Dipolar Gas Simulator/+Plotter/visualizeWavefunction.m +++ b/Dipolar Gas Simulator/+Plotter/visualizeWavefunction.m @@ -12,11 +12,11 @@ function visualizeWavefunction(psi,Params,Transf) %Plotting height = 10; - width = 35; + width = 45; figure(1) clf set(gcf, 'Units', 'centimeters') - set(gcf, 'Position', [2 4 width height]) + set(gcf, 'Position', [2 8 width height]) set(gcf, 'PaperPositionMode', 'auto') subplot(1,3,1) diff --git a/Dipolar Gas Simulator/+Scripts/test.m b/Dipolar Gas Simulator/+Scripts/test.m index 23d5ea2..1b6d93f 100644 --- a/Dipolar Gas Simulator/+Scripts/test.m +++ b/Dipolar Gas Simulator/+Scripts/test.m @@ -7,13 +7,14 @@ OptionsStruct = struct; OptionsStruct.NumberOfAtoms = 1E6; -OptionsStruct.DipolarAngle = pi/2; +OptionsStruct.DipolarPolarAngle = pi/2; +OptionsStruct.DipolarAzimuthAngle = 0; OptionsStruct.ScatteringLength = 86; OptionsStruct.TrapFrequencies = [125, 125, 250]; -OptionsStruct.NumberOfPoints = [64, 64, 48]; +OptionsStruct.NumberOfGridPoints = [64, 64, 48]; OptionsStruct.Dimensions = [40, 40, 20]; OptionsStruct.CutoffType = 'Cylindrical'; -OptionsStruct.PotentialType = 'Harmonic'; +OptionsStruct.TrapPotentialType = 'Harmonic'; OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' OptionsStruct.TimeStep = 50e-06; % in s @@ -32,9 +33,9 @@ pot = Simulator.Potentials(options{:}); %-% Run Simulation %-% [Params, Transf, psi, V, VDk] = sim.runSimulation(calc); -%% - Plot results -% Plotter.visualizeSpace(Transf) - +%% - Plot numerical grid +Plotter.visualizeSpace(Transf) +%% - Plot trap potential Plotter.visualizeTrapPotential(V,Params,Transf) - -% Plotter.visualizeWavefunction(psi,Params,Transf) \ No newline at end of file +%% - Plot initial wavefunction +Plotter.visualizeWavefunction(psi,Params,Transf) \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/Calculator.m b/Dipolar Gas Simulator/+Simulator/@Calculator/Calculator.m index 2195bb8..16522d3 100644 --- a/Dipolar Gas Simulator/+Simulator/@Calculator/Calculator.m +++ b/Dipolar Gas Simulator/+Simulator/@Calculator/Calculator.m @@ -2,15 +2,23 @@ classdef Calculator < handle & matlab.mixin.Copyable properties (Access = private) - CalculatorDefaults = struct('OrderParameter', 1, ... + CalculatorDefaults = struct('ChemicalPotential', 1, ... + 'EnergyComponents', 1, ... + 'NormalizedResiduals', 1, ... + 'OrderParameter', 1, ... 'PhaseCoherence', 1, ... + 'TotalEnergy', 1, ... 'CutoffType', 'Cylindrical'); end properties (Access = public) + ChemicalPotential; + EnergyComponents; + NormalizedResiduals; OrderParameter; PhaseCoherence; + TotalEnergy; CutoffType; end @@ -19,14 +27,22 @@ classdef Calculator < handle & matlab.mixin.Copyable methods function this = Calculator(varargin) + this.ChemicalPotential = this.CalculatorDefaults.ChemicalPotential; + this.EnergyComponents = this.CalculatorDefaults.EnergyComponents; + this.NormalizedResiduals = this.CalculatorDefaults.NormalizedResiduals; this.OrderParameter = this.CalculatorDefaults.OrderParameter; this.PhaseCoherence = this.CalculatorDefaults.PhaseCoherence; + this.TotalEnergy = this.CalculatorDefaults.TotalEnergy; this.CutoffType = this.CalculatorDefaults.CutoffType; end function restoreDefaults(this) - this.OrderParameter = this.PotentialDefaults.OrderParameter; + this.ChemicalPotential = this.CalculatorDefaults.ChemicalPotential; + this.EnergyComponents = this.CalculatorDefaults.EnergyComponents; + this.NormalizedResiduals = this.CalculatorDefaults.NormalizedResiduals; + this.OrderParameter = this.CalculatorDefaults.OrderParameter; this.PhaseCoherence = this.CalculatorDefaults.PhaseCoherence; + this.TotalEnergy = this.CalculatorDefaults.TotalEnergy; this.CutoffType = this.CalculatorDefaults.CutoffType; end diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateVDCutoff.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateVDCutoff.m index ef41bf3..3c49b4e 100644 --- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateVDCutoff.m +++ b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateVDCutoff.m @@ -4,27 +4,61 @@ function VDk = calculateVDCutoff(this,Params,Transf,TransfRad) % VDk needs to be multiplied by Cdd % approach is that of Lu, PRA 82, 023622 (2010) -Zcutoff = Params.Lz/2; -alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2)); -alph(1) = pi/2; +% == Calulating the DDI potential in Fourier space with appropriate cutoff == % +% Cylindrical (semianalytic) +% Cylindrical infinite Z, polarized along x (analytic) +% Spherical -% Analytic part of cutoff for slice 0 0))); - addParameter(p, 'DipolarAngle', pi/2,... + addParameter(p, 'DipolarPolarAngle', pi/2,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > -2*pi) && (x < 2*pi))); + addParameter(p, 'DipolarAzimuthAngle', 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),... + 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))); @@ -57,10 +58,11 @@ classdef DipolarGas < handle & matlab.mixin.Copyable p.parse(varargin{:}); this.NumberOfAtoms = p.Results.NumberOfAtoms; - this.DipolarAngle = p.Results.DipolarAngle; + this.DipolarPolarAngle = p.Results.DipolarPolarAngle; + this.DipolarAzimuthAngle = p.Results.DipolarAzimuthAngle; this.ScatteringLength = p.Results.ScatteringLength; this.TrapFrequencies = p.Results.TrapFrequencies; - this.NumberOfPoints = p.Results.NumberOfPoints; + this.NumberOfGridPoints = p.Results.NumberOfGridPoints; this.Dimensions = p.Results.Dimensions; this.SimulationMode = p.Results.SimulationMode; diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m index f7a5e32..d7d5eb0 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m @@ -14,7 +14,7 @@ else VDk = calcObj.calculateVDCutoff(Params,Transf,TransfRad); save(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat')),'VDk'); end -disp('Finished DDI') +fprintf('Computed and saved DDI potential in Fourier space with %s cutoff.', calcObj.CutoffType) % == Setting up the initial wavefunction == % psi = this.setupWavefunction(Params,Transf); diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m index 7bb6005..c23e0c4 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m @@ -11,9 +11,9 @@ 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.NumberOfPoints(1); -Params.Ny = this.NumberOfPoints(2); -Params.Nz = this.NumberOfPoints(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); @@ -28,8 +28,8 @@ l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length Params.N = this.NumberOfAtoms; % Dipole angle -Params.theta = this.DipolarAngle; % pi/2 dipoles along x, theta=0 dipoles along z -Params.phi = 0; +Params.theta = this.DipolarPolarAngle; % pi/2 dipoles along x, theta=0 dipoles along z +Params.phi = this.DipolarAzimuthAngle; % Dipole lengths (units of muB) Params.mu = CONSTANTS.DyMagneticMoment; diff --git a/Dipolar Gas Simulator/+Simulator/@Potentials/Potentials.m b/Dipolar Gas Simulator/+Simulator/@Potentials/Potentials.m index 8421484..4718c85 100644 --- a/Dipolar Gas Simulator/+Simulator/@Potentials/Potentials.m +++ b/Dipolar Gas Simulator/+Simulator/@Potentials/Potentials.m @@ -2,11 +2,12 @@ classdef Potentials < handle & matlab.mixin.Copyable properties (Access = private) - PotentialDefaults = struct('TrapFrequency', 100e3); - + PotentialDefaults = struct('TrapPotentialType', 'Harmonic', ... + 'TrapFrequency', 100e3); end properties (Access = public) + TrapPotentialType; TrapFrequency; end