#include "cbox5.h"void Quad5(Xp,Yp,X1,Y1,X2,Y2,H,G,qx,qy,ux,uy,K)	float Xp,Yp,X1,Y1,X2,Y2,*H,*G;	float *qx,*qy,*ux,*uy;	int K;{    /*[ Program Quad5 :This function computes the integral of several nonsingular functions along the boundary elements using a four points Gauss Quadrature.If K = 0, the off-diagonal coefficients of the H and G matrices are computed; when K = 1, all coefficients needed to compute the potential and fluxes at the interior points are computed.Ra = Radius = Distance from the collocation point to the Gauss Integration points on the boundary element; nx,ny = components of the unit normal to the element; rx,ry,rn = Radius derivatives       ]*/const float Z[] = { 0.0,0.86113631,-0.86113631,0.33998104,-0.33998104};const float W[] = {0.0,0.34785485, 0.34785485,0.65214515, 0.65214515};float Xg[5],Yg[5];float Ax,Bx,Ay,By,HL,nx,ny,Ra,rx,ry,rn;int i;	Ax = (X2 - X1)/2;	Bx = (X2 + X1)/2;	Ay = (Y2 - Y1)/2;	By = (Y2 + Y1)/2;	HL = sqrt(SQ(Ax) + SQ(Ay));	nx =  Ay/HL;	ny = -Ax/HL;	(*G) = 0.0;	(*H) = 0.0;	(*ux)  = 0.0;	(*uy)  = 0.0;	(*qx)  = 0.0;	(*qy)  = 0.0;/*[ Compute G[x][y], H[x][y], qx, qy, ux, and uy ]*/  							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;		  rn = rx*nx + ry*ny;		  if(K > 0)		    {			  (*ux) += rx*W[i]*HL/Ra;			  (*uy) += ry*W[i]*HL/Ra;			  (*qx) -= (( 2.0 * SQ(rx) -1.0)*nx + 2.0 *rx*ry*ny)*			  		  	W[i]*HL/SQ(Ra);			  			  (*qy) -= (( 2.0 * SQ(ry) -1.0)*ny + 2.0 *rx*ry*nx)*			  		  	W[i]*HL/SQ(Ra);			 }			 			  (*G) += log(1/Ra) * W[i]*HL;			  (*H) -= rn*W[i]*HL/Ra;		  }}void Diag5(X1,Y1,X2,Y2,G)float X1,Y1,X2,Y2,*G;{  float Ax,Ay,Li;  Ax = (X2 - X1)/2.0;  Ay = (Y2 - Y1)/2.0;  Li = sqrt(SQ(Ax) + SQ(Ay));  (*G) = 2* Li*(1.0 - log(Li)); }