%  This is the MATLAB version of SUBROUTINE COEFF from TWODD
%  It calculates the xy stress and displacement components
%  at all points (x,y) in an observation array
%  due to unit normal and unit shear displacement 
%  discontinuities at element j(xe,ye).

% Define trig functions and coordinates pertaining to the element
% and create arrays of the same size as the observation array
unit = ones(size(x));
COS2Be = (COSBe.*COSBe-SINBe.*SINBe).*unit;
SIN2Be = (2.*SINBe.*COSBe).*unit;
COSB2e = (COSBe.*COSBe).*unit;
SINB2e = (SINBe.*SINBe).*unit;
xe = xe.*unit;
ye = ye.*unit;	

% Calculate coordinates of observation points in reference frame
% centered at the element and oriented along the element.
% The transformation requires a translation and a rotation
Xb = (x-xe).*COSBe + (y-ye).*SINBe;
Yb = -(x-xe).*SINBe + (y-ye).*COSBe;
% Flag all observation points for which Yb is not 0
i1 = find(Yb);
% Flag any observation points on the element
i2 = find(Yb == 0 & abs(Xb) < a);

% Find squares of the distances from obs. pts. to element ends
R1S = (Xb-a).*(Xb-a) + Yb.*Yb;
R2S = (Xb+a).*(Xb+a) + Yb.*Yb;

%  Calculate intermediate functions F
FL1 = 0.5.*log(R1S);
FL2 = 0.5.*log(R2S);
FB2 = CON.*(FL1-FL2);
% FB3 = 0 for pts colinear with element, CON*pi for pts. on element 
% FB3 = difference of arc tangents for all other pts.
FB3 = zeros(size(x));				
FB3(i1) =  -1.*CON.*(atan((Xb(i1)+a)./Yb(i1)) - atan((Xb(i1)-a)./Yb(i1)));
FB3(i2) = CON.*pi.*ones(size(i2));
FB4 = CON.*(Yb./R1S - Yb./R2S);
FB5 = CON.*((Xb-a)./R1S - (Xb+a)./R2S);
FB6 = CON.*(((Xb-a).^2 -Yb.*Yb)./R1S.^2 - ((Xb+a).^2-Yb.*Yb)./R2S.^2);
FB7 = 2.*CON.*Yb.*((Xb-a)./R1S.^2 - (Xb+a)./R2S.^2);

%  Calculate DISPLACEMENT components due to unit SHEAR disp. discontinuity
Uxs = -PR1.*SINBe.*FB2 + PR2.*COSBe.*FB3 + Yb.*(SINBe.*FB4-COSBe.*FB5);
Uys =  PR1.*COSBe.*FB2 + PR2.*SINBe.*FB3 - Yb.*(COSBe.*FB4+SINBe.*FB5);

%  Calculate DISPLACEMENT components due to unit NORMAL disp. discontinuity
Uxn = -PR1.*COSBe.*FB2 - PR2.*SINBe.*FB3 - Yb.*(COSBe.*FB4+SINBe.*FB5);
Uyn = -PR1.*SINBe.*FB2 + PR2.*COSBe.*FB3 - Yb.*(SINBe.*FB4-COSBe.*FB5);

%  Calculate STRESS components due to unit SHEAR disp. discontinuity
Sxxs = CONS.*(2.*COSB2e.*FB4 + SIN2Be.*FB5 + Yb.*(COS2Be.*FB6-SIN2Be.*FB7));
Syys = CONS.*(2.*SINB2e.*FB4 - SIN2Be.*FB5 - Yb.*(COS2Be.*FB6-SIN2Be.*FB7));
Sxys = CONS.*(SIN2Be.*FB4 - COS2Be.*FB5 + Yb.*(SIN2Be.*FB6 + COS2Be.*FB7));

%  Calculate STRESS components due to unit NORMAL displacement discontinuity
Sxxn = CONS.*(-FB5 + Yb.*(SIN2Be.*FB6 + COS2Be.*FB7));
Syyn = CONS.*(-FB5 - Yb.*(SIN2Be.*FB6 + COS2Be.*FB7));
Sxyn = CONS.*(-Yb.*(COS2Be.*FB6 - SIN2Be.*FB7));

% Multiply the components by the shear and normal displacement 
% discontintuities.  For the first pass through coeff, Sd=Nd=1. 
Uxs =  Sd*Uxs;
Uys =  Sd*Uys;
Uxn =  Nd*Uxn;
Uyn =  Nd*Uyn;
Sxxs = Sd*Sxxs;
Syys = Sd*Syys;
Sxys = Sd*Sxys;
Sxxn = Nd*Sxxn;
Syyn = Nd*Syyn;
Sxyn = Nd*Sxyn;