% example_2D_scatter.m

%

% This script file demonstrates, in 2D,

% how "scattering" information can be

% used to probe spatial structure.

%

% K. Beers. MIT ChE. 12/19/2002

% First, we generate a grid

% x-direction

Px = pi; Nx = 2^8; dx = 2*Px/Nx;

x_val = [0 : dx : (2*Px - dx)];

% y-direction

Py = 1.5*pi; Ny = 2^8; dy = 2*Py/Ny;

y_val = [0 : dy : (2*Py - dy)];

% use meshgrid to produce X and Y

% values at 2-D grid points.  Here

% the (i,j) element of the matrix

% captures the value of a function

% at the position (x(i),y(j))

[X,Y] = meshgrid(x_val,y_val);

% Next, we generate a spatial distribution

% function of scattering sites.  We place

% spherical site densities around a corner

% lattice site and a middle lattice site.

RHO = zeros(size(X));

sigma_sq_center = 0.5; sigma_sq_corner = 0.1;

dr = zeros(1,2);

weight_center = 1.5; weight_corner = 1;

for kx=1:Nx

    for ky=1:Ny

        % center site

        dr(1) = X(kx,ky) - Px;

        dr(2) = Y(kx,ky) - Py;

        dist_sq = dot(dr,dr);

        RHO(kx,ky) = RHO(kx,ky) + weight_center* ...

            exp(-dist_sq/2/sigma_sq_center);

       % corner site

       dr(1) = X(kx,ky) - (2*Px)*round(X(kx,ky)/(2*Px));

       dr(2) = Y(kx,ky) - (2*Py)*round(Y(kx,ky)/(2*Py));

       dist_sq = dot(dr,dr);

       RHO(kx,ky) = RHO(kx,ky) + weight_corner* ...

           exp(-dist_sq/2/sigma_sq_corner);

   end

end

% plot the distribution of sites

figure; contourf(X,Y,RHO);

axis equal; colorbar;

xlabel('x'); ylabel('y');

title('\rho(x,y)');

% Now, we generate the discrete Fourier transform.

% Calculate the sampled frequencies

% x-direction

dqx = pi/Px; qx_crit = pi/dx;

qx_val = [0 : dqx : (2*qx_crit-dqx)];

% y-direction

dqy = pi/Py; qy_crit = pi/dy;

qy_val = [0 : dqy : (2*qy_crit-dqy)];

[QX,QY] = meshgrid(qx_val,qy_val);

% Calculate the Fourier transform and power spectrum

% (the autocorrelation) of the site distribution function

RHO_FT = fft2(RHO);

RHO_PS = real(conj(RHO_FT).*RHO_FT);

% plot the power spectrum

figure; surf(QX,QY,RHO_PS);

xlabel('q_x'); ylabel('q_y');

zlabel('|\rho(q)|^2');

title('Power spectrum of \rho(x,y)');

% contour plot

figure;

z_max = max(max(log10(RHO_PS)));

color_vect = (ceil(z_max) : -0.5 : 0);

contourf(QX,QY,log10(RHO_PS),color_vect);

axis([qx_val(1) 10 qy_val(1) 10]);

xlabel('q_x'); ylabel('q_y');

title('log10 of Power spectrum of \rho(x,y)');

colorbar;

% get Patterson diagram

RHO_Patterson = abs(ifft2(RHO_PS));

figure; contourf(X,Y,RHO_Patterson,15);

axis equal; colorbar;

xlabel('s_x'); ylabel('s_y');

title('Patterson diagram');