z=16.6;
bh=2;
bz=1.02*bh/3;
ny=10;
for n=0:ny;
   dh(n+1)=1.8;
   dz(n+1)=n/1;
   d(n+1)=sqrt(dh(n+1)^2+dz(n+1)^2);
   fh(n+1)=exp(-bh*(d(n+1)/z)^1.33);
   fz(n+1)=exp(-bz*(d(n+1)/z)^1.33);
   ff(n+1)=(dh(n+1)^2*fh(n+1)+dz(n+1)^2*fz(n+1))/d(n+1)^2;
end;
plot(dz(1:ny),ff(1:ny))