Graham Dennis
Calculate the 3D ground state of a Rubidium BEC in a harmonic magnetic trap assuming cylindrical symmetry about the z axis and reflection symmetry about z=0. This permits us to use the cylindrical bessel functions to expand the solution transverse to z and a cosine series to expand the solution along z.
const real omegaz = 2*M_PI*20; const real omegarho = 2*M_PI*200; const real hbar = 1.05457148e-34; const real M = 1.409539200000000e-25; const real g = 9.8; const real scatteringLength = 5.57e-9; const real Uint = 4.0*M_PI*hbar*hbar*scatteringLength/M; const real Nparticles = 5.0e5;

/* offset constants */ const real EnergyOffset = pow(15.0*Nparticles*Uint*omegaz*omegarho*omegarho/(8*M_PI), 2.0/5.0) * pow(M/2.0, 3.0/5.0);

t
V1
real Vtrap = 0.5*M*(omegarho*omegarho*r*r + omegaz*omegaz*z*z);

V1 = -i/hbar*(Vtrap - EnergyOffset);

eta
phi
noise

if ((abs(r) < 0.9e-5) && abs(z) < 0.9e-4) { phi = eta; // This will be automatically normalised later } else { phi = 0.0; }

Ncalc
wavefunction
// Calculate the current normalisation of the wave function. Ncalc = mod2(phi);
300
wavefunction normalisation
// Correct normalisation of the wavefunction phi *= sqrt(Nparticles/Ncalc);
T
T = -0.5*hbar/M*(kr*kr + kz*kz);
wavefunction
potential
dphi_dt = T[phi] - (i*V1 + Uint/hbar*mod2(phi))*phi;
norm_dens
wavefunction normalisation
norm_dens = mod2(phi)/Ncalc;