71 lines
1.8 KiB
Matlab
71 lines
1.8 KiB
Matlab
function [evals, modes] = solveBogoliubovdeGennesIn2D(psi, Params, VDk, VParams, Transf, muchem)
|
|
|
|
wz_tilde = Params.wz / Params.w0;
|
|
gs = Params.gs;
|
|
gdd = Params.gdd;
|
|
gammaQF = Params.gammaQF;
|
|
|
|
KEop = 0.5*(Transf.KX.^2+Transf.KY.^2);
|
|
Ez = (0.25*VParams.sigma^2) + (0.25*wz_tilde^2*VParams.sigma^2);
|
|
muchem_tilde = muchem - Ez;
|
|
|
|
g_pf_2D = 1/(sqrt(2*pi)*VParams.sigma);
|
|
gQF_pf_2D = sqrt(2/5)/(pi^(3/4)*VParams.sigma^(3/2));
|
|
|
|
% eigs only works with column vectors
|
|
psi = psi.';
|
|
KEop = KEop.';
|
|
VDk = VDk.';
|
|
|
|
% Interaction Potential
|
|
frho = fftn(abs(psi).^2);
|
|
Phi = real(ifftn(frho.*VDk));
|
|
|
|
% Operators
|
|
H = @(w) real(ifft(KEop.*fft(w)));
|
|
C = @(w) (((g_pf_2D*gs*abs(psi).^2) + (g_pf_2D*gdd*Phi)).*w) + (gQF_pf_2D*gammaQF*abs(psi).^3.*w);
|
|
muHC = @(w) (-muchem_tilde * w) + H(w) + C(w);
|
|
X = @(w,psi) (psi.*real(ifft(VDk.*fft(psi.*w)))) + (3/2)*(gQF_pf_2D*gammaQF*abs(psi).^3).*w;
|
|
|
|
% Operate in order on g
|
|
BdG = @(g) muHC(muHC(g) + (2.*X(g)));
|
|
|
|
syssize = size(psi);
|
|
opts.v0 = psi(:);
|
|
opts.tol = 1e-16;
|
|
opts.disp = 1;
|
|
opts.issym = 0;
|
|
opts.isreal = 1;
|
|
opts.maxit = 1e4;
|
|
Neigs = syssize;
|
|
|
|
[g,D] = eigs(BdG,syssize,Neigs,'sr',opts);
|
|
evals = diag(D);
|
|
clear D;
|
|
|
|
% Eigenvalues
|
|
evals = sqrt(evals);
|
|
|
|
% Obtain f from g
|
|
f = zeros(size(g));
|
|
for ii = 1:Neigs
|
|
f(:,:,ii) = (1/evals(ii)) * (muHC(g(:,:,ii)) + (2.*X(g(:,:,ii))));
|
|
end
|
|
|
|
% Obtain u and v from f and g
|
|
u = (f + g)/2;
|
|
v = (f - g)/2;
|
|
|
|
% Renormalize to \int |u|^2 - |v|^2 = 1
|
|
for ii=1:Neigs
|
|
normalization = sum(abs(u(:,:,ii)).^2 - abs(v(:,:,ii)).^2);
|
|
u(:,:,ii) = u(:,:,ii)/sqrt(normalization);
|
|
v(:,:,ii) = v(:,:,ii)/sqrt(normalization);
|
|
end
|
|
|
|
modes.u = u'; modes.v = v';
|
|
modes.g = g'; modes.f = f';
|
|
|
|
end
|
|
|