clear
clc
l=3.2; %l is wavelength of the acoustic waves
s=90*pi/180;% angle between wave vectors and the Z direction
s1=0*pi/180; %angle between wave vectors and the x-y plane
a=[1,1,1,1,1,1,1];
t1=sin(s);
t2=cos(s);
t3=5*t1/5;  

ab=6; % number of wave vectors
n=360/ab; % angles between each wave vector

k0=[t1*cos(s1+0*n*pi/180),t1*sin(s1+0*n*pi/180),t2]*2*pi/l;
k1=[t3*cos(s1+1*n*pi/180),t3*sin(s1+1*n*pi/180),t2]*2*pi/(l);
k2=[t1*cos(s1+2*n*pi/180),t1*sin(s1+2*n*pi/180),t2]*2*pi/l; 
k3=[t3*cos(s1+3*n*pi/180),t3*sin(s1+3*n*pi/180),t2]*2*pi/(l);
k4=[t1*cos(s1+4*n*pi/180),t1*sin(s1+4*n*pi/180),t2]*2*pi/(l); 
k5=[t3*cos(s1+5*n*pi/180),t3*sin(s1+5*n*pi/180),t2]*2*pi/(l);

% the following is the amplitude of every wave vector
a0=1;
a1=1;
a2=1;
a3=1;
a4=1;
a5=1;

% the following is setting the simulation area
X=[-12:0.005:12];   
Y=-[-12:0.005:12]';
Z=0;
x=ones(size(X,2),1)*X;
y=Y*ones(1,size(Y,1));
z=Z*ones(size(Y,1),size(X,2));


disp('start simulation......');
% the following is the simulation of the wavefield distribution.
% phase relations between wave vectors should be modulated here
   
        I=a1^2+a2^2+a3^2+a4^2+a0^2+a5^2 ...
          +2*(a0*a1*cos((k0(1)-k1(1))*x+(k0(2)-k1(2))*y+(k0(3)-k1(3))*z+0*pi/180)...
          +a0*a2*cos((k0(1)-k2(1))*x+(k0(2)-k2(2))*y+(k0(3)-k2(3))*z+0*pi/180)...
          +a0*a3*cos((k0(1)-k3(1))*x+(k0(2)-k3(2))*y+(k0(3)-k3(3))*z+0*pi/180)...
          +a0*a4*cos((k0(1)-k4(1))*x+(k0(2)-k4(2))*y+(k0(3)-k4(3))*z+0*pi/180)...
          +a5*a2*cos((k5(1)-k2(1))*x+(k5(2)-k2(2))*y+(k5(3)-k2(3))*z+0*pi/180)...
          +a5*a1*cos((k5(1)-k1(1))*x+(k5(2)-k1(2))*y+(k5(3)-k1(3))*z+0*pi/180)...
          +a5*a0*cos((k5(1)-k0(1))*x+(k5(2)-k0(2))*y+(k5(3)-k0(3))*z+0*pi/180)...
          +a5*a3*cos((k5(1)-k3(1))*x+(k5(2)-k3(2))*y+(k5(3)-k3(3))*z+0*pi/180)...
          +a5*a4*cos((k5(1)-k4(1))*x+(k5(2)-k4(2))*y+(k5(3)-k4(3))*z+0*pi/180)...
          +a4*a1*cos((k4(1)-k1(1))*x+(k4(2)-k1(2))*y+(k4(3)-k1(3))*z+0*pi/180)...
          +a4*a2*cos((k4(1)-k2(1))*x+(k4(2)-k2(2))*y+(k4(3)-k2(3))*z+0*pi/180)...
          +a4*a3*cos((k4(1)-k3(1))*x+(k4(2)-k3(2))*y+(k4(3)-k3(3))*z+0*pi/180)...
          +a1*a2*cos((k1(1)-k2(1))*x+(k1(2)-k2(2))*y+(k1(3)-k2(3))*z+0*pi/180)...
          +a1*a3*cos((k1(1)-k3(1))*x+(k1(2)-k3(2))*y+(k1(3)-k3(3))*z+0*pi/180)...
          +a2*a3*cos((k2(1)-k3(1))*x+(k2(2)-k3(2))*y+(k2(3)-k3(3))*z+0*pi/180));
     
    
disp ('ending simulation and plottingˇ­')
  

figure
  pcolor (I)


%imshow( I,[0 max(I(:))])
 %mesh (x1, y1, I)
shading interp
  axis off
  axis equal




% [DX,DY] = gradient(I,0.5,6);
%  
% contour(X,Y,I)
%  
%  hold on
%  
%  quiver(X,Y,DX,DY)
%  
%  colormap hsv
%  
%  hold off
%   
% imshow( I,[0 10])
% 
% imshow( I,[0 max(I(:))])
%  colorbar('horiz')