c--mribbon.f
c------------------------------------------------------------
      subroutine ribbon(x0,z0,x1,z1,x2,z2,mx,mz,fx,fz,ier)
c
c  Subroutine RIBBON computes the x and z components of magnetic 
c  induction due to a single side of a two-dimensional prism with 
c  polygonal cross section.  The prism is assumed infinitely 
c  extended parallel to the y axis; z axis is vertical down.
c
c  Input parameters:
c    Observation point is (x0,z0).  Coordinates (x1,z1) and 
c    (x2,z2) are two consecutive corners of the polygon taken in 
c    clockwise order around the polygon as viewed with the x 
c    axis to the right.  The x and z components of magnetization 
c    are mx and mz.  Distance units are irrelevant but must be 
c    consistent; magnetization in A/m.
c
c  Output parameters:
c    Components of magnetic field fx and fz, in nT. 
c    Errors are recorded by ier:  
c        ier=0, no errors; 
c        ier=1, two corners are too close (no calculation);
c        ier=2, field point too close to corner (calculation 
c               continues).
c
      real mx,mz
      data pi/3.14159265/,small/1.e-18/,cm/1.e-7/,t2nt/1.e9/
      ier=0
      sx=x2-x1
      sz=z2-z1
      s=sqrt(sx**2+sz**2)
c
c  --  If two corners are too close, return
c
      if (s.lt.small)then
         ier=1
         return
         end if
      sx=sx/s
      sz=sz/s
      qs=mx*sz-mz*sx
      rx1=x1-x0
      rz1=z1-z0
      rx2=x2-x0
      rz2=z2-z0
c
c  --  If field point is too near a corner, signal error
c
      if(abs(rx1).lt.small.and.abs(rz1).lt.small)ier=2
      if(abs(rx2).lt.small.and.abs(rz2).lt.small)ier=2
      if(ier.eq.2)then
         rx1=small
         rz1=small
         rx2=small
         rz2=small
         end if
      r1=sqrt(rx1**2+rz1**2)
      r2=sqrt(rx2**2+rz2**2)
      theta1=atan2(rz1,rx1)
      theta2=atan2(rz2,rx2)
      angle=theta1-theta2
      if (angle.gt.pi)angle=angle-2.*pi
      if (angle.lt.-pi)angle=angle+2.*pi
c
c  --  If field point is too close to side, signal error
c
      if (abs(angle).gt.(.995*pi))ier=2
      flog=alog(r2)-alog(r1)
      factor=-2.*cm*qs*t2nt
      fx=factor*(sx*flog-sz*angle)
      fz=factor*(sz*flog+sx*angle)
      return
      end