subroutine facmag(mx,my,mz,x0,y0,z0,x,y,z,n,fx,fy,fz) c c Subroutine FACMAG computes the magnetic field due to surface c charge on a polygonal face. Repeated calls can build the c field of an arbitrary polyhedron. X axis is directed north, c z axis vertical down. Requires subroutines CROSS, ROT, LINE, c and PLANE. Algorithm from Bott (1963). c c Input parameters: c Observation point is (x0,y0,z0). Polygon corners defined c by arrays x, y, and z of length n. Magnetization given by c mx,my,mz. Polygon limited to 10 corners. Distance units c are irrelevant but must be consistent; magnetization in A/m. c c Output parameters: c Three components of magnetic field (fx,fy,fz), in nT. c real mx,my,mz,nx,ny,nz dimension u(10),v2(10),v1(10),s(10),xk(10),yk(10), & zk(10),xl(10),yl(10),zl(10),x(10),y(10),z(10) data cm/1.e-7/,t2nt/1.e9/,epsilon/1.e-20/ fx=0. fy=0. fz=0. x(n+1)=x(1) y(n+1)=y(1) z(n+1)=z(1) do 1 i=1,n xl(i)=x(i+1)-x(i) yl(i)=y(i+1)-y(i) zl(i)=z(i+1)-z(i) rl=sqrt(xl(i)**2+yl(i)**2+zl(i)**2) xl(i)=xl(i)/rl yl(i)=yl(i)/rl 1 zl(i)=zl(i)/rl call cross(xl(2),yl(2),zl(2),xl(1),yl(1),zl(1),nx,ny,nz,rn) nx=nx/rn ny=ny/rn nz=nz/rn dot=mx*nx+my*ny+mz*nz if(dot.eq.0.)return call plane(x0,y0,z0,x(1),y(1),z(1),x(2),y(2),z(2),x(3), & y(3),z(3),px,py,pz,w) do 2 i=1,n call rot(x(i),y(i),z(i),x(i+1),y(i+1),z(i+1),nx,ny,nz, & px,py,pz,s(i)) if(s(i).eq.0.)go to 2 call line(px,py,pz,x(i),y(i),z(i),x(i+1),y(i+1),z(i+1), & u1,v,w1,v1(i),v2(i),u(i)) rk=sqrt((u1-px)**2+(v-py)**2+(w1-pz)**2) xk(i)=(u1-px)/rk yk(i)=(v-py)/rk 2 zk(i)=(w1-pz)/rk do 3 j=1,n if(s(j).eq.0.)go to 3 us=u(j)**2 v2s=v2(j)**2 v1s=v1(j)**2 a2=v2(j)/u(j) a1=v1(j)/u(j) f2=sqrt(1.+a2*a2) f1=sqrt(1.+a1*a1) rho2=sqrt(us+v2s) rho1=sqrt(us+v1s) r2=sqrt(us+v2s+w**2) r1=sqrt(us+v1s+w**2) if(w.ne.0.)then fu2=(a2/f2)*alog((r2+rho2)/abs(w)) & -.5*alog((r2+v2(j))/(r2-v2(j))) fu1=(a1/f1)*alog((r1+rho1)/abs(w))- & .5*alog((r1+v1(j))/(r1-v1(j))) fv2=(1./f2)*alog((r2+rho2)/abs(w)) fv1=(1./f1)*alog((r1+rho1)/abs(w)) fw2=atan2((a2*(r2-abs(w))),(r2+a2*a2*abs(w))) fw1=atan2((a1*(r1-abs(w))),(r1+a1*a1*abs(w))) fu=dot*(fu2-fu1) fv=-dot*(fv2-fv1) fw=(-w*dot/abs(w))*(fw2-fw1) else fu2=(a2/f2)*(1.+alog((r2+rho2)/epsilon))- & .5*alog((r2+v2(j))/(r2-v2(j))) fu1=(a1/f1)*(1.+alog((r1+rho1)/epsilon))- & .5*alog((r1+v1(j))/(r1-v1(j))) fv2=(1./f2)*(1.+alog((r2+rho2)/epsilon)) fv1=(1./f1)*(1.+alog((r1+rho1)/epsilon)) fu=dot*(fu2-fu1) fv=-dot*(fv2-fv1) fw=0. end if fx=fx-s(j)*(fu*xk(j)+fv*xl(j)+fw*nx) fy=fy-s(j)*(fu*yk(j)+fv*yl(j)+fw*ny) fz=fz-s(j)*(fu*zk(j)+fv*zl(j)+fw*nz) 3 continue fx=fx*cm*t2nt fy=fy*cm*t2nt fz=fz*cm*t2nt return end