(*

The wave function (properly normalized)

		Psi[n,l,m]

with the dependence on the azimuthal angle phi suppressed.

Radial part of the wave function (not normalized)

		RadPsi[n,l]

Angular part of the wave fucntion (not normalized)

		Y[l,m]

Plotting:

	1. Angular wave function	PlotYY[l,m]
	2. Radial wave function		PlotRR[n,l]
	3. Probability density		PlotProb[n,l,m,res]

where 

	n = 1, 2, ... 	 	principal quantum number
	l = 0, 1, ... , n-1	angular quantum number
	m = -l, -l+1, ... , l	magnetic quantum number

	res = 1, 2, ... is the resolution, slow if res>4
 
*)

<<Graphics`Graphics`

Clear[RadPsi]
RadPsi[n_,l_]:=Rx[n,l,2 r/n]
  

Clear[Y,YY]
YY[l_,m_,th_]:=SphericalHarmonicY[l,m,th,0]
Y[l_,m_]:=SphericalHarmonicY[l,m,"\[Theta]","\[Phi]"]


Clear[LL,Rx,Rxn,RR,RRn]
LL[n_,a_,x_]:=LaguerreL[n,a,x]
Rx[n_,l_,x_]:=Exp[-x/2] (x)^l LaguerreL[n-l-1,2 l+1,x]
Rxn[n_,l_,x_]:=Exp[-x/2] (x+0.0000000001)^l LaguerreL[n-l-1,2 l+1,x]
RR[n_,l_,r_]:=(*((2/n)^3 (n-l-1)!/(2 n ((n+l)!)^3))^(1/2) *)Rx[n,l,2 r/n]
RRn[n_,l_,r_]:=N[Rxn[n,l,2 r/n]]


Clear[Psi,Prob]
Psi[n_,l_,m_]:=((2/n)^3 (n-l-1)!/(2 n ((n+l)!)^3))^(1/2) Rx[n,l,2 r/n]*
		 SphericalHarmonicY[l,m,"\[Theta]",0]

Clear[tth]
tth[x_,y_]:=ArcTan[y/x]
tth[0,y_]:=Pi/2
tth[0,0]:=0
tth[0.,y_]:=Pi/2
tth[0.,0.]:=0
Prob[n_,l_,m_,x_,y_]:=Block[{rx,tx,ex},
	rx=Sqrt[x^2+y^2];
	tx=tth[x,y];
	(RRn[n,l,rx] N[YY[l,m,tx]])^2]
Prob[n_,l_,m_,x_,y_]:=Prob[n,l,m,-x,y]/;x<0
Prob[n_,l_,m_,x_,y_]:=Prob[n,l,m,x,-y]/;y<0


Clear[PlotYY]
PlotYY[l_,m_]:=Block[{yx,px},
	yx=YY[l,m,x];
	px=PolarPlot[yx^2,{x,0,2 Pi},Ticks->None, ImageSize->72 4]]

PlotRR[n_,l_]:=Plot[RR[n,l,r],{r,0,5/2 n^2},Ticks->{None,None},
		PlotRange->All]

PlotYYsm[l_,m_]:=Block[{yx,px},
	yx=YY[l,m,x];
	px=PolarPlot[yx^2,{x,0,2 Pi},Ticks->None]]

PlotRRsm[n_,l_]:=Plot[RR[n,l,r],{r,0,5/2 n^2},Ticks->None]


(* Color scale *)

Clear[cff]
cff[x_]:=Hue[1/2-1/2 x]

ProbData[n_,l_,m_,res_]:=Block[{ex1,ex2},
	ex1=5./2 n^2;
	Table[N[Prob[n,l,m,Chop[-ex1+(ex1/10/res) i],
		Chop[-ex1+(ex1/10/res) j]]],
		{i,0,20 res+1},{j,0,20 res+1}]]

PlotProb[n_,l_,m_,res_]:=DensityPlot[
	Prob[n,l,m,x,y],
		{x,-5/2 n^2,5/2 n^2},{y,-5/2 n^2,5/2 n^2},
		ColorFunction->cff,PlotPoints->6 res,
		Mesh->False,
		Frame->False]
	
(*

Clear[PlotProb]
PlotProb[n_,l_,m_,res_]:=Block[{rex,tex,ex1,ex2,xre,tre,ex3},
	rex=Table[N[5/Sqrt[2] n^2/10/res ii],{ii,0,10 res}];
	ex1=Map[N[RR[n,l,#]]^2&,rex,{1}];
	tex=Table[2 Pi/5/res jj,{jj,0,5 res-1}];
	ex2=Map[N[YY[l,m,#]]^2&,tex,{1}];
	sex=Map[N[Sin[#]]&,tex,{1}];
	cex=Map[N[Cos[#]]&,tex,{1}];
	xre=Table[{rex[[i]],ex1[[i]]},{i,Length[rex]}];
	tre=Table[{cex[[j]],sex[[j]],ex2[[j]]},{j,Length[cex]}];
	ex3=Flatten[
		Table[{xre[[i]],tre[[j]]},{i,Length[xre]},{j,Length[tre]}],
			1]]

Clear[toxyp]
toxyp[{{r_,rr_},{c_,s_,y_}}]:={r c, r s, rr y}
		
test:=DensityPlot[Exp[Sqrt[x^2+y^2]] Sin[y]^2,{x,0,2  Pi},{y,0,Pi},
PlotPoints->200,
	ColorFunction->cff,Mesh->False, ImageSize->72 8]

test2:=DensityPlot[(2^(3/2)/(8 Pi)^(1/2) Exp[-Sqrt[x^2+y^2]])^2,
       {x,-5,5},{y,-5,5},
PlotPoints->20,
	ColorFunction->cff,Mesh->False]


Show[GraphicsArray[{%21,%10,%8}],ImageSize->72 {8,2.5},
         PlotLabel->"n=6 l=3 m=0"] *)