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