function [psi] = setupWavefunction(~,Params,Transf)

    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; 
    
    Rx   = 4*ellx; 
    Ry   = 4*elly; 
    Rz   = 4*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 = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
    psi      = psi/sqrt(cur_norm);
    
    % --- Adding some noise --- 
    r     = normrnd(0,1,size(X));
    theta = rand(size(X));
    noise = r.*exp(2*pi*1i*theta);
    psi   = psi + Params.nsf*noise;
    
    % Renormalize wavefunction
    Norm  = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
    psi   = sqrt(Params.N)*psi/sqrt(Norm);

end