#include "cbox11.h"void Stress( Xp,Yp,X1,Y1,X2,Y2,dx11,dy11,dx12,dy12,             dx22,dy22,sx11,sy11,sx12,sy12,sx22,sy22)float Xp,Yp,X1,Y1,X2,Y2,*dx11,*dy11,*dx12,*dy12;float *dx22,*dy22,*sx11,*sy11,*sx12,*sy12,*sx22,*sy22; {	float Xg[5],Yg[5],Z[5],W[5]; 	/*[dimension increased by 1]*/	float Ax,Bx,Ay,By,nx,ny,slope,Perp,sgn,SR,FA,AL,Denom,rx,ry,Ra;	int i;	Z[1] = 0.86113631;	Z[2] = -Z[1];	Z[3] = 0.33998104;	Z[4] = -Z[3];	W[1] = 0.34785485;	W[2] = W[1];	W[3] = 0.65214515;	W[4] = W[3];		Ax = (X2-X1)/2.0;	Bx = (X2+X1)/2.0;	Ay = (Y2-Y1)/2.0;	By = (Y2+Y1)/2.0;	SR = sqrt(SQ(Ax)+SQ(Ay));	nx = (Y2-Y1)/(2*SR);	ny = (X1-X2)/(2*SR);		if( Ax ) 	{		slope = Ay/Ax;		Perp = fabs((slope*Xp-Yp+Y1-slope*X1)/sqrt(SQ(slope)+1));	}	else		Perp = fabs(Xp-X1);        	/*[Determine the direction of the outward normal ]*/	sgn = (X1-Xp)*(Y2-Yp)-(X2-Xp)*(Y1-Yp);	if ( sgn < 0 ) 		Perp = -Perp;        	(*dx11) = 0.0;	(*dy11) = 0.0;	(*dx12) = 0.0;	(*dy12) = 0.0;	(*dx22) = 0.0;	(*dy22) = 0.0;	(*sx11) = 0.0;	(*sy11) = 0.0;	(*sx12) = 0.0;	(*sy12) = 0.0;	(*sx22) = 0.0;	(*sy22) = 0.0;	/*[ Compute displacement and stress coefficients ]*/		FA = 1.0-4.0*nu;	AL = 1.0-2.0*nu;	Denom = 4.0*pi*(1.0-nu); 	           	for(i=1;i<=4;i++)	{		Xg[i] = Ax*Z[i]+Bx;		Yg[i] = Ay*Z[i]+By;		Ra  = sqrt(SQ(Xp-Xg[i])+SQ(Yp-Yg[i]));		rx = (Xg[i]-Xp)/Ra;		ry = (Yg[i]-Yp)/Ra;		(*dx11) += (AL*rx+2*cube(rx))*W[i]*SR/(Denom*Ra);		(*dy11) += (2*SQ(rx)*ry-AL*ry)*W[i]*SR/(Denom*Ra);		(*dx12) += (AL*ry+2*(SQ(rx))*ry)/(Denom*Ra)*W[i]*SR;		(*dy12) += (AL*rx+2*rx*SQ(ry))/(Denom*Ra)*W[i]*SR;		(*dx22) += (2*rx*SQ(ry)-AL*rx)/(Denom*Ra)*W[i]*SR;		(*dy22) += (AL*ry+2*cube(ry))/(Denom*Ra)*W[i]*SR;		(*sx11) += (2*Perp/Ra*(AL*rx+nu*2*rx-4*cube(rx))+                       4*nu*nx*SQ(rx)+AL*(2*nx*SQ(rx)+2*nx)-                       FA*nx)*2*mu/(Denom*SQ(Ra))*W[i]*SR;		(*sy11) += (2*Perp/Ra*(AL*ry-4*SQ(rx)*ry)+                       4*nu*nx*rx*ry+AL*2*ny*SQ(rx)-                       FA*ny)*2*mu/(Denom*SQ(Ra))*W[i]*SR;		(*sx12) += (2*Perp/Ra*(nu*ry-4*SQ(rx)*ry)+2*nu*                       (nx*ry*rx+ny*SQ(rx))+AL*(2*nx*rx*                       ry+ny))*2*mu/(Denom*SQ(Ra))*W[i]*SR;		(*sy12) += (2*Perp/Ra*(nu*rx-4*rx*SQ(ry))+2*nu*                       (nx*SQ(ry)+ny*rx*ry)+AL*(2*ny*rx*ry+                       nx))*2*mu/(Denom*SQ(Ra))*W[i]*SR;		(*sx22) += (2*Perp/Ra*(AL*rx-4*rx*SQ(ry))+4*nu*                       ny*rx*ry+AL*2*nx*SQ(ry)-FA*nx)*2*mu/                       (Denom*SQ(Ra))*W[i]*SR;		(*sy22) += (2*Perp/Ra*(AL*ry+2*nu*ry-4*cube(ry))+                       4*nu*ny*SQ(ry)+AL*(2*ny*SQ(ry)+2*ny)-                       FA*ny)*2*mu/(Denom*SQ(Ra))*W[i]*SR;     }    }