%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% dislocation_array_stress_field_Cai
%
% Description: MATLAB function for calculating the stress field of edge
% dislocation array when each dislocation core is distributed over a 
% finite region.
%
% Input arguments:
%   b              = a row vector for the burgers vector b of the array 
%                    of parallel edge dislocations
%   r              = a position vector for the dislocation in the simulation 
%                    box
%   D              = the spacing between dislocation lines in the infinite 
%                    dislocation array which along the y-axis
%   X, Y           = computational grid (as generated by meshgrid()) on 
%                    which to calculate the stress 
%   G              = shear modulus
%   poisson_ratio  = poisson ratio
%   core_radius    = core radius
%
% Returns values:
%   sigma_xx = xx-component of stress 
%   sigma_yy = yy-component of stress 
%   sigma_xy = xy-component of stress 
%
% Function dislocation_array_stress_field_Cai assumes that the Burgers
% vector lies in the x-y plane with the tangent vector in the z-direction.
% The dislocation array is periodic and inifinite along the y-axis.
% The stress field is the infinite sum of the stress field due to a single
% non-singular edge dislocation with the specified core radius (Cai 2006). 
%
% Notes:
% (a) Reference: W. Cai et. al., J. Mech. Phys. Solids, 54: 561-587, 2006.
%
% File: dislocation_array_stress_field_Cai.m
% Written by: Adele T. Lim, MAE, Princeton University
% Date modified: 29 Aug 2007
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% CHANGE LOG
% - 2007/08/30: Kevin T. Chu
%   - Fixed bug in sigma_yy.
%   - Modified to take physical parameters are arguments.
%   - Cleaned up code and documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [sigma_xx sigma_yy sigma_xy] = ...
  dislocation_array_stress_field_Cai (b, r, D, X, Y, ...
                                      G, poisson_ratio, ...
                                      core_radius)

A = core_radius/D;
X = ( X - r(1) )/D;
Y = ( Y - r(2) )/D;

X_a = sqrt(X.^2 + A^2);

sigma_0_D = 0.5*G/(1-poisson_ratio)/D./( cosh(2*pi*X_a) - cos(2*pi*Y) ).^2;

sigma_xx = -sigma_0_D.*(...
    b(1)*sin(2*pi*Y).*( cosh(2*pi*X_a) - cos(2*pi*Y) ...
                      + 2*pi*X_a.*sinh(2*pi*X_a) ) ...
   -b(2)*2*pi*X.*( cosh(2*pi*X_a).*cos(2*pi*Y) - 1 ) ...
    );

sigma_xy =  sigma_0_D.*(...
    b(1)*2*pi*X.*( cosh(2*pi*X_a).*cos(2*pi*Y) - 1 ) ...
   -b(2)*( sin(2*pi*Y).*( cosh(2*pi*X_a) - cos(2*pi*Y) ...
                         - 2*pi*X_a.*sinh(2*pi*X_a) ) ...
         + 2*pi*A^2./X_a.*sinh(2*pi*X_a).*sin(2*pi*Y) ) ...
    );

sigma_yy = -sigma_0_D.*(...
    b(1)*sin(2*pi*Y).*( cosh(2*pi*X_a) - cos(2*pi*Y) ...
                      - 2*pi*X_a.*sinh(2*pi*X_a) ...
                      + 2*A^2*pi./X_a.*sinh(2*pi*X_a) ) ...
    -b(2)*( X./X_a.*sinh(2*pi*X_a).*( cosh(2*pi*X_a) - cos(2*pi*Y) ) ...
             .*( 2 + A^2./X_a.^2 ) ...
          - 2*pi*X.*( cosh(2*pi*X_a).*cos(2*pi*Y) - 1 ) ...
            .*(1 - A^2./X_a.^2 ) )...
    );