subroutine mbox(x0,y0,z0,x1,y1,z1,x2,y2,mi,md,fi,fd,m,theta,t)
c
c  Subroutine MBOX computes the total field anomaly of an infinitely 
c  extended rectangular prism.  Sides of prism are parallel to x,y,z 
c  axes, and z is vertical down.  Bottom of prism extends to infinity.
c  Two calls to mbox can provide the anomaly of a prism with finite 
c  thickness; e.g., 
c
c     call mbox(x0,y0,z0,x1,y1,z1,x2,y2,mi,md,fi,fd,m,theta,t1)
c     call mbox(x0,y0,z0,x1,y1,z2,x2,y2,mi,md,fi,fd,m,theta,t2)
c     t=t1-t2
c
c  Requires subroutine DIRCOS.  Method from Bhattacharyya (1964).
c
c  Input parameters:
c    Observation point is (x0,y0,z0).  Prism extends from x1 to 
c    x2, y1 to y2, and z1 to infinity in x, y, and z directions, 
c    respectively.  Magnetization defined by inclination mi, 
c    declination md, intensity m.  Ambient field defined by 
c    inclination fi and declination fd.  X axis has declination 
c    theta. Distance units are irrelevant but must be consistent.  
c    Angles are in degrees, with inclinations positive below 
c    horizontal and declinations positive east of true north.
c    Magnetization in A/m.
c
c  Output paramters:
c    Total field anomaly t, in nT.
c
      real alpha(2),beta(2),mi,md,m,ma,mb,mc
      data cm/1.e-7/,t2nt/1.e9/
      call dircos(mi,md,theta,ma,mb,mc)
      call dircos(fi,fd,theta,fa,fb,fc)
      fm1=ma*fb+mb*fa
      fm2=ma*fc+mc*fa
      fm3=mb*fc+mc*fb
      fm4=ma*fa
      fm5=mb*fb
      fm6=mc*fc
      alpha(1)=x1-x0
      alpha(2)=x2-x0
      beta(1)=y1-y0
      beta(2)=y2-y0
      h=z1-z0
      t=0.
      hsq=h**2
      do 1 i=1,2
         alphasq=alpha(i)**2
         do 1 j=1,2
            sign=1.
            if(i.ne.j)sign=-1.
            r0sq=alphasq+beta(j)**2+hsq
            r0=sqrt(r0sq)
            r0h=r0*h
            alphabeta=alpha(i)*beta(j)
            arg1=(r0-alpha(i))/(r0+alpha(i))
            arg2=(r0-beta(j))/(r0+beta(j))
            arg3=alphasq+r0h+hsq
            arg4=r0sq+r0h-alphasq
            tlog=fm3*alog(arg1)/2.+fm2*alog(arg2)/2.
     &            -fm1*alog(r0+h)
            tatan=-fm4*atan2(alphabeta,arg3)
     &            -fm5*atan2(alphabeta,arg4)
     &            +fm6*atan2(alphabeta,r0h)
    1 t=t+sign*(tlog+tatan)
      t=t*m*cm*t2nt
      return
      end