Question: Why Want My Graph Plot, Please Help

This is my program below. When I try to get an output, nothing happens, please help!

> restart:

> with(plots):

readlib(readdata):

lrc:=array(-1..1,-1..1):

lrn:=array(-1..1,-1..1):

urn:=array(-2..2,-2..2):

unr:=array(-2..2,-2..2):

urc:=array(-2..2,-2..2):

Ylm:=array(1..2,-2..2):

Dj:=array(-1..1,-1..1):

Dj2:=array(-2..2,-2..2):

rhoc[-1..1,-1..1,1..21]:

rhocu[-2..2,-2..2,1..21]:

> assume(t,real):

assume(phi,real):

assume(p,real):

assume(m,integer):

> Ylm[1,1]:=-sqrt(3/8/Pi)*sin(phi)*exp(I*t):

Ylm[1,0]:=-sqrt(3/4/Pi)*cos(phi):

Ylm[1,-1]:=sqrt(3/8/Pi)*sin(phi)*exp(-I*t):

Ylm[2,2]:=sqrt(15/32/Pi)*sin(phi)^2*exp(2*I*t):

Ylm[2,1]:=-sqrt(15/8/Pi)*sin(phi)*cos(phi)*exp(I*t):

Ylm[2,0]:=sqrt(5/16/Pi)*(3*cos(phi)^2-1):

Ylm[2,-1]:=sqrt(15/8/Pi)*sin(phi)*cos(phi)*exp(-I*t):

Ylm[2,-2]:=sqrt(15/32/Pi)*sin(phi)^2*exp(-2*I*t):

> for j from -1 to 1 do

for k from -1 to 1 do

for l from 1 to 21 do

rhoc[j,k,l]:=0;

od

od

od:

for j from -2 to 2 do

for k from -2 to 2 do

for l from 1 to 21 do

rhocu[j,k,l]:=0;

od

od

od:

> datafile:=readdata(`HgL100mW0det.dat`,6):

for j from 1 to 21 do

timescale[j]:=op(1,datafile[j]):

rho22[j]:=op(2,datafile[j]):

rho33[j]:=op(3,datafile[j]):

rho31[j]:=op(4,datafile[j]):

rho32R[j]:=op(5,datafile[j]):

rho32I[j]:=op(6,datafile[j]):

od:

datafile:=readdata(`HgU100mW0det.dat`,6):

for j from 1 to 21 do

timescale[j]:=op(1,datafile[j]):

rho55[j]:=op(2,datafile[j]):

rho66[j]:=op(3,datafile[j]):

rho64[j]:=op(4,datafile[j]):

rho65R[j]:=op(5,datafile[j]):

rho65I[j]:=op(6,datafile[j]):

od:

> for j from 1 to 21 do

rhoc[1,1,j]:=rho33[j]:

rhoc[1,0,j]:=rho32R[j]+I*rho32I[j]:

rhoc[1,-1,j]:=rho31[j]:

rhoc[0,1,j]:=rho32R[j]-I*rho32I[j]:

rhoc[0,0,j]:=rho22[j]:

rhoc[0,-1,j]:=-1*rhoc[0,1,j]:

rhoc[-1,1,j]:=rhoc[1,-1,j]:

rhoc[-1,0,j]:=-1*rhoc[1,0,j]:

rhoc[-1,-1,j]:=rho33[j]:

rhocu[1,1,j]:=rho66[j]:

rhocu[1,0,j]:=rho65R[j]+I*rho65I[j]:

rhocu[1,-1,j]:=rho64[j]:

rhocu[0,1,j]:=rho65R[j]-I*rho65I[j]:

rhocu[0,0,j]:=rho55[j]:

rhocu[0,-1,j]:=-1*rhocu[0,1,j]:

rhocu[-1,1,j]:=rhocu[1,-1,j]:

rhocu[-1,0,j]:=-1*rhocu[1,0,j]:

rhocu[-1,-1,j]:=rho66[j]:

od:

> for j from 1 to 21 do

lower[j]:=cat(`Hglower0det`,j,`.ps`);

upper[j]:=cat(`Hgupper0det`,j,`.ps`);

od:

> for j from 1 by 1 to 21 do

nWcol[j]:=

combine(evalc(sum(sum(rhoc[jl,n,j]*Ylm[1,jl]*conjugate(Ylm[1,n]),jl=-1..1),n=-1..1)),trig);

nWcolu[j]:=

combine(evalc(sum(sum(rhocu[jl,n,j]*Ylm[2,jl]*conjugate(Ylm[2,n]),jl=-2..2),n=-2..2)),trig);

od:

> for j from 1 by 1 to 21 do

plotsetup(ps,plotoutput=lower[j],plotoptions=`,color,portrait,noborder`):

sphereplot(((nWcol[j]),t=-Pi..Pi,phi=0..Pi),grid=[30,30],axes=BOXED,

orientation=[-60,70],shading=Z,scaling=CONSTRAINED,tickmarks=[3,3,3]);

plotsetup(ps,plotoutput=upper[j],plotoptions=`,color,portrait,noborder`):

sphereplot(((nWcolu[j]),t=-Pi..Pi,phi=0..Pi),grid=[30,30],axes=BOXED,

orientation=[-60,70],shading=Z,scaling=CONSTRAINED,tickmarks=[3,3,3]);

od;

> sphereplot(((nWcol[8]),t=-Pi..Pi,phi=0..Pi),grid=[30,30],axes=BOXED,

orientation=[-60,70],shading=Z,scaling=CONSTRAINED);

sphereplot(((nWcolu[8]),t=-Pi..Pi,phi=0..Pi),grid=[30,30],axes=BOXED,

orientation=[-60,70],shading=Z,scaling=CONSTRAINED);

>