%  coeffhs
%  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) in an elastic half space.

	global	Uxn Uxs Uyn Uys Sxyn Sxys Syyn Syys Sxxn Sxxs

% Define trig functions and coordinates pertaining to the element
% and create arrays of the same size as the observation array
unit = ones(size(x));
COSB2e = (COSBe.*COSBe).*unit;
SINB2e = (SINBe.*SINBe).*unit;
COS2Be = (COSBe.*COSBe-SINBe.*SINBe).*unit;
SIN2Be = (2.*SINBe.*COSBe).*unit;
COS3Be = (4*COSBe.^3 - 3*COSBe).*unit;
SIN3Be = (3*SINBe - 4*SINBe.^3).*unit;
COS4Be = (8*COSBe.^4 - 8*COSB2e + 1).*unit;
SIN4Be = (4*SINBe.*COSBe - 8*COSBe.*SINBe.^3).*unit;
xe = xe.*unit;
ye = ye.*unit;

% set tolerance to correct roundoff errors of trig functions
bad = find(abs(COS2Be)<1e-10);
COS2Be(bad) = 0;
bad = find(abs(COS3Be)<1e-10);
COS3Be(bad) = 0;
bad = find(abs(SIN3Be)<1e-10);
SIN3Be(bad) = 0;
bad = find(abs(SIN4Be)<1e-10);
SIN4Be(bad) = 0;
clear bad

% 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;		%equation 7.4.2
Yb = -(x-xe).*SINBe + (y-ye).*COSBe;

% Coordinates of the image dislocation
Xbi = (x-xe).*COSBe - (y+ye).*SINBe;		%equation 7.4.6
Ybi = (x-xe).*SINBe + (y+ye).*COSBe;

% Fix roundoff errors in Ybi and Yb from trig function problems
bad = find(abs(Ybi)<1e-10);
bad2 = find(abs(Yb)<1e-10);
Ybi(bad) = 0;
Yb(bad2) = 0;
clear bad bad2

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

% Same thing for the image dislocation
R1Si = (Xbi-a).^2 + Ybi.^2;
R2Si = (Xbi+a).^2 + Ybi.^2;

%  Calculate intermediate functions Fna for the actual dislocation
F2a = CON.*( log(sqrt(R1S)) - log(sqrt(R2S)) );
F3a = evalf3(a,Xb,Yb,CON);
F4a = CON.*(Yb./R1S - Yb./R2S);
F5a = CON.*( (Xb-a)./R1S - (Xb+a)./R2S );
F6a = CON.*( ((Xb-a).^2 - Yb.^2)./R1S.^2 - ((Xb+a).^2 - Yb.^2)./R2S.^2 );
F7a = 2.*CON.*Yb.*( (Xb-a)./R1S.^2 - (Xb+a)./R2S.^2 );

% Calculate intermediate functions Fnifor the image dislocation
F2i = CON.*( log(sqrt(R1Si)) - log(sqrt(R2Si)) );
F3i = evalf3(a,Xbi,Ybi,CON);
F4i = CON.*(Ybi./R1Si - Ybi./R2Si);
F5i = CON.*( (Xbi-a)./R1Si - (Xbi+a)./R2Si );
F6i = CON.*( ((Xbi-a).^2 - Ybi.^2)./R1Si.^2 - ((Xbi+a).^2 - Ybi.^2)./R2Si.^2 );
F7i = 2.*CON.*Ybi.*( (Xbi-a)./R1Si.^2 - (Xbi+a)./R2Si.^2 );
F8i = 2.*CON.*Ybi.*( (Ybi.^2 - 3.*(Xbi-a).^2)./R1Si.^3 -.....
      (Ybi.^2 - 3.*(Xbi+a).^2)./R2Si.^3 );
F9i = 2.*CON.*( ((Xbi-a).*((Xbi-a).^2 - 3.*Ybi.^2))./R1Si.^3 -.....
      ((Xbi+a).*((Xbi+a).^2 - 3.*Ybi.^2))./R2Si.^3 );

%  See equations 7.4.3 and 7.4.4 in Crouch and Starfield
%  Calculate ACTUAL DISPLACEMENT components due to unit SHEAR disp. disc.
UxsA = -PR1.*SINBe.*F2a + PR2.*COSBe.*F3a + Yb.*(SINBe.*F4a - COSBe.*F5a);
UysA =  PR1.*COSBe.*F2a + PR2.*SINBe.*F3a - Yb.*(COSBe.*F4a + SINBe.*F5a);

%  Calculate ACTUAL DISPLACEMENT components due to unit NORMAL disp. disc.
UxnA = -PR1.*COSBe.*F2a - PR2.*SINBe.*F3a - Yb.*(COSBe.*F4a + SINBe.*F5a);
UynA = -PR1.*SINBe.*F2a + PR2.*COSBe.*F3a - Yb.*(SINBe.*F4a - COSBe.*F5a);

%  Calculate ACTUAL STRESS components due to unit SHEAR disp. discontinuity
SxxsA = CONS.*(2.*COSB2e.*F4a + SIN2Be.*F5a + Yb.*(COS2Be.*F6a - SIN2Be.*F7a));
SyysA = CONS.*(2.*SINB2e.*F4a - SIN2Be.*F5a - Yb.*(COS2Be.*F6a - SIN2Be.*F7a));
SxysA = CONS.*(SIN2Be.*F4a - COS2Be.*F5a + Yb.*(SIN2Be.*F6a + COS2Be.*F7a));

%  Calculate ACTUAL STRESS components due to unit NORMAL displacement disc.
SxxnA = CONS.*(-F5a + Yb.*(SIN2Be.*F6a + COS2Be.*F7a));
SyynA = CONS.*(-F5a - Yb.*(SIN2Be.*F6a + COS2Be.*F7a));
SxynA = CONS.*(-Yb.*(COS2Be.*F6a - SIN2Be.*F7a));


%  See equations 7.4.8 and 7.4.9 in Crouch and Starfield
%  Calculate IMAGE AND SUPPLEMENTAL DISPLACEMENT components due to unit SHEAR 
%  displacement discontinuity
UxsIS = PR1.*SINBe.*F2i - PR2.*COSBe.*F3i +....
        (PR3.*(y.*SIN2Be - Yb.*SINBe) + 2.*y.*SIN2Be).*F4i +...
        (PR3.*(y.*COS2Be - Yb.*COSBe) - y.*(1-2.*COS2Be)).*F5i +...
        2.*y.*(y.*SIN3Be - Yb.*SIN2Be).*F6i -....
	2.*y.*(y.*COS3Be - Yb.*COS2Be).*F7i;

UysIS = -PR1.*COSBe.*F2i - PR2.*SINBe.*F3i -...
        (PR3.*(y.*COS2Be - Yb.*COSBe) + y.*(1-2.*COS2Be)).*F4i +...
        (PR3.*(y.*SIN2Be - Yb.*SINBe) - 2.*y.*SIN2Be).*F5i +...
        2.*y.*(y.*COS3Be - Yb.*COS2Be).*F6i +......
	2.*y.*(y.*SIN3Be - Yb.*SIN2Be).*F7i;


%  Calculate IMAGE AND SUPPLEMENTAL DISPLACEMENT components due to unit NORMAL 
%  displacement discontinuity
UxnIS = PR1.*COSBe.*F2i + PR2.*SINBe.*F3i -.....
	(PR3.*(y.*COS2Be - Yb.*COSBe) - y).*F4i +....
	PR3.*(y.*SIN2Be - Yb.*SINBe).*F5i -....
	2.*y.*(y.*COS3Be - Yb.*COS2Be).*F6i -....
	2.*y.*(y.*SIN3Be - Yb.*SIN2Be).*F7i;

UynIS = PR1.*SINBe.*F2i - PR2.*COSBe.*F3i -..... 
	PR3.*(y.*SIN2Be - Yb.*SINBe).*F4i -.....
 	(PR3.*(y.*COS2Be - Yb.*COSBe) + y).*F5i +.....
 	2.*y.*(y.*SIN3Be - Yb.*SIN2Be).*F6i -....
 	2.*y.*(y.*COS3Be - Yb.*COS2Be).*F7i;


%  Calculate IMAGE AND SUPPLEMENTAL STRESS components due to unit SHEAR 
%  displacement discontinuity
SxxsIS = CONS.*(F4i - 3.*(COS2Be.*F4i - SIN2Be.*F5i) +....
	(2.*y.*(COSBe - 3.*COS3Be) + 3.*Yb.*COS2Be).*F6i +.....
	(2.*y.*(SINBe - 3.*SIN3Be) + 3.*Yb.*SIN2Be).*F7i -.....
	2.*y.*(y.*COS4Be - Yb.*COS3Be).*F8i -....
	2.*y.*(y.*SIN4Be - Yb.*SIN3Be).*F9i);

SyysIS = CONS.*(F4i - (COS2Be.*F4i - SIN2Be.*F5i) -....
	(4.*y.*SINBe.*SIN2Be - Yb.*COS2Be).*F6i +.....
 	(4.*y.*SINBe.*COS2Be + Yb.*SIN2Be).*F7i +....
 	2.*y.*(y.*COS4Be - Yb.*COS3Be).*F8i +....
 	2.*y.*(y.*SIN4Be - Yb.*SIN3Be).*F9i);

SxysIS = CONS.*(SIN2Be.*F4i + COS2Be.*F5i +.... 
	(2.*y.*SINBe.*(1+4.*COS2Be) - Yb.*SIN2Be).*F6i +.....
	(2.*y.*COSBe.*(3-4.*COS2Be) + Yb.*COS2Be).*F7i +..... 
	2.*y.*(y.*SIN4Be - Yb.*SIN3Be).*F8i -.....
	2.*y.*(y.*COS4Be - Yb.*COS3Be).*F9i);


%  Calculate IMAGE AND SUPPLEMENTAL STRESS components due to unit NORMAL
%  displacement discontinuity
SxxnIS = CONS.*(F5i + (2.*y.*(SINBe - 2.*SIN3Be) +.... 
	3.*Yb.*SIN2Be).*F6i - (2.*y.*(COSBe - 2.*COS3Be) +....
	3.*Yb.*COS2Be).*F7i - 2.*y.*(y.*SIN4Be - Yb.*SIN3Be).*F8i +....
	2.*y.*(y.*COS4Be - Yb.*COS3Be).*F9i);

SyynIS = CONS.*(F5i - (2.*y.*SINBe - Yb.*SIN2Be).*F6i +.... 
	(2.*y.*COSBe - Yb.*COS2Be).*F7i +....
	2.*y.*(y.*SIN4Be - Yb.*SIN3Be).*F8i -.... 
	2.*y.*(y.*COS4Be - Yb.*COS3Be).*F9i);

SxynIS = CONS.*((4.*y.*SINBe.*SIN2Be + Yb.*COS2Be).*F6i -.... 
	(4.*y.*SINBe.*COS2Be - Yb.*SIN2Be).*F7i -.....
	2.*y.*(y.*COS4Be - Yb.*COS3Be).*F8i -.... 
	2.*y.*(y.*SIN4Be - Yb.*SIN3Be).*F9i);


%Add actual, image, and supplement solutions to get final solution
Uxs = UxsA+UxsIS;
Uys = UysA+UysIS;
Uxn = UxnA+UxnIS;
Uyn = UynA+UynIS;

Sxxs = SxxsA+SxxsIS;
Syys = SyysA+SyysIS;
Sxys = SxysA+SxysIS;
Sxxn = SxxnA+SxxnIS;
Syyn = SyynA+SyynIS;
Sxyn = SxynA+SxynIS;

% 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;

% Clear unneeded parameters
clear F2a F3a F4a F5a F6a F7a F2i F3i F4i F5i F6i F7i F8i F9i
clear unit xe ye Xb Yb Xbi Ybi R1S R2S R1Si R2Si 
clear COS2Be COS3Be COS4Be SIN2Be SIN3Be SIN4Be 
clear SxxnA SyynA SxynA SxxsA SyysA SxysA UxnA UynA UxsA UysA
clear SxxnIS SyynIS SxynIS SxxsIS SyysIS SxysIS UxnIS UynIS UxsIS UysIS