function [area,centroid,Z,fss,fds,fess,feds] = calcp(panel,evalpnts,directions) % Matlab version of calcp, returns potential at evaluation point due % to unit monopole and unit dipole uniformly distributed on a panel. % Follows a left-hand rule (Clockwise ordered points has normal % pointing up). % panel -- vectors of panel vertices in rows of x,y,z (3 or 4 rows supported). % evalpnts -- matrix of evaluation points, rows of x,y,z coordinates % directions -- matrix of derivative directions, rows of x,y,z coordinates % fss = the vector of potentials due to a monopole % fds = the vector of potentials due to a panel normal dipole distribution % area = panel area. % centroid = panel centroid. % Z = panel normal. % fess = the derivative of the monopole potential at evalpnt along direction % fess = the derivative of the dipole potential at evalpnt along direction % % First check the input. [verts, betterbethree] = size(panel); if betterbethree ~= 3 'wrong panel format: should be rows of x,y,z vectors!' return; end if (verts > 4) | (verts < 2) 'wrong panel format: panel can only have 3 or 4 vertices!' return; end if(nargin > 1) [numevals, betterbethree] = size(evalpnts); if betterbethree ~= 3 'wrong evaluation point format: should be rows of x,y,z vectors!' return; end else numevals = 0; end if (nargin > 2) deriv = 1; [numdirections, betterbethree] = size(directions); if betterbethree ~= 3 'wrong direction vector format: should be rows of x,y,z vectors!' return; end if (numdirections ~= 0) & (numdirections ~= numevals) 'number of direction vectors does not match number of evaluation points!' return; end else deriv = 0; end % The PANEL SETUP!!!!!!!!!*************************************************** % Length of each side and the panel area. for i=1:verts if(i < verts) side(i,:) = panel(i+1,:) - panel(i,:); else side(i,:) = panel(1,:) - panel(i,:); end edgeLength(i) = norm(side(i,:)); end % Calculate the panel coordinate system. X = panel(3,:) - panel(1,:); diagLength = norm(X); if(verts == 3) Y = panel(2,:) - panel(1,:); else Y = panel(2,:) - panel(4,:); end % Z-axis is normal to two diags. Z = cross(X, Y); area = 0.5 * norm(Z); % Normalize panel axises. coord(3,:) = Z / norm(Z); coord(1,:) = X / norm(X); X = coord(1,:); Z = coord(3,:); coord(2,:) = cross(Z, X); Y = coord(2,:); % Determine the centroid. vertex1 = panel(2,:) - panel(1,:); if(verts == 4) vertex3 = panel(4,:) - panel(1,:); else vertex3 = panel(3,:) - panel(1,:); end % Project into the panel axes. y1 = sum(vertex1 .* Y); y3 = sum(vertex3 .* Y); x1 = sum(vertex1 .* X); x3 = sum(vertex3 .* X); yc = (y1 + y3)/3.0; xc = (diagLength + ((x1 * y1 - x3 * y3)/(y1 - y3)))/3.0; % Compute the centroid. centroid = panel(1,:) + xc * X + yc * Y; % Put the corners in the newly defined coordinate system. for i=1:verts npanel(i,:) = (coord * (panel(i,:) - centroid).').'; end % Check that panel is in the x-y panel. for i=1:verts if(abs(npanel(i,3)) > (1.0e-8 * diagLength)) 'Coordinate transform failure!!' npanel return; end; end; % Compute the contributions terms for each edge. for i=1:verts if (i==verts) next=1; else next=i+1; end; ct(i) = (npanel(next,1)-npanel(i,1))/edgeLength(i); st(i) = (npanel(next,2)-npanel(i,2))/edgeLength(i); end % Done with the PANEL SETUP!!!!!!!!!************************************ % Done with Setup, now loop through the evaluation points!! for evalindex = 1:numevals % Point eval pt in in new coordinate system and get the z-comp. point = (coord * ((evalpnts(evalindex,:) - centroid).')).'; if(deriv == 1) nrm = (coord * directions(evalindex,:).').'; end zn=point(3); znabs=abs(zn); evalDistance = norm(point); % Once per vertex computation OK = 1; for i=1:verts xc=point(1)-npanel(i,1); yc=point(2)-npanel(i,2); zc=point(3)-npanel(i,3); xmxv(i)=xc; ymyv(i)=yc; fe(i)=xc*xc+zc*zc; r(i)=sqrt(yc*yc+fe(i)); if (r(i) < (1.005*znabs)) OK = 0; end; if(deriv == 1) xri(i) = xmxv(i)/r(i); yri(i) = ymyv(i)/r(i); end; end; % The potential and dipole contributions are made by summing up % a contribution from each edge fs=0; fd=0; if(deriv == 1) fsx = 0; fsy = 0; fdx = 0; fdy = 0; fdz = 0; end for i=1:verts if (i==verts) next=1; else next=i+1; end; % v is the projection of the eval-i edge on the perpend to the side-i: % Exploits the fact that corner points in panel coordinates. v=xmxv(i)*st(i) - ymyv(i)*ct(i); % arg == zero if eval on next-i edge, but then v = 0. % arg=(r(i)+r(next)-edgeLength(i))/(r(i)+r(next)+edgeLength(i)); if(arg == 0) 'in calcp' keyboard; end fln = -log(arg); if (arg>0.0) fs = fs + v * fln; end; if ( deriv == 1 ) if ( arg > 0.0 ) fac = (r(i)+r(next)-edgeLength(i)) * (r(i)+r(next)+edgeLength(i)); fac = v*(edgeLength(i)+ edgeLength(i))/fac; fsx = fsx + (fln*st(i) - fac*(xri(i) + xri(next))); fsy = fsy - (fln*ct(i) + fac*(yri(i) + yri(next))); fdz = fdz - (fac*( 1.0/r(i) + 1.0/r(next))); end end % OK means eval not near a vertex normal, use Hess-Smith: if (OK == 1) s1=v*r(i); c1=znabs*(xmxv(i)*ct(i)+ymyv(i)*st(i)); s2=v*r(next); c2=znabs*(xmxv(next)*ct(i)+ymyv(next)*st(i)); else % Near a vertex normal, use Newman s1=(fe(i)*st(i))-(xmxv(i)*ymyv(i)*ct(i)); c1=znabs*r(i)*ct(i); s2=(fe(next)*st(i))-(xmxv(next)*ymyv(next)*ct(i)); c2=znabs*r(next)*ct(i); end; s12=(s1*c2)-(s2*c1); c12=(c1*c2)+(s1*s2); val=atan2(s12, c12); fd=fd+val; if (deriv == 1) u1 = xmxv(i)*ct(i) + ymyv(i)*st(i); u2 = xmxv(next)*ct(i)+ymyv(next)*st(i); if (OK == 0) % Near a vertex normal. rr = r(i)*r(i); fh1 = xmxv(i)*ymyv(i); fh2 = xmxv(next)*ymyv(next); fac = c1/((c1*c1+s1*s1)*rr ); if(zn < 0.0) fac = -1.0 * fac; end fdx = fdx + ((rr*v+fh1*u1)*fac); fdy = fdy - (fe(i)*u1*fac); rr = r(next)*r(next); fac = c2/((c2*c2+s2*s2)*rr); if(zn < 0.0) fac = -1.0 * fac; end fdx = fdx - ((rr*v+fh2*u2)*fac); fdy = fdy + fe(next)*u2*fac; else fac = zn/(c1*c1+s1*s1); fdx = fdx + (u1*v*xri(i)+r(i)*ymyv(i))*fac; fdy = fdy + (u1*v*yri(i)-r(i)*xmxv(i))*fac; fac = zn/(c2*c2+s2*s2); fdx = fdx - ((u2*v*xri(next)+r(next)*ymyv(next))*fac); fdy = fdy - ((u2*v*yri(next)-r(next)*xmxv(next))*fac); end end; end; if (fd<0.0) fd = fd + 2*pi; end; if (zn < 0) fd=fd*(-1.0); end; fs=fs-zn*fd; if(deriv == 1 ) fsx = fsx - zn*fdx; fsy = fsy - zn*fdy; fes = nrm(1)*fsx + nrm(2)*fsy - nrm(3)*fd; fed = nrm(1)*fdx + nrm(2)*fdy + nrm(3)*fdz; end % No area normalization! fss(evalindex) = fs; fds(evalindex) = fd; if(deriv == 1) fess(evalindex) = fes; feds(evalindex) = fed; end; end;