function [C,B,Ds,Dn,UxN,UyN,UsN,UnN,UsP,UnP,SIGxx,SIGyy,SIGxy,...
Ux,Uy,Sxx,Syy,Sxy,S1,S2,tau,mean,theta] = twoddhs_func(e_dat,r_dat,b_dat,X,Y)
%  Function version of two-dimensional half-space boundary element program
%  All the functions twodd_func calls are internally stored subfunctions
%  with stress boundary conditions for Matlab 4.1
%  This program is modified from Appendix B of the book
%  "Boundary element methods in Solid Mechanics"
%  by S.L. Crouch and A.M. Starfield, 1983
% Last revised on 6/29/03.  By Stephen J. Martel.
% Examples:
%load e_dat
%load r_dat
%load b_dat
%xmin = -1.9;
%xmax = +1.9;
%dx = 0.2;
%ymin = -2.8;
%ymax = 0;
%dy = 0.2;
%xgrid = xmin:dx:xmax;
%ygrid = ymin:dy:ymax;
%[X,Y] = meshgrid(xgrid,ygrid);
%[C,B,Ds,Dn,UxN,UyN,UsN,UnN,UsP,UnP,SIGxx,SIGyy,SIGxy] = twoddhs_func(e_dat,r_dat,b_dat)
%[C,B,Ds,Dn,UxN,UyN,UsN,UnN,UsP,UnP,SIGxx,SIGyy,SIGxy,...
%Ux,Uy,Sxx,Syy,Sxy,S1,S2,tau,mean,theta] = twoddhs_func(e_dat,r_dat,b_dat,X,Y)


%  Read elastic constants, remote stress field data, and boundary data from 
%  data files elastic.dat, remote.dat, and boundary.dat, respectively

% Elastic constants.  PR and E should be in one row;
	PR = e_dat(1);			%  PR = Poisson's ratio;
	E = e_dat(2);			%  E = Young's modulus;

% Ambient field.  PXX, PYY, PXY should be in one row;
	Pxx = r_dat(1);		%  Remote stress sigma xx;
	Pyy = r_dat(2);		%  Remote stress sigma yy;
	Pxy = r_dat(3);		%  Remote stress sigma xy;

% Boundary element geometry and boundary conditions
%  Data in each row of this file pertains to one element
%  Now separate out data types by column to form column vectors
	XBEG = b_dat(:,1);	%  element endpoint coordinate
	YBEG = b_dat(:,2);	%  element endpoint coordinate
	XEND = b_dat(:,3);	%  element endpoint coordinate
	YEND = b_dat(:,4);	%  element endpoint coordinate
	BVs = b_dat(:,5);	%  Shear traction boundary condition
	BVn = b_dat(:,6);	%  Normal traction boundary condition

	NUM = length(BVs);		% NUM = number of elements

%  Define new elastic constants to be used internally
	CON = 1/(4*pi*(1-PR));
	CONS = E/(1+PR);
	PR1 = 1-2*PR;
	PR2 = 2*(1-PR);
	PR3 = 1-PR;

%  Define locations, lengths, orientations, and boundary conditions for boundary elements
%  This information is stored here as column vectors of size (NUMx1)
	x = (XBEG + XEND)/2;				% element midpoints
	y = (YBEG + YEND)/2;				% element midpoints
	XE = x;
	YE = y;
	XD = XEND-XBEG;
	YD = YEND-YBEG;
	A = sqrt(XD.*XD +YD.*YD)/2;			% half-length of element
%  In the lines below, B = beta = orientation of element relative to global x-axis
	SINB = YD./(2*A);					% sine beta
	COSB = XD./(2*A);					% cosine beta
	SIN2B = 2*SINB.*COSB;				% sine 2*beta
	COS2B = COSB.*COSB-SINB.*SINB;		% cosine 2*beta
	SINB2 = SINB.*SINB;					% sine squared of beta
	COSB2 = COSB.*COSB;					% cosine squared of beta

%  Adjust stress boundary conditions to account for remote stresses.
%  First resolve components of the remote stress onto the plane of each element...
	SIGs = (Pyy-Pxx).*SIN2B/2 + Pxy.*COS2B;
	SIGn = Pxx.*SINB2 - Pxy.*SIN2B + Pyy.*COSB2;
%  and then subtract the resolved remote stress from the boundary stresses.
	BVs = BVs - SIGs;
	BVn = BVn - SIGn;
%  The remote stress will be added back at the end of the solution.

%  Compute influence coefficients between boundary elements.
%  C(i,j) = effect at obs pt i due to a load at element j.  
%  The organization of the coefficients here is slightly different 
%  from Crouch & Starfield
%	First dimension the vectors storing the influence coefficients
	nobs=NUM*NUM;
	[Uxsv,Uysv,Uxnv,Uynv,Sxxsv,Syysv,Sxysv,Sxxnv,Syynv,Sxynv] = set1_func(nobs);			
	for i=1:NUM
		first = (i-1)*NUM+1;
		last = i*NUM;
		xe = XE(i);
		ye = YE(i);
		COSBe = COSB(i);
		SINBe = SINB(i);
		a = A(i);
		Sd = 1;
		Nd = 1;
		[Uxs,Uys,Uxn,Uyn,Sxxs,Syys,Sxys,Sxxn,Syyn,Sxyn] = ...
		coeffhs_func(x,y,xe,ye,a,COSBe,SINBe,CON,CONS,PR1,PR2,PR3,Sd,Nd);
		% Append sequentially produced vectors to master column vectors 
		[Uxsv,Uysv,Uxnv,Uynv,Sxxsv,Syysv,Sxysv,Sxxnv,Syynv,Sxynv] = ...
		append1_func(first,last,Uxs,Uys,Uxn,Uyn,Sxxs,Syys,Sxys,Sxxn,Syyn,Sxyn,...
		Uxsv,Uysv,Uxnv,Uynv,Sxxsv,Syysv,Sxysv,Sxxnv,Syynv,Sxynv);
	end

% Calculate trig functions
% and direction cosine matrices for stress transformations
	[SINNB,COSSB,SINN2B,COSS2B,SINNB2,COSSB2] = ...
	dircos1_func(SINB,COSB,SIN2B,COS2B,SINB2,COSB2,NUM);

%  Reshape column vectors into square matrices, then clear column vectors
	dimx = NUM;
	dimy = NUM;
	[Uxs,Uxn,Uys,Uyn,Sxxs,Sxxn,Syys,Syyn,Sxys,Sxyn] = ...
	square1_func(Uxsv,Uxnv,Uysv,Uynv,Sxxsv,Sxxnv,Syysv,Syynv,Sxysv,Sxynv,dimx,dimy);

%  Convert xy stresses returned from coeff to normal and shear tractions on elements
	Cisjs = (Syys-Sxxs).*(SINN2B)/2 + Sxys.*COSS2B;
	Cisjn = (Syyn-Sxxn).*(SINN2B)/2 + Sxyn.*COSS2B;
	Cinjs = Sxxs.*SINNB2 - Sxys.*SINN2B + Syys.*COSSB2;
	Cinjn = Sxxn.*SINNB2 - Sxyn.*SINN2B + Syyn.*COSSB2;

%  Solve system of algebraic equations to get the displacement discontinuities
%  First regroup the influence coefficient and boundary condition submatrices as shown:
%  |   [Cisjs]	[Cisjn]	|	| [Ds]	|	=	| [Bs]	|
%  |   [Cinjs]	[Cinjn] |	| [Dn]	|		| [Bn]	|
	C = [Cisjs, Cisjn; Cinjs, Cinjn];		% C = influence coefficients;
%	clear Cisjs Cisjn Cinjs Cinjn;
	B = [BVs; BVn];							% B= Stress boundary conditions;
%  Solve the matrix system [C][D]=[B] to get the displacement discontinuity vector D
	D= C\B;
%  Separate the D vector into subvectors Ds and Dn
	Ds = D(1:NUM);				% Ds elements are in upper half of D
	Dn = D(NUM+1:2*NUM);		% Dn elements are in lower half of D
%	clear B C D

%  Compute displacements and stresses on boundary elements
%  The stresses obtained should match the boundary conditions

%  First find the stresses and displacements, in an x-y reference frame, by superposition
	UxN = Uxs*Ds + Uxn*Dn;			% x-displacement on negative side of i;
	UyN = Uys*Ds + Uyn*Dn;			% y-displacement on negative side of i;
	SIGxx = Pxx + Sxxs*Ds + Sxxn*Dn;	% Note that remote stress is added back in;
	SIGyy = Pyy + Syys*Ds + Syyn*Dn;	% Note that remote stress is added back in;
	SIGxy = Pxy + Sxys*Ds + Sxyn*Dn;	% Note that remote stress is added back in;

%  Now resolve DISPLACEMENTS in the xy frame into shear and normal components
%  Start by solving for displacements on the NEGATIVE (N) sides of the elements...;
	UsN = UxN.*COSB + UyN.*SINB;
	UnN = -UxN.*SINB + UyN.*COSB;

%  and then add the displacement discontinuity at each element to get 
%  the displacement components on the positive (P) sides of the elements
	UsP = UsN - Ds;
	UnP = UnN - Dn;

%  The last step is to resolve the xy STRESSES to normal and shear stresses on the elements.
	SIGs = (SIGyy-SIGxx).*(SIN2B/2)  +  SIGxy.*COS2B;
	SIGn = SIGxx.*SINB2  -  SIGxy.*SIN2B  +  SIGyy.*COSB2;

%  COMPUTATION OF DISPLACEMENTS AND STRESSES AT SPECIFIED OBSERVATION POINTS I IN BODY
%  Proceed if output on observation grid X,Y is desired, otherwise end
if nargout > 13
	[dimx,dimy] = size(X);
	NUM2 = dimx*dimy;
%  Calculate influence coefficients at gridpoints
	x = X(:);			% Grid nodes in column vector form
	y = Y(:);			% Grid nodes in column vector form
%	First dimension the vectors storing the influence coefficients
	nobs=NUM2*NUM;
	[Uxsv,Uysv,Uxnv,Uynv,Sxxsv,Syysv,Sxysv,Sxxnv,Syynv,Sxynv] = set1_func(nobs);
% Now loop through grid element by element
% The end result will be ten column vectors each having NUM2*NUM rows			
	for i=1:NUM
		first = (i-1)*NUM2+1;
		last = i*NUM2;
		xe = XE(i);
		ye = YE(i);
		COSBe = COSB(i);
		SINBe = SINB(i);
		a = A(i);
		Sd = Ds(i);
		Nd = Dn(i);
		[Uxs,Uys,Uxn,Uyn,Sxxs,Syys,Sxys,Sxxn,Syyn,Sxyn] = ...
		coeffhs_func(x,y,xe,ye,a,COSBe,SINBe,CON,CONS,PR1,PR2,PR3,Sd,Nd);
		% Append sequentially produced vectors to master column vectors 
		[Uxsv,Uysv,Uxnv,Uynv,Sxxsv,Syysv,Sxysv,Sxxnv,Syynv,Sxynv] = ...
		append1_func(first,last,Uxs,Uys,Uxn,Uyn,Sxxs,Syys,Sxys,Sxxn,Syyn,Sxyn,...
		Uxsv,Uysv,Uxnv,Uynv,Sxxsv,Syysv,Sxysv,Sxxnv,Syynv,Sxynv);
	end

% Reshape column vectors into square matrices the size of the
% observation grid, and in the process sum the contributions
% from each element 
	[Ux,Uy,Sxx,Syy,Sxy] = ...
square2_func(Uxsv,Uxnv,Uysv,Uynv,Sxxsv,Sxxnv,Syysv,Syynv,Sxysv,Sxynv,NUM2,NUM,dimx,dimy);

% Add back the remote field
	Sxx = Pxx + Sxx;
	Syy = Pyy + Syy;
	Sxy = Pxy + Sxy;

	
	S1 = sig1(Sxx,Syy,Sxy);
	S2 = sig2(Sxx,Syy,Sxy);
	tau = taumax(Sxx,Syy,Sxy);
	mean = ave(Sxx,Syy,Sxy);
	theta = angp(Sxx,Syy,Sxy); 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Alphabetical listing of local subfunctions called in the primary function twodd_func 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function theta=angp(SIGxx,SIGyy,SIGxy)
% function angp.  Calculates a principal stress orientation
theta = 0.5*atan2(SIGxy,(SIGxx-SIGyy)/2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Uxsv,Uysv,Uxnv,Uynv,Sxxsv,Syysv,Sxysv,Sxxnv,Syynv,Sxynv] = ...
append1_func(first,last,Uxs,Uys,Uxn,Uyn,Sxxs,Syys,Sxys,Sxxn,Syyn,Sxyn,...
Uxsv,Uysv,Uxnv,Uynv,Sxxsv,Syysv,Sxysv,Sxxnv,Syynv,Sxynv)
% Function to save influence coeffients from twodd.m in column vectors
% by appending new column vectors generated by the first set of
% calls to coeff to the previously generated column vectors
Uxsv (first:last) = Uxs;
Uysv (first:last) = Uys;
Uxnv (first:last) = Uxn;
Uynv (first:last) = Uyn;
Sxxsv (first:last) = Sxxs;
Syysv (first:last) = Syys;
Sxysv (first:last) = Sxys;
Sxxnv (first:last) = Sxxn;
Syynv (first:last) = Syyn;
Sxynv (first:last) = Sxyn;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mean=ave(SIGxx,SIGyy,SIGxy)
mean = (SIGxx+SIGyy)/2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Uxs,Uys,Uxn,Uyn,Sxxs,Syys,Sxys,Sxxn,Syyn,Sxyn] = ...
coeffhs_func(x,y,xe,ye,a,COSBe,SINBe,CON,CONS,PR1,PR2,PR3,Sd,Nd)
%  This is the MATLAB half-space function 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.

% 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;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [SINNB,COSSB,SINN2B,COSS2B,SINNB2,COSSB2] = ...
dircos1_func(SINB,COSB,SIN2B,COS2B,SINB2,COSB2,NUM)
% function to calculate direction cosine matrices
% for stress transformations
%axpx = cos(x'x)			= cos(thetap - theta)
%axpy = cos(x'y) = sin(x'x)	= sin(thetap - theta)
%aypx = cos(y'x) = -sin(x'x)
%aypy = cos(y'y) = cos(x'x)
% cos(A-B) = cosA*cosB + sinA*sinB
% sin(A-B) = sinA*cosB - cosA*sinB

%axpx = COSB*(COSB') + SINB*(SINB');
%axpy = COSB*(COSB') + SINB*(SINB');
%aypx = -aypx;
%aypy = axpx;

SINNB = SINB*(ones(1,NUM));
COSSB = COSB*(ones(1,NUM));
SINN2B = SIN2B*(ones(1,NUM));
COSS2B = COS2B*(ones(1,NUM));
SINNB2 = SINB2*(ones(1,NUM));
COSSB2 = COSB2*(ones(1,NUM));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function F3 = evalF3(A,X,Y,CON)
%	function that evaluates arctangent terms in coeffHS
%	checks for observation points that are colinear with
%	actual and image element endpoints as described in
%	Crouch and Starfield (pgs 50-51)
%	i finds all nonzero values of Y
%	j finds all points 'within' the element
%	k finds all points in the plane of, but outside of each element
%	note for index j the solution is +pi 

i = find(Y);
j = find( Y==0 & abs(X)<A );
k = find( Y==0 & abs(X)>A );
F3(i) = atan2(Y(i),(X(i)-A)) - atan2(Y(i),(X(i)+A));
F3(j) = -pi*ones(size(j));
F3(k) = zeros(size(k));

F3 = -CON.*(F3)';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Uxsv,Uysv,Uxnv,Uynv,Sxxsv,Syysv,Sxysv,Sxxnv,Syynv,Sxynv] = set1_func(nobs)
% Function to dimension the vectors storing the influence coefficients
% from the first set of calls to coeff
Uxsv   = zeros(nobs,1);
Uysv   = Uxsv;
Uxnv   = Uxsv;
Uynv   = Uxsv;
Sxxsv  = Uxsv;
Syysv  = Uxsv;
Sxysv  = Uxsv;
Sxxnv  = Uxsv;
Syynv  = Uxsv;
Sxynv  = Uxsv;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function S1=sig1(Sxx,Syy,Sxy)
% function S1=sig1(SIGxx,SIGyy,SIGxy)
% Calculates the most tensile 2-D principal stress magnitude
% Input parameters
% Sxx = sigma xx
% Syy = sigma yy
% Sxy = sigma xy
% Output parameter
% S1 = most tensile principal stress
% Example 
% S1=sig1(4,-4,3)

S1 = (Sxx+Syy)/2 + sqrt( ((Sxx-Syy)/2).^2 + Sxy.^2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function S2=sig2(Sxx,Syy,Sxy)
% function S1=sig2(Sxx,Syy,Sxy)
% Calculates the most tensile 2-D principal stress magnitude
% Input parameters
% Sxx = sigma xx
% Syy = sigma yy
% Sxy = sigma xy
% Output parameter
% S1 = most tensile principal stress
% Example 
% S1=sig2(4,-4,3)

S2 = (Sxx+Syy)/2 - sqrt( ((Sxx-Syy)/2).^2 + Sxy.^2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Uxs,Uxn,Uys,Uyn,Sxxs,Sxxn,Syys,Syyn,Sxys,Sxyn] = ...
square1_func(Uxsv,Uxnv,Uysv,Uynv,Sxxsv,Sxxnv,Syysv,Syynv,Sxysv,Sxynv,dimx,dimy)
% function to turn column vectors into square matrices
% and then clear the unneeded column vectors
Uxs = reshape(Uxsv,dimx,dimy);
Uxn = reshape(Uxnv,dimx,dimy);
Uys = reshape(Uysv,dimx,dimy);
Uyn = reshape(Uynv,dimx,dimy);
Sxxs = reshape(Sxxsv,dimx,dimy);
Sxxn = reshape(Sxxnv,dimx,dimy);
Syys = reshape(Syysv,dimx,dimy);
Syyn = reshape(Syynv,dimx,dimy);
Sxys = reshape(Sxysv,dimx,dimy);
Sxyn = reshape(Sxynv,dimx,dimy);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Ux,Uy,Sxx,Syy,Sxy] = ...
square2_func(Uxsv,Uxnv,Uysv,Uynv,Sxxsv,Sxxnv,Syysv,Syynv,Sxysv,Sxynv,NUM2,NUM,dimx,dimy)
% Function square2 to help turn column vectors into square matrices
% and clears unneeded column vectors.  

% First turn the entering column vectors into matrices, 
% where the number of rows is the number of observation points (NUM2)
% and the number of columns equals the number of elements.
% The column vectors are cleared once they are done with.
Uxs = reshape(Uxsv,NUM2,NUM);
Uxn = reshape(Uxnv,NUM2,NUM);
Uys = reshape(Uysv,NUM2,NUM);
Uyn = reshape(Uynv,NUM2,NUM);
Sxxs = reshape(Sxxsv,NUM2,NUM);
Sxxn = reshape(Sxxnv,NUM2,NUM);	
Syys = reshape(Syysv,NUM2,NUM);	
Syyn = reshape(Syynv,NUM2,NUM);	
Sxys = reshape(Sxysv,NUM2,NUM);
Sxyn = reshape(Sxynv,NUM2,NUM);

% Next, sum the entries in each row to superpose
% the contributions of all the elements.
% This step will yield column vectors  
% See p. 488 of the Matlab manual for this
Uxs = sum(Uxs')';
Uxn = sum(Uxn')';
Uys = sum(Uys')';
Uyn = sum(Uyn')';
Sxxs = sum(Sxxs')';
Sxxn = sum(Sxxn')';
Syys = sum(Syys')';
Syyn = sum(Syyn')';
Sxys = sum(Sxys')';
Sxyn = sum(Sxyn')';

% Now, reshape the column vectors into matrices
% that have the dimensions of the observation grid,
% adding the contributions from the shear and
% normal displacement discontinuities 
Ux = reshape(Uxs,dimx,dimy) + reshape(Uxn,dimx,dimy);
Uy = reshape(Uys,dimx,dimy) + reshape(Uyn,dimx,dimy);
Sxx = reshape(Sxxs,dimx,dimy) + reshape(Sxxn,dimx,dimy);
Syy = reshape(Syys,dimx,dimy) + reshape(Syyn,dimx,dimy);
Sxy = reshape(Sxys,dimx,dimy) + reshape(Sxyn,dimx,dimy);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function tau=taumax(Sxx,Syy,Sxy)
% function tau=taumax(Sxx,Syy,Sxy)
% Calculates maximum shear stress magnitude (2-D)
% Input parameters
% Sxx = sigma xx
% Syy = sigma yy
% Sxy = sigma xy
% Output parameter
% tau = maximum shear stress
% Example 
% tau=taumax(4,-4,3)

tau = sqrt( ((Sxx-Syy)/2).^2 + Sxy.^2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%