plotOrbital[l_, m_] := 
  Module[{psi, plotPoints, plotRange,plot1,plot2,plot3},
         plotPoints = 40;
  plotRange = .8;
  psi = SphericalHarmonicY[l,m,Θ,ϕ];
  Print[ "l = ",l,"   m = ",m, "   psi = ",Re[psi]];
  If[m< 0, psi = Exp[ ((-1)^(-m))I Pi /2]SphericalHarmonicY[l,m,Θ,ϕ]];  (* rotation correction to match wikipedia table *)
  plot1 = SphericalPlot3D[Max[Re[psi],0],{Θ,0,π},{ϕ,0,2 π},ColorFunction->Function[{x,y,z,Θ,ϕ,r},RGBColor[1,0,0]], Mesh->None, PlotPoints->plotPoints, Boxed->False, Axes->None, PlotRange->All];
  plot2 = SphericalPlot3D[-Min[Re[psi],0],{Θ,0,π},{ϕ,0,2 π},ColorFunction->Function[{x,y,z,Θ,ϕ,r},RGBColor[0,1,1]], Mesh->None, PlotPoints->plotPoints, Boxed->False, Axes->None, PlotRange->All];
  plot3 = Show[plot1,plot2,ViewPoint->{ -2.4,0,1.5}];
  plot3
  ];
plotOrbitalUp[l_, m_] := 
  Module[{psi, plotPoints, plotRange,plot1,plot2,plot3,color1,color2,color},
         plotPoints = 40;
  plotRange = .8;
  color1 = RGBColor[.8,.5,.8];
  color2 = RGBColor[.6,.5,1];
  If[m+1/2 > 0, color = color1, color = color2];
  psi = SphericalHarmonicY[l,m,Θ,ϕ]*Exp[I (m+1/2)ϕ]/Exp[I (m)ϕ] ;
  Print[ "l = ",l,"   m = ",m, "   psiUp = ",Re[psi]];
  If[m< 0, psi = Exp[ ((-1)^(-m))I Pi /2]psi];  (* rotation correction to match wikipedia table *)
  plot1 = SphericalPlot3D[Max[Re[psi],0],{Θ,0,π},{ϕ,0,2 π},ColorFunction->Function[{x,y,z,Θ,ϕ,r},color], Mesh->None, PlotPoints->plotPoints, Boxed->False, Axes->None, PlotRange->All];
  plot2 = SphericalPlot3D[-Min[Re[psi],0],{Θ,0,π},{ϕ,0,2 π},ColorFunction->Function[{x,y,z,Θ,ϕ,r},color], Mesh->None, PlotPoints->plotPoints, Boxed->False, Axes->None, PlotRange->All];
  plot3 = Show[plot1,plot2,ViewPoint->{ -2.4,0,1.5}];
  plot3
  ];
plotOrbitalDown[l_, m_] := 
  Module[{psi, plotPoints, plotRange,plot1,plot2,plot3,color1, color2, color},
         plotPoints = 40;
  plotRange = .8;
  color1 = RGBColor[.8,.5,.8];
  color2 = RGBColor[.6,.5,1];
  If[m-1/2 > 0, color = color1, color = color2];
  psi = SphericalHarmonicY[l,m,Θ,ϕ]*Exp[I (m-1/2)ϕ]/Exp[I (m)ϕ] ;
  Print[ "l = ",l,"   m = ",m, "   psiDown = ",Re[psi]];
  If[m< 0, psi = Exp[ ((-1)^(-m))I Pi /2]psi];  (* rotation correction to match wikipedia table *)
  plot1 = SphericalPlot3D[Max[Re[psi],0],{Θ,0,π},{ϕ,0,2 π},ColorFunction->Function[{x,y,z,Θ,ϕ,r},color], Mesh->None, PlotPoints->plotPoints, Boxed->False, Axes->None, PlotRange->All];
  plot2 = SphericalPlot3D[-Min[Re[psi],0],{Θ,0,π},{ϕ,0,2 π},ColorFunction->Function[{x,y,z,Θ,ϕ,r},color], Mesh->None, PlotPoints->plotPoints, Boxed->False, Axes->None, PlotRange->All];
  plot3 = Show[plot1,plot2,ViewPoint->{ -2.4,0,1.5}];
  plot3
  ];

plotOrbitalUp[0,0]

plotOrbitalDown[0,0]

plotOrbitalDown[4,4]

plotOrbitalDown[4,-4]

plotOrbitalUp[4,4]

plotOrbitalUp[4,-4]

plotOrbital[4,4]

plotOrbital[4,-4]

plotOrbital[4,4]
plotOrbitalUp[4,4]

plotOrbital[4,4]
plotOrbital[4,3]

plotOrbital[4,-4]
plotOrbital[4,4]

plotOrbital[3,-3]
plotOrbital[3,3]

plotOrbital[2,-2]
plotOrbital[2,2]

plotOrbital[1,-1]
plotOrbital[1,1]

lmax =3;
majorArray = {};
Do[
(* minorArray = {};
Do[
AppendTo[minorArray,plotOrbital[l,m]];
,{m,-lmax,lmax}];
AppendTo[majorArray,GraphicsRow[minorArray,Spacings->-40,ImageSize->{100*(lmax*2+1),70}]]; *)
minorArray = {}; 
Do[
AppendTo[minorArray,plotOrbitalUp[l,m]];
,{m,-lmax,lmax}];
AppendTo[majorArray,GraphicsRow[minorArray,Spacings->-40,ImageSize->{100*(lmax*2+1),70}]];
minorArray = {};
Do[
AppendTo[minorArray,plotOrbitalDown[l,m]];
,{m,-lmax,lmax}];
AppendTo[majorArray,GraphicsRow[minorArray,Spacings->-40,ImageSize->{100*(lmax*2+1),70}]];
,{l,0,lmax}];

GraphicsColumn[majorArray]

lmax = 2;
majorArray = {};
Do[
(* minorArray = {};
Do[
AppendTo[minorArray,plotOrbital[l,m]];
,{m,-lmax,lmax}];
AppendTo[majorArray,GraphicsRow[minorArray,Spacings->-40,ImageSize->{100*(lmax*2+1),70}]]; *)
minorArray = {}; 
Do[
AppendTo[minorArray,plotOrbitalUp[l,m]];
,{m,-lmax,lmax}];
AppendTo[majorArray,GraphicsRow[minorArray,Spacings->-40,ImageSize->{100*(lmax*2+1),70}]];
minorArray = {};
Do[
AppendTo[minorArray,plotOrbitalDown[l,m]];
,{m,-lmax,lmax}];
AppendTo[majorArray,GraphicsRow[minorArray,Spacings->-40,ImageSize->{100*(lmax*2+1),70}]];
,{l,0,lmax}];

GraphicsColumn[majorArray]

lmax = 2;
majorArray = {};
Do[
minorArray = {};
Do[
AppendTo[minorArray,plotOrbital[l,m]];
,{m,-lmax,lmax}];
AppendTo[majorArray,GraphicsRow[minorArray,Spacings->-40,ImageSize->{100*(lmax*2+1),70}]];
minorArray = {};
Do[
AppendTo[minorArray,plotOrbitalUp[l,m]];
,{m,-lmax,lmax}];
AppendTo[majorArray,GraphicsRow[minorArray,Spacings->-40,ImageSize->{100*(lmax*2+1),70}]];
minorArray = {};
Do[
AppendTo[minorArray,plotOrbitalDown[l,m]];
,{m,-lmax,lmax}];
AppendTo[majorArray,GraphicsRow[minorArray,Spacings->-40,ImageSize->{100*(lmax*2+1),70}]];
,{l,0,lmax}];

GraphicsColumn[majorArray]

lmax = 1;
majorArray = {};
Do[
minorArray = {};
Do[
AppendTo[minorArray,plotOrbital[l,m]];
,{m,-lmax,lmax}];
AppendTo[majorArray,GraphicsRow[minorArray,Spacings->-40,ImageSize->{100*(lmax*2+1),70}]];
minorArray = {};
Do[
AppendTo[minorArray,plotOrbitalUp[l,m]];
,{m,-lmax,lmax}];
AppendTo[majorArray,GraphicsRow[minorArray,Spacings->-40,ImageSize->{100*(lmax*2+1),70}]];
minorArray = {};
Do[
AppendTo[minorArray,plotOrbitalDown[l,m]];
,{m,-lmax,lmax}];
AppendTo[majorArray,GraphicsRow[minorArray,Spacings->-40,ImageSize->{100*(lmax*2+1),70}]];
,{l,0,lmax}];

GraphicsColumn[majorArray]

lmax = 4;
majorArray = {};
Do[
minorArray = {};
Do[
AppendTo[minorArray,plotOrbital[l,m]];
,{m,-lmax,lmax}];
AppendTo[majorArray,GraphicsRow[minorArray,Spacings->-40,ImageSize->{100*(lmax*2+1),70}]];
,{l,0,lmax}];
GraphicsColumn[majorArray]

(* parameterized, good colors *)
l=2;
m=2;
plotOrbital[l,m]

(* parameterized, good colors *)
l=2;
m=2;
plotPoints = 50;
plotRange = .8;
psi = SphericalHarmonicY[l,m,Θ,ϕ];
psi
plot1 = SphericalPlot3D[Max[Re[psi],0],{Θ,0,π},{ϕ,0,2 π},ColorFunction->Function[{x,y,z,Θ,ϕ,r},RGBColor[1,0,0]], Mesh->None, PlotPoints->plotPoints, Boxed->False, Axes->None, PlotRange->All];
plot2 = SphericalPlot3D[-Min[Re[psi],0],{Θ,0,π},{ϕ,0,2 π},ColorFunction->Function[{x,y,z,Θ,ϕ,r},RGBColor[0,1,1]], Mesh->None, PlotPoints->plotPoints, Boxed->False, Axes->None, PlotRange->All];
Show[plot1,plot2,ViewPoint->{1.3, -2.4, 2.}]

(* parameterized, good colors *)
l=1;
m=1;
plotPoints = 50;
plotRange = .8;
psi = SphericalHarmonicY[l,m,Θ,ϕ];
psi
plot1 = SphericalPlot3D[Max[Re[psi],0],{Θ,0,π},{ϕ,0,2 π},ColorFunction->Function[{x,y,z,Θ,ϕ,r},RGBColor[1,0,0]], Mesh->None, PlotPoints->plotPoints, Boxed->False, Axes->None, PlotRange->All, ViewPoint->{1.3, -2.4, 2.},ViewVertical->{0,0,1}];
plot2 = SphericalPlot3D[-Min[Re[psi],0],{Θ,0,π},{ϕ,0,2 π},ColorFunction->Function[{x,y,z,Θ,ϕ,r},RGBColor[0,1,1]], Mesh->None, PlotPoints->plotPoints, Boxed->False, Axes->None, PlotRange->All, ViewPoint->{1.3, -2.4, 2.},ViewVertical->{0,0,1}];
Show[plot1,plot2,Boxed->True,Axes->True ]

l=0;
m=0;
plotPoints = 50;
plotRange = .8;
psi = SphericalHarmonicY[l,m,Θ,ϕ];
psi
plot1 = SphericalPlot3D[Max[Re[psi],0],{Θ,0,π},{ϕ,0,2 π},ColorFunction->Function[{x,y,z,Θ,ϕ,r},RGBColor[1,0,0]], Mesh->None, PlotPoints->plotPoints, Boxed->False, Axes->None, PlotRange->All];
plot2 = SphericalPlot3D[-Min[Re[psi],0],{Θ,0,π},{ϕ,0,2 π},ColorFunction->Function[{x,y,z,Θ,ϕ,r},RGBColor[0,1,1]], Mesh->None, PlotPoints->plotPoints, Boxed->False, Axes->None, PlotRange->All];
Show[plot1,plot2]

l=2;
m=1;
psi = SphericalHarmonicY[l,m,Θ,ϕ];
psi
SphericalPlot3D[Re[SphericalHarmonicY[l,m,Θ,ϕ]],{Θ,0,π},{ϕ,0,2 π},ColorFunction->Function[{x,y,z,Θ,ϕ,r},RGBColor[(Sign[Re[SphericalHarmonicY[l,m,Θ,ϕ]]]+1)/2,(-Sign[Re[SphericalHarmonicY[l,m,Θ,ϕ]]]+1)/2,(-Sign[Re[SphericalHarmonicY[l,m,Θ,ϕ]]]+1)/2]],ColorFunctionScaling->False,Mesh->None, PlotPoints->50, Boxed->False, Axes->None, PlotRange->All]

l=2;
m=1;
SphericalPlot3D[Re[SphericalHarmonicY[l,m,Θ,ϕ]],{Θ,0,π},{ϕ,0,2 π},ColorFunction->Function[{x,y,z,Θ,ϕ,r},RGBColor[(Sign[Re[SphericalHarmonicY[l,m,Θ,ϕ]]]+1)/2,(-Sign[Re[SphericalHarmonicY[l,m,Θ,ϕ]]]+1)/2,(-Sign[Re[SphericalHarmonicY[l,m,Θ,ϕ]]]+1)/2]],ColorFunctionScaling->False,Mesh->None, PlotPoints->50, Boxed->False, Axes->None, PlotRange->All]

l=2;
m=1;
colorsph[Θ_,ϕ_]:=RGBColor[(Sign[Re[SphericalHarmonicY[l,m,Θ,ϕ]]]+1)/2,(-Sign[Re[SphericalHarmonicY[l,m,Θ,ϕ]]]+1)/2,(-Sign[Re[SphericalHarmonicY[l,m,Θ,ϕ]]]+1)/2];
SphericalPlot3D[Re[SphericalHarmonicY[l,m,Θ,ϕ]],{Θ,0,π},{ϕ,0,2 π},ColorFunction->Function[{x,y,z,Θ,ϕ,r},colorsph[Θ,ϕ]],ColorFunctionScaling->False,Mesh->None, PlotPoints->50, Boxed->False, Axes->None, PlotRange->All]

n=2;
l=2;
m = -1;
psi[θ_,ϕ_] := SphericalHarmonicY[l, m, θ,ϕ];
ψSchrodingerAngular[{n_, l_, ml_}, r_, θ_,ϕ_] := SphericalHarmonicY[l, ml, θ,ϕ];
colorsph[Θ_,ϕ_] := RGBColor[(Sign[Re[SphericalHarmonicY[l, m, θ,ϕ]]]+1)/2,0,(-Sign[Re[SphericalHarmonicY[l, m, θ,ϕ]]]+1)/2]; 
SphericalPlot3D[Re[SphericalHarmonicY[l, m, θ,ϕ]],{θ,0,Pi},{ϕ,0,2 Pi},PlotRange->All,PlotPoints->50,Mesh->None, Boxed->False, Axes->None, ColorFunction->Function[{x,y,z,Θ,ϕ,r},colorsph[Θ,ϕ]],ColorFunctionScaling->False]

l=2;
m=-1;
psi = SphericalHarmonicY[l, m, θ,ϕ];
colorsph[Θ_,ϕ_]:=RGBColor[(Sign[Re[psi]]+1)/2,0,(-Sign[Re[psi]]+1)/2];
colorsph[Θ,ϕ]
SphericalPlot3D[Re[psi],{Θ,0,π},{ϕ,0,2 π},ColorFunction->Function[{x,y,z,Θ,ϕ,r},colorsph[Θ,ϕ]],ColorFunctionScaling->False,Mesh->None,PlotPoints->130,PlotPoints->50, Boxed->False, Axes->None]

l=2;
m=1;
colorsph[Θ_,ϕ_]:=RGBColor[(Sign[Re[SphericalHarmonicY[l,m,Θ,ϕ]]]+1)/2,(-Sign[Re[SphericalHarmonicY[l,m,Θ,ϕ]]]+1)/2,(-Sign[Re[SphericalHarmonicY[l,m,Θ,ϕ]]]+1)/2];
SphericalPlot3D[Re[SphericalHarmonicY[l,m,Θ,ϕ]],{Θ,0,π},{ϕ,0,2 π},ColorFunction->Function[{x,y,z,Θ,ϕ,r},colorsph[Θ,ϕ]],ColorFunctionScaling->False,Mesh->None,PlotPoints->230, Boxed->False, Axes->None, PlotRange->All]

l=2;
m=-1;
colorsph[Θ_,ϕ_]:=RGBColor[(Sign[Re[SphericalHarmonicY[l,m,Θ,ϕ]]]+1)/2,0,(-Sign[Re[SphericalHarmonicY[l,m,Θ,ϕ]]]+1)/2];
SphericalPlot3D[Re[SphericalHarmonicY[l,m,Θ,ϕ]],{Θ,0,π},{ϕ,0,2 π},ColorFunction->Function[{x,y,z,Θ,ϕ,r},colorsph[Θ,ϕ]],ColorFunctionScaling->False,Mesh->None,PlotPoints->130,PlotPoints->50, Boxed->False, Axes->None]

l=1;
m=0;
colorsph[Θ_,Φ_]:=RGBColor[(Sign[Re[SphericalHarmonicY[l,m,Θ,Φ]]]+1)/2,0,(-Sign[Re[SphericalHarmonicY[l,m,Θ,Φ]]]+1)/2];
SphericalPlot3D[Re[SphericalHarmonicY[l,m,Θ,Φ]],{Θ,0,π},{Φ,0,2 π},ColorFunction->Function[{x,y,z,Θ,Φ,r},colorsph[Θ,Φ]],ColorFunctionScaling->False,Mesh->None,PlotPoints->130,PlotPoints->50, Boxed->False, Axes->None]

colorsph[Θ_,Φ_]:=RGBColor[(Sign[Re[SphericalHarmonicY[2,1,Θ,Φ]]]+1)/2,0,(-Sign[Re[SphericalHarmonicY[2,1,Θ,Φ]]]+1)/2];
SphericalPlot3D[Re[SphericalHarmonicY[2,1,Θ,Φ]],{Θ,0,π},{Φ,0,2 π},ColorFunction->Function[{x,y,z,Θ,Φ,r},colorsph[Θ,Φ]],ColorFunctionScaling->False,Mesh->None,PlotPoints->130]