MatLab R207b	Requirements to PC: operative memory 2GB.
-----------------------------------------------------------------
1) Auxiliary functions:
--------------------------------------------------------------
function Z = xBesselj(n,X)

% function Z = dxBesselj(X,n) computes the product X^1/2 on a cylindrical 
% function J which has argument X (can be 3D-arrey) and the order n+1/2 (n=1,2,3,4,5).

M = ones(size(X));	k=(2/pi)^(1/2);

if n==1
    Z = k*( sin(X)./X - cos(X) );
elseif n==2
    Z = k*( (3*(X.^(-2))-M).*sin(X) -3*cos(X)./X );
elseif n==3
    Z = k*( (15*(X.^(-3))-6*(X.^(-1))).*sin(X) - (15*(X.^-(2))-M).*cos(X) );
elseif n==4
    Z = k*( (105*(X.^(-4))-45*(X.^(-2))+M).*sin(X) - (105*(X.^(-3))-10*(X.^(-1))).*cos(X) );
end
---------------------------------------------------------------
function y = turn3d(x,dim)

[i,j,k]=size(x);

if dim==1
    n=1:i;    y = x(ones(1,i)*(i+1)-n,:,:);
elseif dim==2
    n=1:j;    y = x(:,ones(1,j)*(j+1)-n,:);
elseif dim==3
    n=1:k;    y = x(:,:,ones(1,k)*(k+1)-n);
end
-----------------------------------------------------------------
2) Script
-----------------------------------------------------------------
r0 = 5; b1= 20; d1 = 20; m=260; n=2.2;

X = cumsum(ones(n*m,m,m),1)-(n*m/2 + 0.5)*ones(n*m,m,m);
Y = cumsum(ones(n*m,m,m),2)-(m/2+0.5)*ones(n*m,m,m);
Z = cumsum(ones(n*m,m,m),3)-(m/2+0.5)*ones(n*m,m,m);
save tmp/U X Y Z;	clear X Y Z;

for x = 80:625;
    r = x/5;
        load tmp/U;
        R1 = ((X-(r+b1/2)*ones(n*m,m,m)).^2 + (Y-(d1/2)*ones(n*m,m,m)).^2 + Z.^2).^(1/2);
        R2 = ((X-(r-b1/2)*ones(n*m,m,m)).^2 + (Y-(d1/2)*ones(n*m,m,m)).^2 + Z.^2).^(1/2);
        clear X Y Z;
        G = 0.25*( sign( R1 - r0*ones(n*m,m,m) ) +  ones(n*m,m,m) ) ...
            .*( sign( R2 - r0*ones(n*m,m,m) ) +  ones(n*m,m,m) );
        G = G.*turn3d(G,1);         G = G.*turn3d(G,2);
        G = G.*permute(G,[1,3,2]);
        W = xBesselj(1,R1);         clear R1 ;       
        W = W + xBesselj(1,R2);     clear R2;
        W = W + turn3d(W,1);        W = W + turn3d(W,2);
        W = G.*( W + permute(W,[1 3 2]) ).^2;   clear G;
        sW(x) = sum(sum(sum(W)))/(2*m^3);       clear W;
        save tmp/W ;
end
-----------------------------------------------------------------------------
3) Plot
-----------------------------------------------------------------------------
n=size(sW,2);	x=2*(1:n)/5;
F= diff(sW); 	F=cat(2,F(1,1),F); 	F(1,79)=0;  F(1,80)=0;  
[b,a]=butter(5,0.02);		Ff= filter(b,a,F);

plot(x,sW,'-',x,20*F,'-',x,400*Ff,'-',x,zeros(1,n),':');  % Plots
legend('W','20*F','400*Ff');