32 lines
1.3 KiB
Mathematica
32 lines
1.3 KiB
Mathematica
|
function muchem = calculateChemicalPotential(~,psi,Params,VParams,Transf,VDk,V)
|
||
|
|
|
||
|
|
g_eff = Params.gs*VParams.nu/(2^(1+1/VParams.nu)*VParams.ell*gamma(1/VParams.nu));
|
||
|
|
gamma_eff = Params.gammaQF*2^(1/VParams.nu-1.5)*5^(-1/VParams.nu)*VParams.ell*gamma(1+1/VParams.nu)*( VParams.nu/(VParams.ell*gamma(1/VParams.nu)) )^(5/2);
|
||
|
|
EVar = VParams.nu^2*gamma(2-1/VParams.nu)/(8*VParams.ell^2*gamma(1/VParams.nu)) + 0.5*Params.gz*VParams.ell^2*gamma(3/VParams.nu)/gamma(1/VParams.nu);
|
||
|
|
|
||
|
|
% Parameters
|
||
|
|
normfac = Params.Lx*Params.Ly/numel(psi);
|
||
|
|
KEop = 0.5*(Transf.KX.^2+Transf.KY.^2);
|
||
|
|
|
||
|
|
% DDIs
|
||
|
|
frho = fftn(abs(psi).^2);
|
||
|
|
Phi = real(ifftn(frho.*VDk));
|
||
|
|
Eddi = (Params.gdd*Phi.*abs(psi).^2)/(sqrt(2*pi)*VParams.ell_eff);
|
||
|
|
|
||
|
|
%Kinetic energy
|
||
|
|
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
|
||
|
|
Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky/(2*pi)^2;
|
||
|
|
|
||
|
|
%Potential energy
|
||
|
|
Epot = V.*abs(psi).^2;
|
||
|
|
|
||
|
|
%Contact interactions
|
||
|
|
Eint = g_eff*abs(psi).^4;
|
||
|
|
|
||
|
|
%Quantum fluctuations
|
||
|
|
Eqf = gamma_eff*abs(psi).^5;
|
||
|
|
|
||
|
|
%Total energy
|
||
|
|
muchem = Ekin + EVar*Params.N + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy; %
|
||
|
|
muchem = muchem / Params.N; %Only use if psi is normalized to N
|
||
|
|
end
|