How to derive beam width and beam divergence from Fourier transform of initial wavefront

Hi, I wonder, how the divergence of the laser beam with non-gaussian profile can be computed. Namely, divergence of the product of such beam combining:

One can assume wavefront from each subaperture as flat. To obtain the diffraction pattern I try

ms=0.5;
[mm nn]=meshgrid(-400:ms:400);
pupil=exp(1i*59052493488.5*100.2)*double(sqrt((mm).^2+(nn-6).^2)<=1.5|sqrt((mm).^2+(nn-12).^2)<=1.5|sqrt((mm).^2+(nn+6).^2)<=1.5|sqrt((mm).^2+(nn+12).^2)<=1.5|sqrt((mm-12).^2+(nn-6).^2)<=1.5|sqrt((mm-12).^2+(nn+6).^2)<=1.5|sqrt((mm+12).^2+(nn+6).^2)<=1.5|sqrt((mm+12).^2+(nn-6).^2)<=1.5|sqrt((mm-12).^2+(nn).^2)<=1.5|sqrt((mm+12).^2+(nn).^2)<=1.5|sqrt((mm).^2+(nn).^2)<=1.5|sqrt((mm-6).^2+(nn+3).^2)<=1.5|sqrt((mm-6).^2+(nn+9).^2)<=1.5|sqrt((mm-6).^2+(nn-3).^2)<=1.5|sqrt((mm-6).^2+(nn-9).^2)<=1.5|sqrt((mm+6).^2+(nn-9).^2)<=1.5|sqrt((mm+6).^2+(nn-3).^2)<=1.5|sqrt((mm+6).^2+(nn+9).^2)<=1.5|sqrt((mm+6).^2+(nn+3).^2)<=1.5);
       pupil2=ifftshift(pupil);
       psf2=fftshift(ifft2(pupil2));
           pupil3=fftshift(fft2(psf2));
           psf3=ifft2(ifftshift(pupil3));
figure(1);
imagesc(abs(psf3));
figure(2);
imagesc(imag(psf3));

abs(psf3) don't depend on distance of propagation, but phases pattern is changing due to different values in complex exponential multiplier. I need to know the divergence of central peak. I suppose it can't be described as gaussian beam.

rp-photonics.com

article about non gaussian beam divergence

I've read that divergence can be retrieved from root mean squared width of fourier spectrum, but I don't exactly understand the relation between width of fourier spectrum and width of the beam.

R = sqrt(sum(sum(imag(psf3).^2))/prod(size(psf3)));

I assumed root mean squared width with the formula above. Probably it should be applied to cropped image of diffraction phase pattern which contains central peak. Maybe root mean squared value need to be found for image with one phase portrait subtracted from image with another phase portrait, and I would use that value further. For example first image is the fourier transform with the initial phase multiplier exp(1i*101), and another with exp(1i*102). But I don't know if it is right or not, and don't understand how connect fourier transform with beam diameter at all. Help me, please.