function [ result ] = atorbit(xlocal,ylocal,zlocal,Zlocal,laplace,type)
%ATORBIT Atomic orbital centered at specified coordinates
%
% xlocal  ... x-coordinate of the center of the atom (in m)
% ylocal  ... y-coordinate (in m)
% zlocal  ... z-coordinate (in m)
% Zlocal  ... atomic number of the atom
% laplace ... boolean. true - laplace operator is appled to atomic orbital
%                      false - "pure" atomic orbital is returned
% type    ... string specifying the orbital (1s, 2s, 2px, 2py, 2pz)
		
	result = [];
	
    %[x_mesh,xAtom_mesh] = meshgrid(x,xAtom);
    %xlocal = x - xAtom;
    %[y_mesh,yAtom_mesh] = meshgrid(y,yAtom);
    %ylocal = y - yAtom;
    %[z_mesh,zAtom_mesh] = meshgrid(z,zAtom);
    %zlocal = z - zAtom;
    
    %Zlocal = repmat(Z,length(x),1);

    a0 = 5.2917721E-11; % Bohr radius in meters
        
    r = sqrt(xlocal.^2 + ylocal.^2 + zlocal.^2);
    if ~laplace
        switch lower(type)
            case '1s'
                result = sqrt(Zlocal.^3./(pi*a0.^3)).*exp(-Zlocal.*r/a0);
            case '2s'
                result= 1/4*sqrt(Zlocal.^3./(2*pi*a0.^3)).*(2-Zlocal.*r/a0).*exp(-Zlocal.*r/(2*a0));
            case '2px'
                result= 1/4*sqrt(Zlocal.^5./(2*pi*a0.^5)).*xlocal.*exp(-Zlocal.*r/(2*a0));
            case '2py'
                result= 1/4*sqrt(Zlocal.^5./(2*pi*a0.^5)).*ylocal.*exp(-Zlocal.*r/(2*a0));

            case '2pz'
                result= 1/4*sqrt(Zlocal.^5./(2*pi*a0.^5)).*zlocal.*exp(-Zlocal.*r/(2*a0));

        end
    else
        switch(lower(type))

            case '1s' 
                 result = sqrt(Zlocal.^5./(pi*a0.^5)).*(Zlocal/a0 - 2./r).*exp(-Zlocal.*r/a0);
            case '2s'
                result = 1/4*sqrt(Zlocal.^5./(2*pi*a0.^5)).*(-Zlocal.^2.*r/(4*a0^2) + 5*Zlocal/(2*a0) - 4./r).*exp(-Zlocal.*r/(2*a0));
            case '2px'
                result = 1/4*sqrt(Zlocal.^7./(2*pi*a0.^7)).*(Zlocal.*r/(4*a0)- 2).*xlocal./r.*exp(-Zlocal.*r/(2*a0));
            case '2py'
                result = 1/4*sqrt(Zlocal.^7./(2*pi*a0.^7)).*(Zlocal.*r/(4*a0)- 2).*ylocal./r.*exp(-Zlocal.*r/(2*a0));
            case '2pz'
                result = 1/4*sqrt(Zlocal.^7./(2*pi*a0.^7)).*(Zlocal.*r/(4*a0)- 2).*zlocal./r.*exp(-Zlocal.*r/(2*a0));
        end
    end
end