94 lines
3.5 KiB
Matlab
94 lines
3.5 KiB
Matlab
function [psi,V,VDk,gammaQF] = Initialize(Params,Transf)
|
|
|
|
format long
|
|
X = Transf.X; Y = Transf.Y; Z = Transf.Z;
|
|
m = Params.m;
|
|
|
|
Zcutoff = Params.Lz/2;
|
|
|
|
% == Potential == %
|
|
V = 0.5*(Params.gx.*X.^2+Params.gy.*Y.^2+Params.gz*Z.^2);
|
|
|
|
% == Calulating the DDIs == %
|
|
% For a cylindrical cutoff, we first construct a kr grid based on the 3D parameters using Bessel quadrature
|
|
loadDDI = 1;
|
|
|
|
if loadDDI == 1
|
|
VDk = load(sprintf('./Data/VDk_M.mat'));
|
|
VDk = VDk.VDk;
|
|
else
|
|
Params.Lr = 0.5*min(Params.Lx,Params.Ly);
|
|
Params.Nr = max(Params.Nx,Params.Ny);
|
|
[TransfRad] = setup_space_radial(Params); %morder really doesn't matter
|
|
VDk = VDcutoff(TransfRad.kr,TransfRad.kz,TransfRad.Rmax,Zcutoff);
|
|
|
|
disp('Calculated radial grid and cutoff')
|
|
|
|
% VDk = interp2(DDI.kz,DDI.kr,DDI.VDk,Transf.kz,Transf.kr,'spline');
|
|
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(VDk',2),VDk']';
|
|
VDk = interpn(KR,KZ,fullVDK,KR3D,KZ3D,'spline',-1/3); %Interpolating the radial VDk onto a new grid
|
|
VDk = fftshift(VDk);
|
|
save(sprintf('./Data/VDk_M.mat'),'VDk');
|
|
end
|
|
disp('Finished DDI')
|
|
|
|
|
|
%XXXX
|
|
% alph=acos((Transf.KZ*cos(Params.theta)+Transf.KX*sin(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));
|
|
% % func2 = @(n,u) (8*besselj(n,u) - 2*u.*besselj(n+1,u))./u.^2;
|
|
%
|
|
% 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;
|
|
%XXXX
|
|
|
|
|
|
% % == Setting up the initial wavefunction == %
|
|
loadstate = 0;
|
|
|
|
if loadstate == 1
|
|
loadnumber = 1250;
|
|
% load(sprintf('./Data/Seed/psi_%i.mat',loadnumber),'psi');
|
|
load(sprintf('./Data/Run_004/psi_gs.mat'),'psi');
|
|
|
|
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
|
|
psi = sqrt(Params.N)*psi/sqrt(Norm);
|
|
|
|
else
|
|
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 = 4*sqrt(2)*ellx; Ry = 4*sqrt(2)*elly; Rz = sqrt(2)*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);
|
|
|
|
psiz = exp(-(Z-Z0).^2/Rz^2)/sqrt(ellz*sqrt(pi));
|
|
psi2d = load(sprintf('./Data/Seed/psi_2d_SS.mat'),'psiseed_2d'); psi2d = psi2d.psiseed_2d;
|
|
psi = psiz.*repmat(psi2d,[1 1 length(Transf.z)]);
|
|
|
|
%add some noise
|
|
r = normrnd(0,1,size(X));
|
|
theta = rand(size(X));
|
|
noise = r.*exp(2*pi*1i*theta);
|
|
psi = psi + 0.00*noise;
|
|
|
|
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
|
|
Norm
|
|
psi = sqrt(Params.N)*psi/sqrt(Norm);
|
|
end |