#include "cbox1.h"void Quad1(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;/*[ Quad1 computes the integral of several nonsingular functions along the boundary elements using a four points Gauss Quadrature.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, as defined in section 5.3 ]*/{const float Z[] =  {0.0,0.861136311594053,-0.861136311594053,0.339981043584856,-0.339981043584856};const float W[] = {0.0,0.347854845137454,0.347854845137454,0.652145154862546,0.652145154862546};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 coefficients ]*/  						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 Diag1(X1,Y1,X2,Y2,G)float X1,Y1,X2,Y2,*G;{  float Ax,Ay,HL;  Ax = (X2 - X1)/2.0;  Ay = (Y2 - Y1)/2.0;  HL = sqrt(SQ(Ax) + SQ(Ay));  (*G) = 2* HL*(1.0 - log(HL)); }