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

Latest working version.

Browse Source
This commit is contained in:
Karthik 2024-06-13 11:30:12 +02:00
parent 7cab430c05
commit cfe15aa13d
10 changed files with 121 additions and 58 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
Dipolar Gas Simulator/+Plotter/visualizeSpace.m
View File

@ -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

9
Dipolar Gas Simulator/+Plotter/visualizeTrapPotential.m
View File

@ -1,5 +1,4 @@
function visualizeTrapPotential(V,Params,Transf)
set(0,'defaulttextInterpreter','latex')
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
@ -12,11 +11,11 @@ function visualizeTrapPotential(V,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)
@ -25,6 +24,10 @@ function visualizeTrapPotential(V,Params,Transf)
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)'}, ...

4
Dipolar Gas Simulator/+Plotter/visualizeWavefunction.m
View File

@ -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)

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

@ -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)
%% - Plot initial wavefunction
Plotter.visualizeWavefunction(psi,Params,Transf)

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

@ -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

74
Dipolar Gas Simulator/+Simulator/@Calculator/calculateVDCutoff.m
View File

@ -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<z<Zmax, 0<r<Inf Ronen, PRL 98, 030406 (2007)
cossq = cos(alph).^2;
VDk = cossq-1/3;
sinsq = 1 - cossq;
VDk = VDk + exp(-Zcutoff*sqrt(Transf.KX.^2+Transf.KY.^2)).*( sinsq .* cos(Zcutoff * Transf.KZ) - sqrt(sinsq.*cossq).*sin(Zcutoff * Transf.KZ) );
switch this.CutoffType
case 'Cylindrical' %Cylindrical (semianalytic)
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;
% Nonanalytic part
% For a cylindrical cutoff, we need to construct a kr grid based on the 3D parameters using Bessel quadrature
VDkNon = this.calculateNumericalHankelTransform(TransfRad.kr, TransfRad.kz, TransfRad.Rmax, Zcutoff);
% Analytic part of cutoff for slice 0<z<Zmax, 0<r<Inf Ronen, PRL 98, 030406 (2007)
cossq = cos(alph).^2;
VDk = cossq-1/3;
sinsq = 1 - cossq;
VDk = VDk + exp(-Zcutoff*sqrt(Transf.KX.^2+Transf.KY.^2)).*( sinsq .* cos(Zcutoff * Transf.KZ) - sqrt(sinsq.*cossq).*sin(Zcutoff * Transf.KZ) );
% Interpolating the nonanalytic part onto 3D grid
fullkr = [-flip(TransfRad.kr)',TransfRad.kr'];
[KR,KZ] = ndgrid(fullkr,TransfRad.kz);
[KX3D,KY3D,KZ3D] = ndgrid(ifftshift(Transf.kx),ifftshift(Transf.ky),ifftshift(Transf.kz));
KR3D = sqrt(KX3D.^2 + KY3D.^2);
fullVDK = [flip(VDkNon',2),VDkNon']';
VDkNon = interpn(KR,KZ,fullVDK,KR3D,KZ3D,'spline',0); %Last argument is -1/3 for full VDk. 0 for nonanalytic piece?
VDkNon = fftshift(VDkNon);
% Nonanalytic part
% For a cylindrical cutoff, we need to construct a kr grid based on the 3D parameters using Bessel quadrature
VDkNon = this.calculateNumericalHankelTransform(TransfRad.kr, TransfRad.kz, TransfRad.Rmax, Zcutoff);
VDk = VDk + VDkNon;
% Interpolating the nonanalytic part onto 3D grid
fullkr = [-flip(TransfRad.kr)',TransfRad.kr'];
[KR,KZ] = ndgrid(fullkr,TransfRad.kz);
[KX3D,KY3D,KZ3D] = ndgrid(ifftshift(Transf.kx),ifftshift(Transf.ky),ifftshift(Transf.kz));
KR3D = sqrt(KX3D.^2 + KY3D.^2);
fullVDK = [flip(VDkNon',2),VDkNon']';
VDkNon = interpn(KR,KZ,fullVDK,KR3D,KZ3D,'spline',0); %Last argument is -1/3 for full VDk. 0 for nonanalytic piece?
VDkNon = fftshift(VDkNon);
VDk = VDk + VDkNon;
case 'CylindricalInfiniteZ' %Cylindrical infinite Z, polarized along x -- PRA 107, 033301 (2023)
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;
rhoc = max([abs(Transf.x),abs(Transf.y)]);
KR = sqrt(Transf.KX.^2+Transf.KY.^2);
func = @(n,u,v) v.^2./(u.^2+v.^2).*(v.*besselj(n,u).*besselk(n+1,v) - u.*besselj(n+1,u).*besselk(n,v));
VDk = -0.5*func(0,KR*rhoc,abs(Transf.KZ)*rhoc) + (Transf.KX.^2./KR.^2 - 0.5).*func(2,KR*rhoc,abs(Transf.KZ)*rhoc);
VDk = (1/3)*(3*(cos(alph).^2)-1) - VDk;
VDk(KR==0) = -1/3 + 1/2*abs(Transf.KZ(KR==0))*rhoc.*besselk(1,abs(Transf.KZ(KR==0))*rhoc);
VDk(Transf.KZ==0) = 1/6 + (Transf.KX(Transf.KZ==0).^2-Transf.KY(Transf.KZ==0).^2)./(KR(Transf.KZ==0).^2).*(1/2 - besselj(1,KR(Transf.KZ==0)*rhoc)./(KR(Transf.KZ==0)*rhoc));
VDk(1,1,1) = 1/6;
case 'Spherical' %Spherical
Rcut = min(Params.Lx/2,Params.Ly/2,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;
K = sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
VDk = (cos(alph).^2-1/3).*(1 + 3*cos(Rcut*K)./(Rcut^2.*K.^2) - 3*sin(Rcut*K)./(Rcut^3.*K.^3));
otherwise
disp('Choose a valid DDI cutoff type!')
return
end
end

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

@ -2,13 +2,12 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
properties (Access = public)
NumberOfAtoms;
DipolarAngle;
DipolarPolarAngle;
DipolarAzimuthAngle
ScatteringLength;
TrapFrequencies;
NumberOfPoints;
NumberOfGridPoints;
Dimensions;
CutoffType;
PotentialType;
SimulationMode;
TimeStep;
@ -31,13 +30,15 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
p.KeepUnmatched = true;
addParameter(p, 'NumberOfAtoms', 1E6,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 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;

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

@ -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);

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

@ -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;

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

@ -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