INSTRUCTIONS:

Save everthing below the first line starting with (* into a file,
say, orbital.m. In Mathematica, load this file as input:

	<< orbital.m <SHIFT-Enter>

Then, you can type in commands like

	orbit2p[] <SHIFT-Enter>

to get the default 2p surface, and so on. Please read the
usage notes included as comments in the file below for more
details.

(* ================================================================ *)
(*       Mathematica program to generate atomic orbitals            *)
(*                                                                  *)	
(*                       written by                                 *)
(*                                                                  *)
(*                 Dr. Pui-Kan Catherine Kong                       *)
(*                 Department of Mathematics                        *)
(*                 Carson-Newman College                            *)
(*                 Jefferson City, TN 37760                         *)
(*                                                                  *)
(*	              kong@cncacc.cn.edu			    *)
(* ================================================================ *)

BeginPackage["orbital`"]

Unprotect[larger, smaller, cp, cn, orbit1s, orbit2s, orbit3s, orbitns,
	  orbit2p, orbit3p, orbitnp, orbit3dz2, orbit3dx2y2, orbit3d,
	  orbital]

Clear[larger, smaller, cp, cn, orbit1s, orbit2s, orbit3s, orbitns,
	orbit2p, orbit3p, orbitnp, orbit3dz2, orbit3dx2y2, orbit3d,
	orbital]

larger::usage = "larger[x] is a function used to round off a number
x to 5 decimal places."

smaller::usage = "smaller[x] is a function used to truncate a number
x to 5 decimal places."

cp::usage = "cp = GrayLevel[0.5] defines the positive surface"

cn::usage = "cn = GrayLevel[0.9] defines the negative surface"

(* ================================================================ *)

orbit1s::usage = "orbit1s[scale] produces an isosurface for the
1s orbital. \"scale\" is the ratio between 0.1 and 0.9 of the
smallest absolute maximum value of the radial function. The
default is orbit1s[0.5]."

orbit2s::usage = "orbit2s[scale] produces an isosurface for the
2s orbital. \"scale\" is a number between 0.1 and 0.9. The default
is orbit2s[0.5]."

orbit3s::usage = "orbit3s[scale] produces an isosurface for the
3s orbital. \"scale\" is a number between 0.1 and 0.9. The default
is orbit3s[.5]."

orbitns::usage = "orbitns[n,scale] produces an isosurface for the
orbitals of 1s, 2s or 3s according to the value of \"n\". The
default is orbitns[1,0.5] which is 1s."

(* ================================================================ *)

orbit2p::usage = "orbit2p[scale, type] produces an isosurface
for the orbitals of 2px, 2py or 2pz. The default is orbit2p[.5,0]
which is 2pz. \n \n
	\t - \"scale\" is a number between 0.1 and 0.9 \n \n
	\t - \"type\" = 1 gives 2px \n
	\t \t  \t   -1 gives 2py \n
	\t \t  \t    0 or otherwise gives 2pz \n"

orbit3p::usage = "orbit3p[scale, type] produces an isosurface
for the orbitals of 3px, 3py or 3pz. The default is orbit3p[0.6,0]
which is 3pz. \n \n
	\t - \"scale\" is a number between 0.1 and 0.9 \n \n
	\t - \"type\" = 1 gives 3px \n
	\t \t  \t   -1 gives 3py \n
	\t \t  \t    0 or otherwise gives 3pz \n"

orbitnp::usage = "orbitnp[n,scale,type] produces an isosurface
for the orbitals of 2p or 3p according to the value of \"n\".
The default is orbitnp[2,0.5,0] which is 2pz. \n \n
\t - \"scale\" is a number between 0.1 and 0.9 \n \n
	\t - \"n\" = 2 : \"type\" = 1 gives 2px \n
	\t             \"type\" =-1 gives 2py \n
	\t             \"type\" = 0 gives 2pz \n \n
	\t - \"n\" = 3 : \"type\" = 1 gives 3px \n
	\t             \"type\" =-1 gives 3py \n
	\t             \"type\" = 0 gives 3pz \n "

(* ================================================================ *)

orbit3dz2::usage = "orbit3dz2[scale] produces an isosurface for the
3dz2 orbital. \"scale\" is a number between 0.1 and 0.9. The default
is orbit3dz2[0.3]."

orbit3dx2y2::usage = "orbit3dx2y2[scale] produces an isosurface for 
the 3dx2y2 orbital. \"scale\" is a number between 0.1 and 0.9. The
default is orbit3dx2y2[0.5]."

orbit3d::usage = "orbit3d[scale,type] produces an isosurface for the
3dz2 or 3dx2y2. The default is orbit3d[0.3,0] which is 3dz2.\n \n
	\t - \"scale\" is a number between 0.1 and 0.9 \n \n
	\t - \"type\" = 1 gives 3dx2y2 \n
	\t \t \t    0 or otherwise gives 3dz2 \n"

(* ================================================================ *)

orbital::usage = "orbital[n,l,scale,type] produces an isosurface for
the orbital according to the values of \"n\", \"l\" and \"type\". The
default is orbital[2,1,0.3,0] which is 2pz. \n \n
	\t - \"scale\" is a number between 0.1 and 0.9 \n \n
	\t - \"l\" = 0 : \"n\" = 1 gives 1s \n
	\t             \"n\" = 2 gives 2s \n
	\t             \"n\" = 3 gives 3s \n \n
	\t - \"l\" = 1 : \"n\" = 2 : \"type\" = 1 gives 2px \n
	\t                       \"type\" =-1 gives 2py \n
	\t                       \"type\" = 0 gives 2pz \n \n
	\t           : \"n\" = 3 : \"type\" = 1 gives 3px \n
	\t                       \"type\" =-1 gives 3py \n
	\t                       \"type\" = 0 gives 3pz \n \n
	\t - \"l\" = 2 : \"type\" = 0 gives 3dz2 \n
	\t             \"type\" = 1 gives 3dx2y2 \n"

(* ---------------------------------------------------------------- *)

Begin["`Private`"]

larger[x_]:=N[(Floor[x*100000]+1)/100000]

smaller[x_]:=N[Floor[x*100000]/100000]

(*cp=RGBColor[.6,0,.5]*)
(*cn=RGBColor[0,.6,.5]*)
cp=GrayLevel[0.5];  (* positive *)
cn=GrayLevel[0.9];  (* negative *)

(* ---------------------------------------------------------------- *)

orbit1s[scale_:.5]:=
	Module[{a,g,gmax,fx,fy,fz,k,lab,
		root},
		
		g[x_]:=Exp[-x];
		gmax=g[.5];
		k=If[scale<=.1,.1*gmax,
		     If[scale>=.9,.9*gmax,
			scale*gmax]];
	
		root=FindRoot[g[x]==k,{x,2}];
		a=larger[root[[1,2]]];

		fx[x_,y_]:=a*Sin[x]*Cos[y];
		fy[x_,y_]:=a*Sin[x]*Sin[y];
		fz[x_]:=a*Cos[x];

		lab=StringForm["1s"];

ParametricPlot3D[{fx[x,y],fy[x,y],fz[x],cp},
		 {x,0,Pi},{y,0,2Pi},PlotLabel->lab,Lighting->False]

	       ] (* End of Module for orbit1s *)


(* ---------------------------------------------------------------- *)

orbit2s[scale_:.5]:=
	Module[{a,b,c,g,gmax,fx,fy,fz,k,lab,
		root,a1,r,dx,dy},
		
		g[x_]:=(2-x)*Exp[-x/2];
		gmax=-1*g[4];
		k=If[scale<=.1,.1*gmax,
		     If[scale>=.9,.9*gmax,
			scale*gmax]];
	
		root=FindRoot[g[x]==k,{x,1}];
		a=larger[root[[1,2]]];
		root=FindRoot[g[x]==-k,{x,3}];
		b=larger[root[[1,2]]];
		root=FindRoot[g[x]==-k,{x,5}];
		c=smaller[root[[1,2]]];

		a1=Pi/(c-b);
		r[x_]:=b+x/a1;

		fx[x_,y_]:=Sin[x]*Cos[y];
		fy[x_,y_]:=Sin[x]*Sin[y];
		fz[x_]:=Cos[x];

		lab=StringForm["2s"];
		dx=Pi/10;
		dy=Pi/10;

ParametricPlot3D[{{a*fx[x,2y],a*fy[x,2y],a*fz[x],cp},
		  {b*fx[x,y],b*fy[x,y],b*fz[x],cn},
		  {c*fx[x,y],c*fy[x,y],c*fz[x],cn},
		  {r[x]*Cos[y-Pi/2],0,r[x]*Sin[y-Pi/2],cn},
		  {r[x]*Cos[-y-Pi/2],0,r[x]*Sin[-y-Pi/2],cn}},
		  {x,0,Pi},{y,0,Pi},PlotLabel->lab,Lighting->False,
		  PlotRange->{{-c-1,c+1},{-a-1,c+1},{-c-1,c+1}}]

	       ] (* End of Module for orbit2s *)


(* ---------------------------------------------------------------- *)

orbit3s[scale_:.5]:=
	Module[{a,b,c,d,e,g,gmax,fx,fy,fz,k,lab,
		root,a1,a2,r1,r2,dx,dy},
		
		g[x_]:=(27-18x+2x^2)*Exp[-x/3];
		gmax=g[15/2+3*Sqrt[7]/2];
		k=If[scale<=.1,.1*gmax,
		     If[scale>=.9,.9*gmax,
			scale*gmax]];
	
		root=FindRoot[g[x]==k,{x,1}];
		a=larger[root[[1,2]]];
		root=FindRoot[g[x]==-k,{x,3}];
		b=larger[root[[1,2]]];
		root=FindRoot[g[x]==-k,{x,5}];
		c=smaller[root[[1,2]]];
		root=FindRoot[g[x]==k,{x,9}];
		d=larger[root[[1,2]]];
		root=FindRoot[g[x]==k,{x,15}];
		e=smaller[root[[1,2]]];

		a1=Pi/(c-b);
		a2=Pi/(e-d);
		r1[x_]:=b+x/a1;
		r2[x_]:=d+x/a2;

		fx[x_,y_]:=Sin[x]*Cos[y];
		fy[x_,y_]:=Sin[x]*Sin[y];
		fz[x_]:=Cos[x];

		lab=StringForm["3s"];
		dx=Pi/10;
		dy=Pi/10;

ParametricPlot3D[{{a*fx[x,2y],a*fy[x,2y],a*fz[x],cp},
		  {b*fx[x,y],b*fy[x,y],b*fz[x],cn},
		  {c*fx[x,y],c*fy[x,y],c*fz[x],cn},
		  {d*fx[x,y],d*fy[x,y],d*fz[x],cp},
		  {e*fx[x,y],e*fy[x,y],e*fz[x],cp},
		  {r1[x]*Cos[y-Pi/2],0,r1[x]*Sin[y-Pi/2],cn},
		  {r1[x]*Cos[-y-Pi/2],0,r1[x]*Sin[-y-Pi/2],cn},
		  {r2[x]*Cos[y-Pi/2],0,r2[x]*Sin[y-Pi/2],cp},
		  {r2[x]*Cos[-y-Pi/2],0,r2[x]*Sin[-y-Pi/2],cp}},
		  {x,0,Pi},{y,0,Pi},PlotLabel->lab,Lighting->False,
		  PlotRange->{{-e-1,e+1},{-a-1,e+1},{-e-1,e+1}}]

	       ] (* End of Module for orbit3s *)

(* ---------------------------------------------------------------- *)

orbitns[n_:1,scale_:.5]:= If[n==2, orbit2s[scale],
		If[n==3, orbit3s[scale], orbit1s[scale]]]

		(* End of orbitns *)

(* ---------------------------------------------------------------- *)
(* ================================================================ *)

orbit2p[scale_:.5,type_:0]:=
	Module[{a,b,g,gmax,f,fx,fy,fz,k,labx,laby,
		labz,root},
		
		g[x_]:=x*Exp[-x/2];
		gmax=g[2];		
		k=If[scale<=.1,.1*gmax,
		     If[scale>=.9,.9*gmax,
			scale*gmax]];
	
		root=FindRoot[g[x]==k,{x,0}];
		a=larger[root[[1,2]]];
		root=FindRoot[g[x]==k,{x,4}];
		b=smaller[root[[1,2]]];
		k=If[N[g[a]]>=N[g[b]],g[b],g[a]];

		f[x_]:=Sqrt[1-k^2/g[x]^2];

		fx[x_,y_]:=x*f[x]*Cos[y];
		fy[x_,y_]:=x*f[x]*Sin[y];
		fz[x_]:=k*Exp[x/2];

		labx=StringForm["2px"];
		laby=StringForm["2py"];
		labz=StringForm["2pz"];

		If[type==1,
		
ParametricPlot3D[{{fz[x],fx[x,y],fy[x,y],cp},
		 {-fz[x],fx[x,y],fy[x,y],cn}},
		 {x,a,b},{y,0,2Pi},PlotLabel->labx,Lighting->False,
		 PlotRange->{{-b-.1,b+.1},{-b+a/b,b-a/b},{-b+a/b,b-a/b}}],

If[type==-1,

ParametricPlot3D[{{fx[x,y],fz[x],fy[x,y],cp},
		 {fx[x,y],-fz[x],fy[x,y],cn}},
		 {x,a,b},{y,0,2Pi},PlotLabel->laby,Lighting->False,
		 PlotRange->{{-b+a/b,b-a/b},{-b-.1,b+.1},{-b+a/b,b-a/b}}],

ParametricPlot3D[{{fx[x,y],fy[x,y],fz[x],cp},
		 {fx[x,y],fy[x,y],-fz[x],cn}},
		 {x,a,b},{y,0,2Pi},PlotLabel->labz,Lighting->False,
		 PlotRange->{{-b+a/b,b-a/b},{-b+a/b,b-a/b},{-b-.1,b+.1}}]

] (* End of If *)

		  ] (* End of If *)

	       ] (* End of Module for orbit2p *)

(* ---------------------------------------------------------------- *)

orbit3p[scale_:.6,type_:0]:=
	Module[{a,b,c,d,g,gmax,root,p,n,
		hp,hn,a1,a2,fx1,fy1,fz1,
		fx2,fy2,fz2,labx,laby,labz},

		g[x_]:=x*(6-x)*Exp[-x/3];

		gmax=-1*g[6+3*Sqrt[2]];
		p=If[scale<=.1,.1*gmax,
		     If[scale>=.9,.9*gmax,
			scale*gmax]];

		n=-p;

		root=FindRoot[g[x]==p,{x,0}];
		a=larger[root[[1,2]]];
		root=FindRoot[g[x]==p,{x,5.75}];
		b=smaller[root[[1,2]]];
		p=If[N[g[a]]>=N[g[b]],g[b],g[a]];

		root=FindRoot[g[x]==n,{x,6}];
		c=larger[root[[1,2]]];
		root=FindRoot[g[x]==n,{x,12}];
		d=smaller[root[[1,2]]];

		a1=(b-a)/(d-c);
		d=(b-a)/a1+c;
		n=If[N[g[c]]<=N[g[d]],g[d],g[c]];

		hp[x_]:=Sqrt[1-p^2/(g[x])^2];
		hn[x_]:=Sqrt[1-n^2/(g[x])^2];


		fx1[x_,y_]:=x*hp[x]*Cos[y];
		fy1[x_,y_]:=x*hp[x]*Sin[y];
		fz1[x_]:=x*p/g[x];

		fx2[x_,y_]:=((x-a)/a1+c)*hn[(x-a)/a1+c]*Cos[y];
		fy2[x_,y_]:=((x-a)/a1+c)*hn[(x-a)/a1+c]*Sin[y];
		fz2[x_]:=((x-a)/a1+c)*n/g[(x-a)/a1+c];

		labx=StringForm["3px"];
		laby=StringForm["3py"];
		labz=StringForm["3pz"];

		If[type==1,

ParametricPlot3D[{{fz1[x],fx1[x,y],fy1[x,y],cp},
		  {fz2[x],fx2[x,y],fy2[x,y],cn},
		  {-fz1[x],fx1[x,y],fy1[x,y],cn},
		  {-fz2[x],fx2[x,y],fy2[x,y],cp}},
		  {x,a,b},{y,0,2Pi},PlotLabel-> labx,Lighting->False,
		  PlotRange->{{-d-.1,d+.1},{-d+c/d,d-c/d},{-d+c/d,d-c/d}}],

If[type==-1,

ParametricPlot3D[{{fx1[x,y],fz1[x],fy1[x,y],cp},
		  {fx2[x,y],fz2[x],fy2[x,y],cn},
		  {fx1[x,y],-fz1[x],fy1[x,y],cn},
		  {fx2[x,y],-fz2[x],fy2[x,y],cp}},
		  {x,a,b},{y,0,2Pi},PlotLabel-> laby,Lighting->False,
		  PlotRange->{{-d+c/d,d-c/d},{-d-.1,d+.1},{-d+c/d,d-c/d}}],

ParametricPlot3D[{{fx1[x,y],fy1[x,y],fz1[x],cp},
		  {fx2[x,y],fy2[x,y],fz2[x],cn},
		  {fx1[x,y],fy1[x,y],-fz1[x],cn},
		  {fx2[x,y],fy2[x,y],-fz2[x],cp}},
		  {x,a,b},{y,0,2Pi},PlotLabel -> labz,Lighting->False,
		  PlotRange->{{-d+c/d,d-c/d},{-d+c/d,d-c/d},{-d-.1,d+.1}}]

] (* End of If *)

		 ] (* End of If *)

	       ] (* End of Module for orbit3p *)

(* ---------------------------------------------------------------- *)

orbitnp[n_:2,scale_:.5,type_:0]:=
	If[n==3,orbit3p[scale,type],orbit2p[scale,type]]
	
	(* End of orbitnp *)

(* ---------------------------------------------------------------- *)
(* ================================================================ *)

orbit3dz2[scale_:.3]:=
	Module[{g,p,n,root,a,b,c,d,
		fx1,fy1,fz1,fx2,fy2,fz2,gmax,
		a1,a2,lab},

       		g[x_]:=x^2*Exp[-x/3];
		gmax=g[6];
		p=If[scale<=.1,.1*gmax,
		     If[scale>=.9,.9*gmax,
			scale*gmax]];

		n=p;

		root=FindRoot[2*g[x]==p,{x,6}];
		c=larger[root[[1,2]]];
		root=FindRoot[2*g[x]==p,{x,10}];
		d=smaller[root[[1,2]]];
		p=2*If[N[g[c]]>=N[g[d]],g[d],g[c]];

		fx1[x_,y_]:=x*Sqrt[1/3(2-p/g[x])]*Cos[y];
		fy1[x_,y_]:=x*Sqrt[1/3(2-p/g[x])]*Sin[y];
		fz1[x_]:=x*Sqrt[1/3(1+p/g[x])];

		root=FindRoot[g[x]==n,{x,6}];
		a=larger[root[[1,2]]];
		root=FindRoot[g[x]==n,{x,10}];
		b=smaller[root[[1,2]]];
		n=If[N[g[a]]>=N[g[b]],g[b],g[a]];
		a1=(d-c)/(b-a);
		a2=c - a*a1;
		
		fx2[x_,y_]:=(x-a2)/a1*Sqrt[1/3(2+n/g[(x-a2)/a1])]*Cos[y];
		fy2[x_,y_]:=(x-a2)/a1*Sqrt[1/3(2+n/g[(x-a2)/a1])]*Sin[y];
		fz2[x_]:=(x-a2)/a1*Sqrt[1/3(1-n/g[(x-a2)/a1])];

		lab=StringForm["3dz2"];


ParametricPlot3D[{{fx1[x,y],fy1[x,y],fz1[x],cp},
		  {fx1[x,y],fy1[x,y],-fz1[x],cp},
		  {fx2[x,y],fy2[x,y],fz2[x],cn},
		  {fx2[x,y],fy2[x,y],-fz2[x],cn}},
		  {x,c,d},{y,0,2Pi},PlotLabel -> lab,Lighting->False]

]  (* End of orbit3dz2 *)

(* ---------------------------------------------------------------- *)

orbit3dx2y2[scale_:.5]:=
	Module[{a,b,g,gmax,k,root,h,u,
		fx,fy,fz,labxy},

		g[x_]:=x^2*Exp[-x/3];

		gmax=.9*g[6];
		k=If[scale<=.1,.1*gmax,
		     If[scale>=.9,.9*gmax,
			scale*gmax]];

		root=FindRoot[g[x]==k,{x,6}];
		a=larger[root[[1,2]]];
		root=FindRoot[g[x]==k,{x,10}];
		b=smaller[root[[1,2]]];
		k=If[N[g[a]]>=N[g[b]],g[b],g[a]];

		u[x_,y_]:=y/2*ArcCos[k/g[x]];
		h[x_,y_]:=k/(g[x]*Cos[2*u[x,y]]);

		fx[x_,y_]:=x*Cos[u[x,y]]*Sqrt[h[x,y]];
		fy[x_,y_]:=x*Sin[u[x,y]]*Sqrt[h[x,y]];
		fz[x_,y_]:=x*Sqrt[1-h[x,y]];

		labxy=StringForm["3dx2-y2"];

ParametricPlot3D[{{fx[x,y],fy[x,y],fz[x,y],cp},
		  {-fx[x,y],fy[x,y],fz[x,y],cp},
		  {fx[x,y],fy[x,y],-fz[x,y],cp},
		  {-fx[x,y],fy[x,y],-fz[x,y],cp},
		  {fy[x,y],-fx[x,y],fz[x,y],cn},
		  {fy[x,y],fx[x,y],fz[x,y],cn},
		  {fy[x,y],-fx[x,y],-fz[x,y],cn},
		  {fy[x,y],fx[x,y],-fz[x,y],cn}},
		  {x,a,b},{y,-.9999999,.9999999},
		  PlotLabel->labxy,Lighting->False]

]  (* End of orbit3dx2y2*)

(* ---------------------------------------------------------------- *)

orbit3d[scale_:.3,type_:0]:=
	If[type==1, orbit3dx2y2[scale], orbit3dz2[scale]]

	    (* End of orbit3d *)

(* ---------------------------------------------------------------- *)
(* ================================================================ *)

orbital[n_:2,l_:1,scale_:.3,type_:0]:=
	If[l==2, orbit3d[scale,type],
		If[l==0, orbitns[n,scale,type],
		 	 orbitnp[n,scale,type]]]

	(* End of orbital *)

(* ---------------------------------------------------------------- *)

End[] 	(* orbital`Private` *)

Protect[larger, smaller, cn, cp, orbit1s, orbit2s, orbit3s, orbitns,
	orbit2p, orbit3p, orbitnp, orbit3dz2, orbit3dx2y2, orbit3d,
	orbital]

EndPackage[]	(* orbital` *) 



(* ---------------------------------------------------------------- *)