rh=0.3;
rs=0.25;
o=1.21;
c=340;
f=600;
a=pi/30;
k=2*pi*f;
N=60;
dt=2*pi/N;
df=2*pi/N;
sumprs=0;
for n=0:N
    for m=-n:n
        sumPrh=0;
for t=-pi

i/30

i
    for f=-pi

i/30

i
         x=rh*sin(f)*cos(t);
         y=rh*sin(f)*sin(t);
         z=rh*cos(f);
         r1=sqrt((x+0.2)^2+y^2+z^2);
         r2=sqrt(x^2+(y-0.2)^2+z^2);
         p1=-a*i*2*pi*f*o*a*exp(i*k*(r1-a))/((1-i*k*a)*r1);
         p2=-a*i*2*pi*f*o*a*exp(i*k*(r2-a))/((1-i*k*a)*r2);
         p=p1+p2; 
         L=legendre(n,cos(t));                   
         G=sqrt((2*n+1)*factorial(n-m)/(4*pi*factorial(n+m)));
         Y=G*L*exp(i*m*f);
         Prh=p*conj(Y)*sin(t)*dt*df;
         sumPrh=sumPrh+Prh;
    end
end
H=besselh(n,k*rs)/besselh(n,k*rh);
t1=-pi

i/30

i;
f1=-pi

i/30

i;
[tt,ff]=meshgrid(t1,f1);
xx=rh*sin(ff).*cos(tt);
yy=rh*sin(ff).*sin(tt);
zz=rh*cos(ff);
L1=legendre(n,cos(tt));                   
G1=sqrt((2*n+1)*factorial(n-m)/(4*pi*factorial(n+m)));
Y1=G1*L1.*exp(i*m*ff);
prs=H*sumPrh*Y1;
sumprs=sumprs+prs;
    end
end
h=abs(sumprs);
surface(xx,yy,zz,h)
??? Error using ==> times
Number of array dimensions must match for binary array op.