%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This function numerically computes the stress field for a collection of 
% arrays of dislocations with the specified Burgers vectors and positions.  
%
% Arguments:
% - X, Y:                    computational grids (as generated by meshgrid())
% - num_dislocation_arrays:  number of dislocation arrays
% - burgers_vectors:         matrix containing Burgers vectors for 
%                            dislocation arrays.  The i-th row should
%                            contain the Burgers vector for the i-th
%                            dislocation array
% - positions:               matrix containing positions for dislocation
%                            arrays.  The i-th row should contain the (x,y)
%                            coordinates of the i-th dislocation array
% - numerical_core_radius:   radius to use for numerical core
% - delta_function_type:     type of delta function to use when smearing
%                            core.  Choices are:
%                              1: first-order trigonometric approximation
%                                 to delta function (see Osher & Fedkiw 2003)
%                              2: second-order trigonometric approximation
%                                 to delta function 
%                              3: second-order trigonometric approximation
%                                 to delta function 
%                              4: modified Lorentzian approximation
%                                 to delta function 
% - G:                       shear modulus
% - poisson_ratio:           Poisson ratio
%
% Return values:
% - sigma_xx:                xx-component of stress field
% - sigma_yy:                yy-component of stress field
% - sigma_xy:                xy-component of stress field
% - d_z:                     dislocation line field for last 
%                            dislocation line array 
%
%
% NOTES:
% - The line direction of all dislocations is assumed to be in the (0,0,1) 
%   direction.  
% - Each array is oriented such that the normal to the plane of the 
%   dislocation array is in the (1,0,0) direction and all dislocation 
%   arrays are periodically repeated in the x-direction.  
% - The stress fields for the dislocation configuration are computed using
%   the Fourier transform formulation in Xiang et. al.  
% - The grid is assumed to be uniform in each spatial coordinate direction.
%
% Kevin T. Chu
% MAE, Princeton University
% 08/2007
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [sigma_xx, sigma_yy, sigma_xy, d_z] = ...
  compute_numerical_stress_fields(X, Y, ...
                                  num_dislocation_arrays, ...
                                  burgers_vectors, positions, ...
                                  numerical_core_radius, ...
                                  delta_function_type, ...
                                  G, poisson_ratio)

% check arguments
if (num_dislocation_arrays ~= size(burgers_vectors,1))
  error('Number of Burgers vectors does not match number of dislocation arrays');
end
if (num_dislocation_arrays ~= size(positions,1))
  error('Number of positions does not match number of dislocation arrays');
end


% compute Nx and Ny
[Nx,Ny] = size(X);

% compute Lx and Ly
Lx = Nx*(X(2,1)-X(1,1));
Ly = Ny*(Y(1,2)-Y(1,1));

% initialize solution for stress
sigma_xx = zeros(size(X));
sigma_yy = zeros(size(X));
sigma_xy = zeros(size(X));


% compute stress field by summing up the stress field from 
% individual dislocations
for line_num = 1:num_dislocation_arrays

  % get Burgers vector and position of current dislocation line
  b = burgers_vectors(line_num,:);
  pos = positions(line_num,:);

  switch (delta_function_type) 

    case 1 
      % compute delta function using first-order delta function
      delta_X = zeros(size(X));
      x_in_core = find(abs(X - pos(1)) <= numerical_core_radius); 
      delta_X(x_in_core) = 0.5/numerical_core_radius ...
        * (1+cos(pi*(X(x_in_core)-pos(1))/numerical_core_radius));

      delta_Y = zeros(size(Y));
      y_in_core = find(abs(Y - pos(2)) <= numerical_core_radius); 
      delta_Y(y_in_core) = 0.5/numerical_core_radius ...
        * (1+cos(pi*(Y(y_in_core)-pos(2))/numerical_core_radius));

    case 2
      % compute delta functions using second-order delta function
 
     delta_X = zeros(size(X));
     x_in_core = find(abs(X - pos(1)) <= numerical_core_radius); 
     delta_X(x_in_core) = -1/6/numerical_core_radius ...
        * (1+cos(pi*(X(x_in_core)-pos(1))/numerical_core_radius));
      x_in_core = find(abs(X - pos(1)) <= 0.5*numerical_core_radius); 
      delta_X(x_in_core) = 4/3/numerical_core_radius ...
        * (1+cos(2*pi*(X(x_in_core)-pos(1))/numerical_core_radius))...
                          -1/6/numerical_core_radius ...
        * (1+cos(pi*(X(x_in_core)-pos(1))/numerical_core_radius));

      delta_Y = zeros(size(Y));
      y_in_core = find(abs(Y - pos(2)) <= numerical_core_radius); 
      delta_Y(y_in_core) = -1/6/numerical_core_radius ...
        * (1+cos(pi*(Y(y_in_core)-pos(2))/numerical_core_radius));
      y_in_core = find(abs(Y - pos(2)) <= 0.5*numerical_core_radius); 
      delta_Y(y_in_core) = 4/3/numerical_core_radius ...
        * (1+cos(2*pi*(Y(y_in_core)-pos(2))/numerical_core_radius))...
                          -1/6/numerical_core_radius ...
        * (1+cos(pi*(Y(y_in_core)-pos(2))/numerical_core_radius));

    case 3
      % compute delta functions using third-order delta function
 
     delta_X = zeros(size(X));
     x_in_core = find(abs(X - pos(1)) <= numerical_core_radius); 
     delta_X(x_in_core) = 1/48/numerical_core_radius ...
        * (1+cos(pi*(X(x_in_core)-pos(1))/numerical_core_radius));
      x_in_core = find(abs(X - pos(1)) <= 0.5*numerical_core_radius); 
      delta_X(x_in_core) = -16/15/numerical_core_radius ...
        * (1+cos(2*pi*(X(x_in_core)-pos(1))/numerical_core_radius))...
                          +1/48/numerical_core_radius ...
        * (1+cos(pi*(X(x_in_core)-pos(1))/numerical_core_radius));
      x_in_core = find(abs(X - pos(1)) <= numerical_core_radius/3); 
      delta_X(x_in_core) = 243/80/numerical_core_radius ...
        * (1+cos(3*pi*(X(x_in_core)-pos(1))/numerical_core_radius))...
                           -16/15/numerical_core_radius ...
        * (1+cos(2*pi*(X(x_in_core)-pos(1))/numerical_core_radius))...
                          +1/48/numerical_core_radius ...
        * (1+cos(pi*(X(x_in_core)-pos(1))/numerical_core_radius));

     delta_Y = zeros(size(Y));
     y_in_core = find(abs(Y - pos(2)) <= numerical_core_radius); 
     delta_Y(y_in_core) = 1/48/numerical_core_radius ...
        * (1+cos(pi*(Y(y_in_core)-pos(1))/numerical_core_radius));
      y_in_core = find(abs(Y - pos(2)) <= 0.5*numerical_core_radius); 
      delta_Y(y_in_core) = -16/15/numerical_core_radius ...
        * (1+cos(2*pi*(Y(y_in_core)-pos(2))/numerical_core_radius))...
                          +1/48/numerical_core_radius ...
        * (1+cos(pi*(Y(y_in_core)-pos(2))/numerical_core_radius));
      y_in_core = find(abs(Y - pos(2)) <= numerical_core_radius/3); 
      delta_Y(y_in_core) = 243/80/numerical_core_radius ...
        * (1+cos(3*pi*(Y(y_in_core)-pos(2))/numerical_core_radius))...
                           -16/15/numerical_core_radius ...
        * (1+cos(2*pi*(Y(y_in_core)-pos(2))/numerical_core_radius))...
                          +1/48/numerical_core_radius ...
        * (1+cos(pi*(Y(y_in_core)-pos(2))/numerical_core_radius));


    case 4
      % compute delta functions using a modified Lorentzian function

      delta_X = 1/(2*(atan(pi)/pi-1/(1+pi^2))) ...
              * 1/numerical_core_radius ...
              * ( 1./(1+(pi*(X-pos(1))/numerical_core_radius).^2) ...
                - 1/(1+pi^2) );
      x_outside_core = find(abs(X - pos(1)) > numerical_core_radius); 
      delta_X(x_outside_core) = 0;
      delta_Y = 1/(2*(atan(pi)/pi-1/(1+pi^2))) ...
              * 1/numerical_core_radius ...
              * ( 1./(1+(pi*(Y-pos(2))/numerical_core_radius).^2)  ...
                - 1/(1+pi^2) );
      y_outside_core = find(abs(Y - pos(2)) > numerical_core_radius); 
      delta_Y(y_outside_core) = 0;


    otherwise
      error('Unsupported order for delta function');
  end

  % compute dislocation line field
  d_z = delta_X.*delta_Y;

  % compute FFT of d_z
  d_z_fft = fft2(d_z);

  % compute stress field from current dislocation in Fourier space
  k_x = 0:Nx-1;
  k_x(ceil((Nx+1)/2):end) = k_x(ceil((Nx+1)/2):end) - Nx; 
  k_x = 2*pi*k_x/Lx;
  k_y = 0:Ny-1;
  k_y(ceil((Ny+1)/2):end) = k_y(ceil((Ny+1)/2):end) - Ny; 
  k_y = 2*pi*k_y/Ly;
  [K_y,K_x] = meshgrid(k_y, k_x);
  norm_K_sq = K_x.^2 + K_y.^2;
  sigma_xx_fft = -2*G*i/(1-poisson_ratio)*K_y.^2./norm_K_sq.^2 ...
               .*(K_x*b(2) - K_y*b(1)).*d_z_fft;
  sigma_yy_fft = -2*G*i/(1-poisson_ratio)*K_x.^2./norm_K_sq.^2 ...
               .*(K_x*b(2) - K_y*b(1)).*d_z_fft;
  sigma_xy_fft = 2*G*i/(1-poisson_ratio)*K_x.*K_y./norm_K_sq.^2 ...
               .*(K_x*b(2) - K_y*b(1)).*d_z_fft;
  
  % correct DC component of stress fields in Fourier space
  sigma_xx_fft(1,1) = 0;
  sigma_yy_fft(1,1) = 0;
  sigma_xy_fft(1,1) = 0;

  % compute stress contribution in real space
  sigma_xx_cur = ifft2(sigma_xx_fft);
  sigma_yy_cur = ifft2(sigma_yy_fft);
  sigma_xy_cur = ifft2(sigma_xy_fft);

  % check that imaginary part is small
  if ( max(max(abs(imag(sigma_xx_cur)))) > 1e-10 ...
     | max(max(abs(imag(sigma_yy_cur)))) > 1e-10 ...
     | max(max(abs(imag(sigma_xy_cur)))) > 1e-10 )
     max(max(abs(imag(sigma_xx_cur))))
     max(max(abs(imag(sigma_yy_cur))))
     max(max(abs(imag(sigma_xy_cur))))
    warning('Imaginary component of stress nonzero!!');
  end

  sigma_xx = sigma_xx + real(sigma_xx_cur);
  sigma_yy = sigma_yy + real(sigma_yy_cur);
  sigma_xy = sigma_xy + real(sigma_xy_cur);

end

