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