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

latest commit

Browse Source
This commit is contained in:
Karthik 2024-06-13 02:02:44 +02:00
parent a1e0ad544e
commit 7cab430c05
14 changed files with 381 additions and 70 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
Dipolar Gas Simulator/+Helper/PhysicsConstants.m
View File

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

94
Dipolar Gas Simulator/+Plotter/visualizeDDIPotential.m Normal file
View File

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

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

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

95
Dipolar Gas Simulator/+Plotter/visualizeTrapPotential.m Normal file
View File

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

92
Dipolar Gas Simulator/+Plotter/visualizeWavefunction.m Normal file
View File

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

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

@ -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)
% Plotter.visualizeSpace(Transf)
Plotter.visualizeTrapPotential(V,Params,Transf)
% Plotter.visualizeWavefunction(psi,Params,Transf)

4
Dipolar Gas Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m
View File

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

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

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

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

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

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

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

3
Dipolar Gas Simulator/+Simulator/@DipolarGas/runSimulation.m
View File

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

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

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

12
Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m
View File

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

14
Dipolar Gas Simulator/+Simulator/@DipolarGas/setupWavefunction.m
View File

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