one way to plot 3D hydrogen atom orbitals

function create()
  variable n=5,l=3,m=0
  variable size1=400,size2=400
  wave w=calcparametric(N,L,M,size1,size2)
  execute "Gizmo0()"
end
 
function/wave calcparametric(N,L,M,size1,size2)
  variable N,L,M,size1,size2
  make/d/o/n=((size1)^2*(size2)) wavex,wavey,wavez,ff
  variable i,j,k,dt,df,psi2,rr,theta,phi,t=0
  variable xx,yy,zz  
  dt=pi/size1
  df=2*dt
  for(i=0;i<size1;i+=1)
    theta=i*dt
    for(j=0;j<size1;j+=1)
      phi=j*df
      for(k=0;k<size2;k+=1)
        rr=k*0.05      
        variable temp1=factorial(n+l),temp2=factorial(n-l-1),temp3=laguerrea(n-l-1,2*l+1,2*rr/n/0.529)
        psi2=magsqr(sqrt((2/n/0.529)^3/2/n/(temp1)^3*temp2)*exp(-rr/n/0.529)*(2*rr/n/0.529)^l*temp3*sphericalHarmonics(l,m,theta,phi))
        ff[t]=psi2
        wavex[t]=rr*sin(theta)*cos(phi)
        wavey[t]=rr*sin(theta)*sin(phi)
        wavez[t]=rr*cos(theta)
        t+=1
      endfor
    endfor
  endfor
  concatenate/np=1 {wavex,wavey,wavez,ff},dest
 
  Variable pee, qew, kappa
  variable u,v,z
  make/n=((size2)/10,(size2)/10,(size2)/10)/d/o kcube=0
  setscale/p x,-0.05*size2,1,kcube
  setscale/p y,-0.05*size2,1,kcube
  setscale/p z,-0.05*size2,1,kcube
  duplicate/o kcube, BlueJacobian
  BlueJacobian = 0
  for (j=0; j<DimSize(dest,0); j+=1)
            pee = round((dest[j][0] - DimOffset(kCube,0))/DimDelta(kCube,0))
            qew = round((dest[j][1] - DimOffset(kCube,1))/DimDelta(kCube,1))
            kappa=round((dest[j][2] - DimOffset(kCube,2))/DimDelta(kCube,2))
            if (pee >= 0 && pee < DimSize(kCube,0) && qew >= 0 && qew < DimSize(kCube,1)&& kappa >= 0 && kappa < DimSize(kCube,2))
                kCube[pee][qew][kappa] += dest[j][3]
                BlueJacobian[pee][qew][kappa] += 1
            endif
  endfor
  for(u=0;u<DimSize(kCube,0);u+=1)
       for(v=0;v<DimSize(kCube,1);v+=1)
            for(z=0;z<dimsize(kcube,2);z+=1)
                kCube[u][v][z] /= max(1,BlueJacobian[u][v][z])
            endfor
        endfor
  endfor
  killwaves dest,BlueJacobian,wavex,wavey,wavez,ff
  return kcube
end
 
Window Gizmo0() : GizmoPlot
    PauseUpdate; Silent 1       // building window...
    // Building Gizmo 8 window...
    NewGizmo/W=(21,44,330,320)
    ModifyGizmo startRecMacro=700
    ModifyGizmo scalingOption=63
    AppendToGizmo isoSurface=root:kcube,name=isoSurface0
    ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ surfaceColorType,1}
    ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ lineColorType,0}
    ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ lineWidthType,0}
    ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ fillMode,2}
    ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ lineWidth,1}
    ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ isoValue,1e-08}
    ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ frontColor,1,0,0,1}
    ModifyGizmo ModifyObject=isoSurface0,objectType=isoSurface,property={ backColor,0,0,1,1}
    ModifyGizmo setDisplayList=0, object=isoSurface0
    ModifyGizmo autoscaling=1
    ModifyGizmo currentGroupObject=""
    ModifyGizmo showInfo
    ModifyGizmo infoWindow={559,7,1392,321}
    ModifyGizmo endRecMacro
    ModifyGizmo SETQUATERNION={-0.681597,0.118549,-0.123937,0.711347}
EndMacro