function atomic_orbital(x,y,z,mu_components,meshsize,atomsize)
%  atomic_orbital(x,y,z, coefficients, meshsize, atomsize)
%
%       Plots angular component of atomic orbital 
%       or hybrid orbital on atom.
%
%       x,y,z = cartesian coordinates of atom 
%       coefficients = atomic orbital coefficients for s, p and d orbitals
%                    = [s, px, py, pz, dx2-y2, dz2, dxy, dxz, dyz];
%       meshsize     = fineness of spherical grid for plotting 
%       atomsize     = scale factor for radial plot.
%

currentfig=figure('Position',[-0,-0,650,650]);
orbitaltransparency=1;

s=mu_components(1);
px=mu_components(2);
py=mu_components(3);
pz=mu_components(4);
dx2=mu_components(5);
dz2=mu_components(6);
dxy=mu_components(7);
dxz=mu_components(8);
dyz=mu_components(9);


J=1.5*[s,px,py,pz,dx2,dz2,dxy,dxz,dyz];
atomsize=2.3*atomsize;
meshsize=round(3.8*meshsize);
orbsize=J*transpose(J); 
[Xj,Yj,Zj]=sphere(meshsize);
[phi,theta,sink]=cart2sph(Xj,Yj,Zj);

rs=(1/(4*pi))^0.5*s*ones(size(phi));
rpx=px*0.5*(3/pi)^0.5*cos(theta).*cos(phi);
rpy=py*0.5*(3/pi)^0.5*cos(theta).*sin(phi);
rpz=pz*0.5*(3/pi)^0.5*sin(theta);

rdz2=dz2*0.25*(5/(pi))^0.5*(3*sin(theta).*sin(theta)-ones(size(theta)));
rdxz=dxz*0.5*(15/(pi))^0.5*(sin(theta).*cos(theta).*cos(phi));
rdyz=dyz*0.5*(15/(pi))^0.5*(sin(theta).*cos(theta).*sin(phi));
rdxy=dxy*(15/(16*pi))^0.5*cos(theta).*cos(theta).*sin(2*phi);
rdx2=dx2*(15/(16*pi))^0.5*cos(theta).*cos(theta).*cos(2*phi);

thetaj=pi/2*ones(size(phi))-theta;
thetaj=abs(thetaj);
phij=phi;

rnew=(rs+(rpx+rpy+rpz)+(rdx2+rdz2+rdxy+rdxz+rdyz));

rnew=1.5*atomsize*rnew;  
max(max(max(rnew)));
[Xnew,Ynew,Znew]=sph2cart(phi,theta,abs(rnew));

C=rnew;
size(C);
d=.5*ones(size(C));
alpha=ones(size(C));
alphacm=ones(size(C));
m=size(C);
for n=1:m(1)
  for o=1:m(2)
  cmH(n,o,:)=[0.6,0.6,0.6];
  cmHe(n,o,:)=[.1,.1,.1];
  cm(n,o,:)=[0.0,140/255,0.0];
  cmS(n,o,:)=[0.7,0.7,.7];
  cmF(n,o,:)=[0.8,0.7,0.7];
  if C(n,o) < 0 
           d(n,o)=-1;
           alpha(n,o)=1;
           alphacm(n,o)=1.0;
           cm(n,o,:)=[0,90/255,1];
           cmH(n,o,:)=[0.0 0.0 0.0];
           cmHe(n,o,:)=[.7,.7,.7];
           cmF(n,o,:)=[0,100/255,0];
           cmS(n,o,:)=[0,0,1];

  end

  end
end
Xnew=Xnew+ones(size(Xnew))*x;
Ynew=Ynew+ones(size(Ynew))*y;
Znew=Znew+ones(size(Znew))*z;
hold on
ec=[0,0,0];
if rnew < .4 
        ec='none';
end
gotit=0;

hold on
g=surf(Xnew,Ynew,Znew,d,'AmbientStrength',.8,'FaceAlpha',1,'EdgeColor','none','Linewidth',0.5,'EdgeAlpha',1,'AlphaData',alphacm,'AlphaDataMapping','none','FaceColor','interp','CData',cm);
colormap([.2,.20,.20;.8 .8 .8]);
gotit=1;

%g=surf(Xnew,Ynew,Znew,d,'AmbientStrength',.8,'FaceAlpha','flat','EdgeColor','interp','AlphaData',0*alpha,'AlphaDataMapping','none','FaceColor','interp','CData',cmHe);
%colormap([.2,.20,.20;.8 .8 .8]);
%gotit=1;
