Latest working version.
This commit is contained in:
parent
570995f3f5
commit
6bb354de7a
@ -28,11 +28,10 @@ options = Helper.convertstruct2cell(OptionsStruct)
|
||||
clear OptionsStruct
|
||||
|
||||
sim = Simulator.DipolarGas(options{:});
|
||||
calc = Simulator.Calculator(options{:});
|
||||
pot = Simulator.Potentials(options{:});
|
||||
|
||||
%-% Run Simulation %-%
|
||||
[Params, Transf, psi, V, VDk] = sim.runSimulation(calc);
|
||||
[Params, Transf, psi, V, VDk] = sim.runSimulation();
|
||||
|
||||
%% - Plot numerical grid
|
||||
Plotter.visualizeSpace(Transf)
|
||||
|
||||
@ -27,13 +27,23 @@ 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;
|
||||
|
||||
p = inputParser;
|
||||
p.KeepUnmatched = true;
|
||||
addParameter(p, 'CutoffType', this.CalculatorDefaults.CutoffType,...
|
||||
@(x) any(strcmpi(x,{'Cylindrical','CylindricalInfiniteZ', 'Spherical'})));
|
||||
|
||||
p.parse(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 = p.Results.CutoffType;
|
||||
this.CalculatorDefaults.CutoffType = this.CutoffType;
|
||||
|
||||
end
|
||||
|
||||
function restoreDefaults(this)
|
||||
|
||||
@ -1,28 +1,29 @@
|
||||
function muchem = calculateChemicalPotential(psi,Params,Transf,VDk,V)
|
||||
function muchem = calculateChemicalPotential(~,psi,Params,Transf,VDk,V)
|
||||
|
||||
%Parameters
|
||||
normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
|
||||
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||
%Parameters
|
||||
normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
|
||||
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||
|
||||
% DDIs
|
||||
frho=fftn(abs(psi).^2);
|
||||
Phi=real(ifftn(frho.*VDk));
|
||||
% DDIs
|
||||
frho=fftn(abs(psi).^2);
|
||||
Phi=real(ifftn(frho.*VDk));
|
||||
|
||||
Eddi = (Params.gdd*Phi.*abs(psi).^2);
|
||||
Eddi = (Params.gdd*Phi.*abs(psi).^2);
|
||||
|
||||
%Kinetic energy
|
||||
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
|
||||
Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
|
||||
%Kinetic energy
|
||||
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
|
||||
Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
|
||||
|
||||
%Potential energy
|
||||
Epot = V.*abs(psi).^2;
|
||||
%Potential energy
|
||||
Epot = V.*abs(psi).^2;
|
||||
|
||||
%Contact interactions
|
||||
Eint = Params.gs*abs(psi).^4;
|
||||
%Contact interactions
|
||||
Eint = Params.gs*abs(psi).^4;
|
||||
|
||||
%Quantum fluctuations
|
||||
Eqf = Params.gammaQF*abs(psi).^5;
|
||||
%Quantum fluctuations
|
||||
Eqf = Params.gammaQF*abs(psi).^5;
|
||||
|
||||
%Total energy
|
||||
muchem = Ekin + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy*Transf.dz; %
|
||||
muchem = muchem / Params.N;
|
||||
%Total energy
|
||||
muchem = Ekin + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy*Transf.dz; %
|
||||
muchem = muchem / Params.N;
|
||||
end
|
||||
@ -1,35 +1,35 @@
|
||||
function E = calculateEnergyComponents(psi,Params,Transf,VDk,V)
|
||||
function E = calculateEnergyComponents(~,psi,Params,Transf,VDk,V)
|
||||
|
||||
%Parameters
|
||||
%Parameters
|
||||
|
||||
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||
normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
|
||||
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||
normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
|
||||
|
||||
% DDIs
|
||||
frho = fftn(abs(psi).^2);
|
||||
Phi = real(ifftn(frho.*VDk));
|
||||
% DDIs
|
||||
frho = fftn(abs(psi).^2);
|
||||
Phi = real(ifftn(frho.*VDk));
|
||||
|
||||
Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2;
|
||||
E.Eddi = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2;
|
||||
E.Eddi = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
|
||||
% EddiTot = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
% EddiTot = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
|
||||
%Kinetic energy
|
||||
% psik = ifftshift(fftn(fftshift(psi)))*normfac;
|
||||
%Kinetic energy
|
||||
% psik = ifftshift(fftn(fftshift(psi)))*normfac;
|
||||
|
||||
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
|
||||
E.Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
|
||||
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
|
||||
E.Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
|
||||
|
||||
% Potential energy
|
||||
Epot = V.*abs(psi).^2;
|
||||
E.Epot = trapz(Epot(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
% Potential energy
|
||||
Epot = V.*abs(psi).^2;
|
||||
E.Epot = trapz(Epot(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
|
||||
%Contact interactions
|
||||
Eint = 0.5*Params.gs*abs(psi).^4;
|
||||
E.Eint = trapz(Eint(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
%Contact interactions
|
||||
Eint = 0.5*Params.gs*abs(psi).^4;
|
||||
E.Eint = trapz(Eint(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
|
||||
%Quantum fluctuations
|
||||
Eqf = 0.4*Params.gammaQF*abs(psi).^5;
|
||||
E.Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
%Quantum fluctuations
|
||||
Eqf = 0.4*Params.gammaQF*abs(psi).^5;
|
||||
E.Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
|
||||
% plot(Transf.x,abs(psi(:,end/2,end/2+1)).^2)
|
||||
end
|
||||
|
||||
@ -1,24 +1,25 @@
|
||||
function res = calculateNormalizedResiduals(psi,Params,Transf,VDk,V,muchem)
|
||||
function res = calculateNormalizedResiduals(~,psi,Params,Transf,VDk,V,muchem)
|
||||
|
||||
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||
|
||||
% DDIs
|
||||
frho=fftn(abs(psi).^2);
|
||||
Phi=real(ifftn(frho.*VDk));
|
||||
% DDIs
|
||||
frho=fftn(abs(psi).^2);
|
||||
Phi=real(ifftn(frho.*VDk));
|
||||
|
||||
Eddi = Params.gdd*Phi.*psi;
|
||||
Eddi = Params.gdd*Phi.*psi;
|
||||
|
||||
%Kinetic energy
|
||||
Ekin = ifftn(KEop.*fftn(psi));
|
||||
%Kinetic energy
|
||||
Ekin = ifftn(KEop.*fftn(psi));
|
||||
|
||||
%Potential energy
|
||||
Epot = V.*psi;
|
||||
%Potential energy
|
||||
Epot = V.*psi;
|
||||
|
||||
%Contact interactions
|
||||
Eint = Params.gs*abs(psi).^2.*psi;
|
||||
%Contact interactions
|
||||
Eint = Params.gs*abs(psi).^2.*psi;
|
||||
|
||||
%Quantum fluctuations
|
||||
Eqf = Params.gammaQF*abs(psi).^3.*psi;
|
||||
%Quantum fluctuations
|
||||
Eqf = Params.gammaQF*abs(psi).^3.*psi;
|
||||
|
||||
%Total energy
|
||||
res = trapz(abs(Ekin(:) + Epot(:) + Eint(:) + Eddi(:) + Eqf(:) - muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz)/trapz(abs(muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz);
|
||||
%Total energy
|
||||
res = trapz(abs(Ekin(:) + Epot(:) + Eint(:) + Eddi(:) + Eqf(:) - muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz)/trapz(abs(muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz);
|
||||
end
|
||||
@ -1,39 +1,39 @@
|
||||
function VDkSemi = calculateNumericalHankelTransform(this,kr,kz,Rmax,Zmax,Nr)
|
||||
|
||||
% accuracy inputs for numerical integration
|
||||
if(nargin==5)
|
||||
Nr = 5e4;
|
||||
end
|
||||
|
||||
Nz = 64;
|
||||
farRmultiple = 2000;
|
||||
|
||||
% midpoint grids for the integration over 0<z<Zmax, Rmax<r<farRmultiple*Rmax (i.e. starts at Rmax)
|
||||
dr=(farRmultiple-1)*Rmax/Nr;
|
||||
r = ((1:Nr)'-0.5)*dr+Rmax;
|
||||
dz=Zmax/Nz;
|
||||
z = ((1:Nz)-0.5)*dz;
|
||||
[R, Z] = ndgrid(r,z);
|
||||
Rsq = R.^2 + Z.^2;
|
||||
|
||||
% real space interaction to be transformed
|
||||
igrandbase = (1 - 3*Z.^2./Rsq)./Rsq.^(3/2);
|
||||
|
||||
% do the Hankel/Fourier-Bessel transform numerically
|
||||
% prestore to ensure each besselj is only calculated once
|
||||
% cell is faster than (:,:,krn) slicing
|
||||
Nkr = numel(kr);
|
||||
besselr = cell(Nkr,1);
|
||||
for krn = 1:Nkr
|
||||
besselr{krn} = repmat(r.*besselj(0,kr(krn)*r),1,Nz);
|
||||
end
|
||||
|
||||
VDkSemi = zeros([Nkr,numel(kz)]);
|
||||
for kzn = 1:numel(kz)
|
||||
igrandbasez = repmat(cos(kz(kzn)*z),Nr,1) .* igrandbase;
|
||||
for krn = 1:Nkr
|
||||
igrand = igrandbasez.*besselr{krn};
|
||||
VDkSemi(krn,kzn) = VDkSemi(krn,kzn) - sum(igrand(:))*dz*dr;
|
||||
end
|
||||
function VDkSemi = calculateNumericalHankelTransform(~,kr,kz,Rmax,Zmax,Nr)
|
||||
|
||||
% accuracy inputs for numerical integration
|
||||
if(nargin==5)
|
||||
Nr = 5e4;
|
||||
end
|
||||
|
||||
Nz = 64;
|
||||
farRmultiple = 2000;
|
||||
|
||||
% midpoint grids for the integration over 0<z<Zmax, Rmax<r<farRmultiple*Rmax (i.e. starts at Rmax)
|
||||
dr=(farRmultiple-1)*Rmax/Nr;
|
||||
r = ((1:Nr)'-0.5)*dr+Rmax;
|
||||
dz=Zmax/Nz;
|
||||
z = ((1:Nz)-0.5)*dz;
|
||||
[R, Z] = ndgrid(r,z);
|
||||
Rsq = R.^2 + Z.^2;
|
||||
|
||||
% real space interaction to be transformed
|
||||
igrandbase = (1 - 3*Z.^2./Rsq)./Rsq.^(3/2);
|
||||
|
||||
% do the Hankel/Fourier-Bessel transform numerically
|
||||
% prestore to ensure each besselj is only calculated once
|
||||
% cell is faster than (:,:,krn) slicing
|
||||
Nkr = numel(kr);
|
||||
besselr = cell(Nkr,1);
|
||||
for krn = 1:Nkr
|
||||
besselr{krn} = repmat(r.*besselj(0,kr(krn)*r),1,Nz);
|
||||
end
|
||||
|
||||
VDkSemi = zeros([Nkr,numel(kz)]);
|
||||
for kzn = 1:numel(kz)
|
||||
igrandbasez = repmat(cos(kz(kzn)*z),Nr,1) .* igrandbase;
|
||||
for krn = 1:Nkr
|
||||
igrand = igrandbasez.*besselr{krn};
|
||||
VDkSemi(krn,kzn) = VDkSemi(krn,kzn) - sum(igrand(:))*dz*dr;
|
||||
end
|
||||
end
|
||||
end
|
||||
@ -1,4 +1,4 @@
|
||||
function [m_Order] = calculateOrderParameter(psi,Transf,Params,VDk,V,T,muchem)
|
||||
function [m_Order] = calculateOrderParameter(~,psi,Transf,Params,VDk,V,T,muchem)
|
||||
|
||||
NumRealiz = 100;
|
||||
|
||||
|
||||
@ -1,18 +1,19 @@
|
||||
function [PhaseC] = calculatePhaseCoherence(psi,Transf,Params)
|
||||
function [PhaseC] = calculatePhaseCoherence(~,psi,Transf,Params)
|
||||
|
||||
norm = sum(sum(sum(abs(psi).^2,1),2),3)*Transf.dx*Transf.dy*Transf.dz;
|
||||
psi = psi/sqrt(norm);
|
||||
norm = sum(sum(sum(abs(psi).^2,1),2),3)*Transf.dx*Transf.dy*Transf.dz;
|
||||
psi = psi/sqrt(norm);
|
||||
|
||||
NumGlobalShifts = 800;
|
||||
betas = []; phishift = [];
|
||||
for jj = 1:NumGlobalShifts
|
||||
phishift(jj) = -pi + 2*pi*(jj-1)/NumGlobalShifts;
|
||||
betas(jj) = sum(sum(sum(abs(angle(psi*exp(-1i*phishift(jj)))).*abs(psi).^2)));
|
||||
NumGlobalShifts = 800;
|
||||
betas = []; phishift = [];
|
||||
for jj = 1:NumGlobalShifts
|
||||
phishift(jj) = -pi + 2*pi*(jj-1)/NumGlobalShifts;
|
||||
betas(jj) = sum(sum(sum(abs(angle(psi*exp(-1i*phishift(jj)))).*abs(psi).^2)));
|
||||
end
|
||||
[minbeta,minidx] = min(betas);
|
||||
|
||||
psi = psi*exp(-1i*phishift(minidx));
|
||||
phi = angle(psi);
|
||||
|
||||
avgphi = sum(sum(sum(phi.*abs(psi).^2,1),2),3)*Transf.dx*Transf.dy*Transf.dz;
|
||||
PhaseC = sum(sum(sum(abs(angle(psi)-avgphi).*abs(psi).^2)))*Transf.dx*Transf.dy*Transf.dz;
|
||||
end
|
||||
[minbeta,minidx] = min(betas);
|
||||
|
||||
psi = psi*exp(-1i*phishift(minidx));
|
||||
phi = angle(psi);
|
||||
|
||||
avgphi = sum(sum(sum(phi.*abs(psi).^2,1),2),3)*Transf.dx*Transf.dy*Transf.dz;
|
||||
PhaseC = sum(sum(sum(abs(angle(psi)-avgphi).*abs(psi).^2)))*Transf.dx*Transf.dy*Transf.dz;
|
||||
@ -1,31 +1,32 @@
|
||||
function E = calculateTotalEnergy(psi,Params,Transf,VDk,V)
|
||||
function E = calculateTotalEnergy(~,psi,Params,Transf,VDk,V)
|
||||
|
||||
%Parameters
|
||||
%Parameters
|
||||
|
||||
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||
normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
|
||||
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||
normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
|
||||
|
||||
% DDIs
|
||||
frho = fftn(abs(psi).^2);
|
||||
Phi = real(ifftn(frho.*VDk));
|
||||
% DDIs
|
||||
frho = fftn(abs(psi).^2);
|
||||
Phi = real(ifftn(frho.*VDk));
|
||||
|
||||
Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2;
|
||||
Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2;
|
||||
|
||||
% EddiTot = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
% EddiTot = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
|
||||
%Kinetic energy
|
||||
% psik = ifftshift(fftn(fftshift(psi)))*normfac;
|
||||
%Kinetic energy
|
||||
% psik = ifftshift(fftn(fftshift(psi)))*normfac;
|
||||
|
||||
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
|
||||
Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
|
||||
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
|
||||
Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
|
||||
|
||||
% Potential energy
|
||||
Epot = V.*abs(psi).^2;
|
||||
% Potential energy
|
||||
Epot = V.*abs(psi).^2;
|
||||
|
||||
%Contact interactions
|
||||
Eint = 0.5*Params.gs*abs(psi).^4;
|
||||
%Contact interactions
|
||||
Eint = 0.5*Params.gs*abs(psi).^4;
|
||||
|
||||
%Quantum fluctuations
|
||||
Eqf = 0.4*Params.gammaQF*abs(psi).^5;
|
||||
%Quantum fluctuations
|
||||
Eqf = 0.4*Params.gammaQF*abs(psi).^5;
|
||||
|
||||
E = Ekin + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
E = Ekin + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
end
|
||||
@ -15,6 +15,8 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
|
||||
EnergyTolerance;
|
||||
MinimumTimeStep;
|
||||
|
||||
Calculator;
|
||||
|
||||
%Flags
|
||||
|
||||
DebugMode;
|
||||
@ -46,6 +48,8 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
|
||||
@(x) assert(isnumeric(x) && isvector(x) && all(x > 0)));
|
||||
addParameter(p, 'SimulationMode', 'ImaginaryTimeEvolution',...
|
||||
@(x) any(strcmpi(x,{'ImaginaryTimeEvolution','RealTimeEvolution'})));
|
||||
addParameter(p, 'CutoffType', 'Cylindrical',...
|
||||
@(x) any(strcmpi(x,{'Cylindrical','CylindricalInfiniteZ', 'Spherical'})));
|
||||
addParameter(p, 'TimeStep', 5E-4,...
|
||||
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
|
||||
addParameter(p, 'SimulationTime', 2e6,...
|
||||
@ -81,6 +85,8 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
|
||||
this.DoSave = p.Results.SaveData;
|
||||
this.SaveDirectory = p.Results.SaveDirectory;
|
||||
|
||||
this.Calculator = Simulator.Calculator('CutoffType', p.Results.CutoffType);
|
||||
|
||||
switch this.SimulationMode
|
||||
case "ImaginaryTimeEvolution"
|
||||
% Development In progress
|
||||
|
||||
@ -1,22 +1,20 @@
|
||||
function [psi,V,VDk] = initialize(this,calcObj,Params,Transf,TransfRad)
|
||||
function [psi,V,VDk] = initialize(this,Params,Transf,TransfRad)
|
||||
format long
|
||||
X = Transf.X; Y = Transf.Y; Z = Transf.Z;
|
||||
|
||||
format long
|
||||
X = Transf.X; Y = Transf.Y; Z = Transf.Z;
|
||||
% == Trap potential == %
|
||||
V = 0.5*(Params.gx.*X.^2+Params.gy.*Y.^2+Params.gz*Z.^2);
|
||||
|
||||
% == Trap potential == %
|
||||
V = 0.5*(Params.gx.*X.^2+Params.gy.*Y.^2+Params.gz*Z.^2);
|
||||
|
||||
% == Calculating the DDIs == %
|
||||
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,TransfRad);
|
||||
save(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat')),'VDk');
|
||||
end
|
||||
fprintf('Computed and saved DDI potential in Fourier space with %s cutoff.', calcObj.CutoffType)
|
||||
|
||||
% == Setting up the initial wavefunction == %
|
||||
psi = this.setupWavefunction(Params,Transf);
|
||||
% == Calculating the DDIs == %
|
||||
if isfile(strcat(this.SaveDirectory, '/VDk_M.mat'))
|
||||
VDk = load(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat')));
|
||||
VDk = VDk.VDk;
|
||||
else
|
||||
VDk = this.Calculator.calculateVDCutoff(Params,Transf,TransfRad);
|
||||
save(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat')),'VDk');
|
||||
end
|
||||
fprintf('Computed and saved DDI potential in Fourier space with %s cutoff.', this.Calculator.CutoffType)
|
||||
|
||||
% == Setting up the initial wavefunction == %
|
||||
psi = this.setupWavefunction(Params,Transf);
|
||||
end
|
||||
@ -1,4 +1,4 @@
|
||||
function [Params, Transf, psi,V,VDk] = runSimulation(this,calcObj)
|
||||
function [Params, Transf, psi,V,VDk] = runSimulation(this)
|
||||
% --- Obtain simulation parameters ---
|
||||
[Params] = this.setupParameters();
|
||||
|
||||
@ -9,15 +9,15 @@ function [Params, Transf, psi,V,VDk] = runSimulation(this,calcObj)
|
||||
% --- Initialize ---
|
||||
mkdir(sprintf(this.SaveDirectory))
|
||||
|
||||
[psi,V,VDk] = this.initialize(calcObj,Params,Transf,TransfRad);
|
||||
[psi,V,VDk] = this.initialize(Params,Transf,TransfRad);
|
||||
|
||||
Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = [];
|
||||
t_idx = 1; %Start at t = 0;
|
||||
Observ.res_idx = 1;
|
||||
|
||||
% --- Job Settings ---
|
||||
% njob = 6;
|
||||
% mkdir(sprintf('./Data/Run_%03i',njob))
|
||||
njob = 6;
|
||||
mkdir(sprintf('./Data/Run_%03i',njob))
|
||||
|
||||
% --- Run Simulation ---
|
||||
% [psi] = this.solver(psi,Params,Transf,VDk,V,njob,t_idx,Observ);
|
||||
|
||||
@ -1,92 +1,91 @@
|
||||
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.
|
||||
% 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);
|
||||
|
||||
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);
|
||||
% Dimensions (in units of l0)
|
||||
Params.Lx = this.Dimensions(1);
|
||||
Params.Ly = this.Dimensions(2);
|
||||
Params.Lz = this.Dimensions(3);
|
||||
|
||||
% Dimensions (in units of l0)
|
||||
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
|
||||
|
||||
% Masses
|
||||
Params.m = CONSTANTS.Dy164Mass;
|
||||
l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length
|
||||
% Atom numbers
|
||||
Params.N = this.NumberOfAtoms;
|
||||
|
||||
% Atom numbers
|
||||
Params.N = this.NumberOfAtoms;
|
||||
% Dipole angle
|
||||
Params.theta = this.DipolarPolarAngle; % pi/2 dipoles along x, theta=0 dipoles along z
|
||||
Params.phi = this.DipolarAzimuthAngle;
|
||||
|
||||
% Dipole angle
|
||||
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;
|
||||
|
||||
% Dipole lengths (units of muB)
|
||||
Params.mu = CONSTANTS.DyMagneticMoment;
|
||||
% Scattering lengths
|
||||
Params.as = this.ScatteringLength*a0;
|
||||
|
||||
% Scattering lengths
|
||||
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);
|
||||
|
||||
% Trapping frequencies
|
||||
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
|
||||
|
||||
% 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 = 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
|
||||
|
||||
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 ================ %
|
||||
|
||||
% ================ Parameters defined by those above ================ %
|
||||
% Contact interaction strength (units of l0/m)
|
||||
Params.gs = 4*pi*Params.as/l0;
|
||||
|
||||
% Contact interaction strength (units of l0/m)
|
||||
Params.gs = 4*pi*Params.as/l0;
|
||||
% Dipole lengths
|
||||
Params.add = mu0*Params.mu^2*Params.m/(12*pi*hbar^2);
|
||||
|
||||
% Dipole lengths
|
||||
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
|
||||
|
||||
% DDI strength
|
||||
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;
|
||||
|
||||
% Trap gamma
|
||||
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;
|
||||
if eps_dd == 0
|
||||
Q5 = 1;
|
||||
elseif eps_dd == 1
|
||||
Q5 = 3*sqrt(3)/2;
|
||||
else
|
||||
yeps = (1-eps_dd)/(3*eps_dd);
|
||||
Q5 = (3*eps_dd)^(5/2)*( (8+26*yeps+33*yeps^2)*sqrt(1+yeps) + 15*yeps^3*log((1+sqrt(1+yeps))/sqrt(yeps)) )/48;
|
||||
Q5 = real(Q5);
|
||||
end
|
||||
|
||||
% == Calculate LHY correction == %
|
||||
eps_dd = Params.add/Params.as;
|
||||
if eps_dd == 0
|
||||
Q5 = 1;
|
||||
elseif eps_dd == 1
|
||||
Q5 = 3*sqrt(3)/2;
|
||||
else
|
||||
yeps = (1-eps_dd)/(3*eps_dd);
|
||||
Q5 = (3*eps_dd)^(5/2)*( (8+26*yeps+33*yeps^2)*sqrt(1+yeps) + 15*yeps^3*log((1+sqrt(1+yeps))/sqrt(yeps)) )/48;
|
||||
Q5 = real(Q5);
|
||||
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.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;
|
||||
end
|
||||
@ -1,33 +1,33 @@
|
||||
function [Transf] = setupSpace(this,Params)
|
||||
Transf.Xmax = 0.5*Params.Lx;
|
||||
Transf.Ymax = 0.5*Params.Ly;
|
||||
Transf.Zmax = 0.5*Params.Lz;
|
||||
Transf.Xmax = 0.5*Params.Lx;
|
||||
Transf.Ymax = 0.5*Params.Ly;
|
||||
Transf.Zmax = 0.5*Params.Lz;
|
||||
|
||||
Nz = Params.Nz; Nx = Params.Nx; Ny = Params.Ny;
|
||||
Nz = Params.Nz; Nx = Params.Nx; Ny = Params.Ny;
|
||||
|
||||
% Fourier grids
|
||||
x = linspace(-0.5*Params.Lx,0.5*Params.Lx-Params.Lx/Params.Nx,Params.Nx);
|
||||
Kmax = pi*Params.Nx/Params.Lx;
|
||||
kx = linspace(-Kmax,Kmax,Nx+1);
|
||||
kx = kx(1:end-1); dkx = kx(2)-kx(1);
|
||||
kx = fftshift(kx);
|
||||
% Fourier grids
|
||||
x = linspace(-0.5*Params.Lx,0.5*Params.Lx-Params.Lx/Params.Nx,Params.Nx);
|
||||
Kmax = pi*Params.Nx/Params.Lx;
|
||||
kx = linspace(-Kmax,Kmax,Nx+1);
|
||||
kx = kx(1:end-1); dkx = kx(2)-kx(1);
|
||||
kx = fftshift(kx);
|
||||
|
||||
y = linspace(-0.5*Params.Ly,0.5*Params.Ly-Params.Ly/Params.Ny,Params.Ny);
|
||||
Kmax = pi*Params.Ny/Params.Ly;
|
||||
ky = linspace(-Kmax,Kmax,Ny+1);
|
||||
ky = ky(1:end-1); dky = ky(2)-ky(1);
|
||||
ky = fftshift(ky);
|
||||
y = linspace(-0.5*Params.Ly,0.5*Params.Ly-Params.Ly/Params.Ny,Params.Ny);
|
||||
Kmax = pi*Params.Ny/Params.Ly;
|
||||
ky = linspace(-Kmax,Kmax,Ny+1);
|
||||
ky = ky(1:end-1); dky = ky(2)-ky(1);
|
||||
ky = fftshift(ky);
|
||||
|
||||
z = linspace(-0.5*Params.Lz,0.5*Params.Lz-Params.Lz/Params.Nz,Params.Nz);
|
||||
Kmax = pi*Params.Nz/Params.Lz;
|
||||
kz = linspace(-Kmax,Kmax,Nz+1);
|
||||
kz = kz(1:end-1); dkz = kz(2)-kz(1);
|
||||
kz = fftshift(kz);
|
||||
z = linspace(-0.5*Params.Lz,0.5*Params.Lz-Params.Lz/Params.Nz,Params.Nz);
|
||||
Kmax = pi*Params.Nz/Params.Lz;
|
||||
kz = linspace(-Kmax,Kmax,Nz+1);
|
||||
kz = kz(1:end-1); dkz = kz(2)-kz(1);
|
||||
kz = fftshift(kz);
|
||||
|
||||
[Transf.X,Transf.Y,Transf.Z]=ndgrid(x,y,z);
|
||||
[Transf.KX,Transf.KY,Transf.KZ]=ndgrid(kx,ky,kz);
|
||||
Transf.x = x; Transf.y = y; Transf.z = z;
|
||||
Transf.kx = kx; Transf.ky = ky; Transf.kz = kz;
|
||||
Transf.dx = x(2)-x(1); Transf.dy = y(2)-y(1); Transf.dz = z(2)-z(1);
|
||||
Transf.dkx = dkx; Transf.dky = dky; Transf.dkz = dkz;
|
||||
[Transf.X,Transf.Y,Transf.Z]=ndgrid(x,y,z);
|
||||
[Transf.KX,Transf.KY,Transf.KZ]=ndgrid(kx,ky,kz);
|
||||
Transf.x = x; Transf.y = y; Transf.z = z;
|
||||
Transf.kx = kx; Transf.ky = ky; Transf.kz = kz;
|
||||
Transf.dx = x(2)-x(1); Transf.dy = y(2)-y(1); Transf.dz = z(2)-z(1);
|
||||
Transf.dkx = dkx; Transf.dky = dky; Transf.dkz = dkz;
|
||||
end
|
||||
@ -1,50 +1,49 @@
|
||||
function [Transf] = setupSpaceRadial(this,Params,morder)
|
||||
function [Transf] = setupSpaceRadial(~,Params,morder)
|
||||
|
||||
Params.Lr = 0.5*min(Params.Lx,Params.Ly);
|
||||
Params.Nr = max(Params.Nx,Params.Ny);
|
||||
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;
|
||||
Zmax = 0.5*Params.Lz;
|
||||
Rmax = Params.Lr;
|
||||
Nz = Params.Nz;
|
||||
Nr = Params.Nr;
|
||||
|
||||
if(nargin==2)
|
||||
morder=0; %only do Bessel J0
|
||||
end
|
||||
if(nargin==2)
|
||||
morder=0; %only do Bessel J0
|
||||
end
|
||||
|
||||
% Fourier grids
|
||||
z=linspace(-Zmax,Zmax,Nz+1);
|
||||
z=z(1:end-1);
|
||||
dz=z(2)-z(1);
|
||||
Kmax=Nz*2*pi/(4*Zmax);
|
||||
kz=linspace(-Kmax,Kmax,Nz+1);
|
||||
kz=kz(1:end-1);
|
||||
% Fourier grids
|
||||
z=linspace(-Zmax,Zmax,Nz+1);
|
||||
z=z(1:end-1);
|
||||
dz=z(2)-z(1);
|
||||
Kmax=Nz*2*pi/(4*Zmax);
|
||||
kz=linspace(-Kmax,Kmax,Nz+1);
|
||||
kz=kz(1:end-1);
|
||||
|
||||
% Hankel grids and transform
|
||||
H = hankelmatrix(morder,Rmax,Nr);
|
||||
r=H.r(:);
|
||||
kr=H.kr(:);
|
||||
T = diag(H.J/H.kmax)*H.T*diag(Rmax./H.J)*dz*(2*pi);
|
||||
Tinv = diag(H.J./Rmax)*H.T'*diag(H.kmax./H.J)/dz/(2*pi);
|
||||
wr=H.wr;
|
||||
wk=H.wk;
|
||||
% H.T'*diag(H.J/H.vmax)*H.T*diag(Rmax./H.J)
|
||||
% Hankel grids and transform
|
||||
H = hankelmatrix(morder,Rmax,Nr);
|
||||
r=H.r(:);
|
||||
kr=H.kr(:);
|
||||
T = diag(H.J/H.kmax)*H.T*diag(Rmax./H.J)*dz*(2*pi);
|
||||
Tinv = diag(H.J./Rmax)*H.T'*diag(H.kmax./H.J)/dz/(2*pi);
|
||||
wr=H.wr;
|
||||
wk=H.wk;
|
||||
% H.T'*diag(H.J/H.vmax)*H.T*diag(Rmax./H.J)
|
||||
|
||||
[Transf.R,Transf.Z]=ndgrid(r,z);
|
||||
[Transf.KR,Transf.KZ]=ndgrid(kr,kz);
|
||||
Transf.T=T;
|
||||
Transf.Tinv=Tinv;
|
||||
Transf.r=r;
|
||||
Transf.kr=kr;
|
||||
Transf.z=z;
|
||||
Transf.kz=kz;
|
||||
Transf.wr=wr;
|
||||
Transf.wk=wk;
|
||||
Transf.Rmax=Rmax;
|
||||
Transf.Zmax=Zmax;
|
||||
Transf.dz=z(2)-z(1);
|
||||
Transf.dkz=kz(2)-kz(1);
|
||||
%b1=Transf;
|
||||
[Transf.R,Transf.Z]=ndgrid(r,z);
|
||||
[Transf.KR,Transf.KZ]=ndgrid(kr,kz);
|
||||
Transf.T=T;
|
||||
Transf.Tinv=Tinv;
|
||||
Transf.r=r;
|
||||
Transf.kr=kr;
|
||||
Transf.z=z;
|
||||
Transf.kz=kz;
|
||||
Transf.wr=wr;
|
||||
Transf.wk=wk;
|
||||
Transf.Rmax=Rmax;
|
||||
Transf.Zmax=Zmax;
|
||||
Transf.dz=z(2)-z(1);
|
||||
Transf.dkz=kz(2)-kz(1);
|
||||
|
||||
function s_HT = hankelmatrix(order,rmax,Nr,eps_roots)
|
||||
%HANKEL_MATRIX: Generates data to use for Hankel Transforms
|
||||
@ -103,9 +102,9 @@ function s_HT = hankelmatrix(order,rmax,Nr,eps_roots)
|
||||
% See also bessel_zeros, besselj
|
||||
|
||||
if (~exist('eps_roots', 'var')||isemtpy(eps_roots))
|
||||
s_HT.eps_roots = 1e-6;
|
||||
s_HT.eps_roots = 1e-6;
|
||||
else
|
||||
s_HT.eps_roots = eps_roots;
|
||||
s_HT.eps_roots = eps_roots;
|
||||
end
|
||||
|
||||
s_HT.order = order;
|
||||
@ -166,146 +165,146 @@ mu3 = mu^3;
|
||||
mu4 = mu^4;
|
||||
|
||||
if (d<3)
|
||||
p = 7*mu - 31;
|
||||
p0 = mu - 1;
|
||||
if ((1+p)==p)
|
||||
p1 = 0;
|
||||
q1 = 0;
|
||||
else
|
||||
p1 = 4*(253*mu2 - 3722*mu+17869)*p0/(15*p);
|
||||
q1 = 1.6*(83*mu2 - 982*mu + 3779)/p;
|
||||
end
|
||||
p = 7*mu - 31;
|
||||
p0 = mu - 1;
|
||||
if ((1+p)==p)
|
||||
p1 = 0;
|
||||
q1 = 0;
|
||||
else
|
||||
p1 = 4*(253*mu2 - 3722*mu+17869)*p0/(15*p);
|
||||
q1 = 1.6*(83*mu2 - 982*mu + 3779)/p;
|
||||
end
|
||||
else
|
||||
p = 7*mu2 + 82*mu - 9;
|
||||
p0 = mu + 3;
|
||||
if ((p+1)==1)
|
||||
p1 = 0;
|
||||
q1 = 0;
|
||||
else
|
||||
p1 = (4048*mu4 + 131264*mu3 - 221984*mu2 - 417600*mu + 1012176)/(60*p);
|
||||
q1 = 1.6*(83*mu3 + 2075*mu2 - 3039*mu + 3537)/p;
|
||||
end
|
||||
p = 7*mu2 + 82*mu - 9;
|
||||
p0 = mu + 3;
|
||||
if ((p+1)==1)
|
||||
p1 = 0;
|
||||
q1 = 0;
|
||||
else
|
||||
p1 = (4048*mu4 + 131264*mu3 - 221984*mu2 - 417600*mu + 1012176)/(60*p);
|
||||
q1 = 1.6*(83*mu3 + 2075*mu2 - 3039*mu + 3537)/p;
|
||||
end
|
||||
end
|
||||
|
||||
if (d==1)|(d==4)
|
||||
t = .25;
|
||||
t = .25;
|
||||
else
|
||||
t = .75;
|
||||
t = .75;
|
||||
end
|
||||
tt = 4*t;
|
||||
|
||||
if (d<3)
|
||||
pp1 = 5/48;
|
||||
qq1 = -5/36;
|
||||
pp1 = 5/48;
|
||||
qq1 = -5/36;
|
||||
else
|
||||
pp1 = -7/48;
|
||||
qq1 = 35/288;
|
||||
pp1 = -7/48;
|
||||
qq1 = 35/288;
|
||||
end
|
||||
|
||||
y = .375*pi;
|
||||
if (a>=3)
|
||||
bb = a^(-2/3);
|
||||
bb = a^(-2/3);
|
||||
else
|
||||
bb = 1;
|
||||
bb = 1;
|
||||
end
|
||||
a1 = 3*a - 8;
|
||||
% psi = (.5*a + .25)*pi;
|
||||
|
||||
for s=1:n
|
||||
if ((a==0)&(s==1)&(d==3))
|
||||
x = 0;
|
||||
j = 0;
|
||||
else
|
||||
if (s>=a1)
|
||||
b = (s + .5*a - t)*pi;
|
||||
c = .015625/(b^2);
|
||||
x = b - .125*(p0 - p1*c)/(b*(1 - q1*c));
|
||||
else
|
||||
if (s==1)
|
||||
switch (d)
|
||||
case (1)
|
||||
x = -2.33811;
|
||||
case (2)
|
||||
x = -1.17371;
|
||||
case (3)
|
||||
x = -1.01879;
|
||||
otherwise
|
||||
x = -2.29444;
|
||||
end
|
||||
else
|
||||
x = y*(4*s - tt);
|
||||
v = x^(-2);
|
||||
x = -x^(2/3) * (1 + v*(pp1 + qq1*v));
|
||||
end
|
||||
u = x*bb;
|
||||
v = fi(2/3 * (-u)^1.5);
|
||||
w = 1/cos(v);
|
||||
xx = 1 - w^2;
|
||||
c = sqrt(u/xx);
|
||||
if (d<3)
|
||||
x = w*(a + c*(-5/u - c*(6 - 10/xx))/(48*a*u));
|
||||
else
|
||||
x = w*(a + c*(7/u + c*(18 - 14/xx))/(48*a*u));
|
||||
end
|
||||
end
|
||||
j = 0;
|
||||
if ((a==0)&(s==1)&(d==3))
|
||||
x = 0;
|
||||
j = 0;
|
||||
else
|
||||
if (s>=a1)
|
||||
b = (s + .5*a - t)*pi;
|
||||
c = .015625/(b^2);
|
||||
x = b - .125*(p0 - p1*c)/(b*(1 - q1*c));
|
||||
else
|
||||
if (s==1)
|
||||
switch (d)
|
||||
case (1)
|
||||
x = -2.33811;
|
||||
case (2)
|
||||
x = -1.17371;
|
||||
case (3)
|
||||
x = -1.01879;
|
||||
otherwise
|
||||
x = -2.29444;
|
||||
end
|
||||
else
|
||||
x = y*(4*s - tt);
|
||||
v = x^(-2);
|
||||
x = -x^(2/3) * (1 + v*(pp1 + qq1*v));
|
||||
end
|
||||
u = x*bb;
|
||||
v = fi(2/3 * (-u)^1.5);
|
||||
w = 1/cos(v);
|
||||
xx = 1 - w^2;
|
||||
c = sqrt(u/xx);
|
||||
if (d<3)
|
||||
x = w*(a + c*(-5/u - c*(6 - 10/xx))/(48*a*u));
|
||||
else
|
||||
x = w*(a + c*(7/u + c*(18 - 14/xx))/(48*a*u));
|
||||
end
|
||||
end
|
||||
j = 0;
|
||||
|
||||
while ((j==0)|((j<5)&(abs(w/x)>e)))
|
||||
xx = x^2;
|
||||
x4 = x^4;
|
||||
a2 = aa - xx;
|
||||
r0 = bessr(d, a, x);
|
||||
j = j+1;
|
||||
if (d<3)
|
||||
u = r0;
|
||||
w = 6*x*(2*a + 1);
|
||||
p = (1 - 4*a2)/w;
|
||||
q = (4*(xx-mu) - 2 - 12*a)/w;
|
||||
else
|
||||
u = -xx*r0/a2;
|
||||
v = 2*x*a2/(3*(aa+xx));
|
||||
w = 64*a2^3;
|
||||
q = 2*v*(1 + mu2 + 32*mu*xx + 48*x4)/w;
|
||||
p = v*(1 + (40*mu*xx + 48*x4 - mu2)/w);
|
||||
end
|
||||
w = u*(1 + p*r0)/(1 + q*r0);
|
||||
x = x+w;
|
||||
end
|
||||
z(s) = x;
|
||||
end
|
||||
xx = x^2;
|
||||
x4 = x^4;
|
||||
a2 = aa - xx;
|
||||
r0 = bessr(d, a, x);
|
||||
j = j+1;
|
||||
if (d<3)
|
||||
u = r0;
|
||||
w = 6*x*(2*a + 1);
|
||||
p = (1 - 4*a2)/w;
|
||||
q = (4*(xx-mu) - 2 - 12*a)/w;
|
||||
else
|
||||
u = -xx*r0/a2;
|
||||
v = 2*x*a2/(3*(aa+xx));
|
||||
w = 64*a2^3;
|
||||
q = 2*v*(1 + mu2 + 32*mu*xx + 48*x4)/w;
|
||||
p = v*(1 + (40*mu*xx + 48*x4 - mu2)/w);
|
||||
end
|
||||
w = u*(1 + p*r0)/(1 + q*r0);
|
||||
x = x+w;
|
||||
end
|
||||
z(s) = x;
|
||||
end
|
||||
end
|
||||
|
||||
function FI = fi(y)
|
||||
c1 = 1.570796;
|
||||
if (~y)
|
||||
FI = 0;
|
||||
elseif (y>1e5)
|
||||
FI = c1;
|
||||
else
|
||||
if (y<1)
|
||||
p = (3*y)^(1/3);
|
||||
pp = p^2;
|
||||
p = p*(1 + pp*(pp*(27 - 2*pp) - 210)/1575);
|
||||
else
|
||||
p = 1/(y + c1);
|
||||
pp = p^2;
|
||||
p = c1 - p*(1 + pp*(2310 + pp*(3003 + pp*(4818 + pp*(8591 + pp*16328))))/3465);
|
||||
end
|
||||
pp = (y+p)^2;
|
||||
r = (p - atan(p+y))/pp;
|
||||
FI = p - (1+pp)*r*(1 + r/(p+y));
|
||||
end
|
||||
c1 = 1.570796;
|
||||
if (~y)
|
||||
FI = 0;
|
||||
elseif (y>1e5)
|
||||
FI = c1;
|
||||
else
|
||||
if (y<1)
|
||||
p = (3*y)^(1/3);
|
||||
pp = p^2;
|
||||
p = p*(1 + pp*(pp*(27 - 2*pp) - 210)/1575);
|
||||
else
|
||||
p = 1/(y + c1);
|
||||
pp = p^2;
|
||||
p = c1 - p*(1 + pp*(2310 + pp*(3003 + pp*(4818 + pp*(8591 + pp*16328))))/3465);
|
||||
end
|
||||
pp = (y+p)^2;
|
||||
r = (p - atan(p+y))/pp;
|
||||
FI = p - (1+pp)*r*(1 + r/(p+y));
|
||||
end
|
||||
return
|
||||
|
||||
function Jr = bessr(d,a,x)
|
||||
switch (d)
|
||||
case (1)
|
||||
Jr = besselj(a, x)./besselj(a+1, x);
|
||||
case (2)
|
||||
Jr = bessely(a, x)./bessely(a+1, x);
|
||||
case (3)
|
||||
Jr = a./x - besselj(a+1, x)./besselj(a, x);
|
||||
otherwise
|
||||
Jr = a./x - bessely(a+1, x)./bessely(a, x);
|
||||
end
|
||||
switch (d)
|
||||
case (1)
|
||||
Jr = besselj(a, x)./besselj(a+1, x);
|
||||
case (2)
|
||||
Jr = bessely(a, x)./bessely(a+1, x);
|
||||
case (3)
|
||||
Jr = a./x - besselj(a+1, x)./besselj(a, x);
|
||||
otherwise
|
||||
Jr = a./x - bessely(a+1, x)./bessely(a, x);
|
||||
end
|
||||
return
|
||||
@ -1,25 +1,25 @@
|
||||
function [psi] = setupWavefunction(this,Params,Transf)
|
||||
function [psi] = setupWavefunction(~,Params,Transf)
|
||||
|
||||
X = Transf.X; Y = Transf.Y; Z = Transf.Z;
|
||||
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;
|
||||
ellz = sqrt(Params.hbar/(Params.m*Params.wz))/Params.l0;
|
||||
ellx = sqrt(Params.hbar/(Params.m*Params.wx))/Params.l0;
|
||||
elly = sqrt(Params.hbar/(Params.m*Params.wy))/Params.l0;
|
||||
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;
|
||||
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);
|
||||
cur_norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
|
||||
psi = psi/sqrt(cur_norm);
|
||||
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);
|
||||
|
||||
% add some noise
|
||||
r = normrnd(0,1,size(X));
|
||||
theta = rand(size(X));
|
||||
noise = r.*exp(2*pi*1i*theta);
|
||||
psi = psi + 0.01*noise;
|
||||
% add some noise
|
||||
r = normrnd(0,1,size(X));
|
||||
theta = rand(size(X));
|
||||
noise = r.*exp(2*pi*1i*theta);
|
||||
psi = psi + 0.01*noise;
|
||||
|
||||
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
|
||||
psi = sqrt(Params.N)*psi/sqrt(Norm);
|
||||
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
|
||||
psi = sqrt(Params.N)*psi/sqrt(Norm);
|
||||
|
||||
end
|
||||
@ -1,90 +1,89 @@
|
||||
function [psi] = solver(this,psi,Params,Transf,VDk,V,njob,t_idx,Observ)
|
||||
function [psi] = propagate(this,psi,Params,Transf,VDk,V,njob,t_idx,Observ)
|
||||
set(0,'defaulttextInterpreter','latex')
|
||||
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
|
||||
|
||||
set(0,'defaulttextInterpreter','latex')
|
||||
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
|
||||
dt=-1j*abs(this.TimeStep);
|
||||
|
||||
dt=-1j*abs(this.TimeStep);
|
||||
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||
Observ.residual = 1; Observ.res = 1;
|
||||
|
||||
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||
Observ.residual = 1; Observ.res = 1;
|
||||
muchem = this.Calculator.ChemicalPotential(psi,Params,Transf,VDk,V);
|
||||
AdaptIdx = 0;
|
||||
|
||||
muchem = Simulator.ChemicalPotential(psi,Params,Transf,VDk,V);
|
||||
AdaptIdx = 0;
|
||||
while t_idx < Params.cut_off
|
||||
%kin
|
||||
psi = fftn(psi);
|
||||
psi = psi.*exp(-0.5*1i*dt*KEop);
|
||||
psi = ifftn(psi);
|
||||
|
||||
while t_idx < Params.cut_off
|
||||
%kin
|
||||
psi = fftn(psi);
|
||||
psi = psi.*exp(-0.5*1i*dt*KEop);
|
||||
psi = ifftn(psi);
|
||||
%DDI
|
||||
frho = fftn(abs(psi).^2);
|
||||
Phi = real(ifftn(frho.*VDk));
|
||||
|
||||
%DDI
|
||||
frho = fftn(abs(psi).^2);
|
||||
Phi = real(ifftn(frho.*VDk));
|
||||
%Real-space
|
||||
psi = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gdd*Phi + Params.gammaQF*abs(psi).^3 - muchem));
|
||||
|
||||
%Real-space
|
||||
psi = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gammaQF*abs(psi).^3 + Params.gdd*Phi - muchem));
|
||||
%kin
|
||||
psi = fftn(psi);
|
||||
psi = psi.*exp(-0.5*1i*dt*KEop);
|
||||
psi = ifftn(psi);
|
||||
|
||||
%kin
|
||||
psi = fftn(psi);
|
||||
psi = psi.*exp(-0.5*1i*dt*KEop);
|
||||
psi = ifftn(psi);
|
||||
%Renorm
|
||||
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
|
||||
psi = sqrt(Params.N)*psi/sqrt(Norm);
|
||||
|
||||
%Renorm
|
||||
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
|
||||
psi = sqrt(Params.N)*psi/sqrt(Norm);
|
||||
muchem = this.Calculator.ChemicalPotential(psi,Params,Transf,VDk,V);
|
||||
|
||||
muchem = Simulator.ChemicalPotential(psi,Params,Transf,VDk,V);
|
||||
if mod(t_idx,1000) == 0
|
||||
|
||||
if mod(t_idx,1000) == 0
|
||||
%Change in Energy
|
||||
E = this.Calculator.TotalEnergy(psi,Params,Transf,VDk,V);
|
||||
E = E/Norm;
|
||||
Observ.EVec = [Observ.EVec E];
|
||||
|
||||
%Change in Energy
|
||||
E = Simulator.TotalEnergy(psi,Params,Transf,VDk,V);
|
||||
E = E/Norm;
|
||||
Observ.EVec = [Observ.EVec E];
|
||||
%Chemical potential
|
||||
Observ.mucVec = [Observ.mucVec muchem];
|
||||
|
||||
%Chemical potential
|
||||
Observ.mucVec = [Observ.mucVec muchem];
|
||||
%Normalized residuals
|
||||
res = this.Calculator.NormalizedResiduals(psi,Params,Transf,VDk,V,muchem);
|
||||
Observ.residual = [Observ.residual res];
|
||||
|
||||
%Normalized residuals
|
||||
res = Simulator.NormalizedResiduals(psi,Params,Transf,VDk,V,muchem);
|
||||
Observ.residual = [Observ.residual res];
|
||||
Observ.res_idx = Observ.res_idx + 1;
|
||||
|
||||
Observ.res_idx = Observ.res_idx + 1;
|
||||
save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
|
||||
|
||||
save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
|
||||
|
||||
%Adaptive time step -- Careful, this can quickly get out of control
|
||||
relres = abs(Observ.residual(Observ.res_idx)-Observ.residual(Observ.res_idx-1))/Observ.residual(Observ.res_idx);
|
||||
if relres <1e-5
|
||||
if AdaptIdx > 4 && abs(dt) > Params.mindt
|
||||
dt = dt / 2;
|
||||
fprintf('Time step changed to '); disp(dt);
|
||||
AdaptIdx = 0;
|
||||
elseif AdaptIdx > 4 && abs(dt) < Params.mindt
|
||||
break
|
||||
%Adaptive time step -- Careful, this can quickly get out of control
|
||||
relres = abs(Observ.residual(Observ.res_idx)-Observ.residual(Observ.res_idx-1))/Observ.residual(Observ.res_idx);
|
||||
if relres <1e-5
|
||||
if AdaptIdx > 4 && abs(dt) > Params.mindt
|
||||
dt = dt / 2;
|
||||
fprintf('Time step changed to '); disp(dt);
|
||||
AdaptIdx = 0;
|
||||
elseif AdaptIdx > 4 && abs(dt) < Params.mindt
|
||||
break
|
||||
else
|
||||
AdaptIdx = AdaptIdx + 1;
|
||||
end
|
||||
else
|
||||
AdaptIdx = AdaptIdx + 1;
|
||||
AdaptIdx = 0;
|
||||
end
|
||||
else
|
||||
AdaptIdx = 0;
|
||||
end
|
||||
if any(isnan(psi(:)))
|
||||
disp('NaNs encountered!')
|
||||
break
|
||||
end
|
||||
t_idx=t_idx+1;
|
||||
end
|
||||
if any(isnan(psi(:)))
|
||||
disp('NaNs encountered!')
|
||||
break
|
||||
end
|
||||
t_idx=t_idx+1;
|
||||
end
|
||||
|
||||
%Change in Energy
|
||||
E = Simulator.TotalEnergy(psi,Params,Transf,VDk,V);
|
||||
E = E/Norm;
|
||||
Observ.EVec = [Observ.EVec E];
|
||||
|
||||
% Phase coherence
|
||||
[PhaseC] = Simulator.PhaseCoherence(psi,Transf,Params);
|
||||
Observ.PCVec = [Observ.PCVec PhaseC];
|
||||
|
||||
Observ.res_idx = Observ.res_idx + 1;
|
||||
save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
|
||||
|
||||
%Change in Energy
|
||||
E = this.Calculator.TotalEnergy(psi,Params,Transf,VDk,V);
|
||||
E = E/Norm;
|
||||
Observ.EVec = [Observ.EVec E];
|
||||
|
||||
% Phase coherence
|
||||
[PhaseC] = this.Calculator.PhaseCoherence(psi,Transf,Params);
|
||||
Observ.PCVec = [Observ.PCVec PhaseC];
|
||||
|
||||
Observ.res_idx = Observ.res_idx + 1;
|
||||
save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
|
||||
end
|
||||
Loading…
Reference in New Issue
Block a user