c ****** 1 ******* 2 ******* 3 ******* 4 ******* 5 ******* 6 ******* 7 *
c PROGRAM TO FIND THE "CUT VOLUME" V0 GIVEN r0, dr0 AND
c m1 x1 + m2 x2 + m3 x3 = alpha
c ****** 1 ******* 2 ******* 3 ******* 4 ******* 5 ******* 6 ******* 7 *
      DOUBLE PRECISION FUNCTION FL3D(m1,m2,m3,alpha,r0,dr0)
c***
      INCLUDE 'undefined.h'
      DOUBLE PRECISION m1,m2,m3,alpha,r0,dr0
      DOUBLE PRECISION al,al0,n1,n2,n3,b1,b2,b3,b12,bm,tmp,pr
      INTRINSIC DMAX1,DMIN1,DABS
c***
c     (1) move origin to r0 along r ;  (2) reflect parallelepiped;
c     (3) limit alpha (0<= al0 <=0.5); (4) order coefficients: b1<b2<b3;
c     (5) calculate volume (NOTE: it is assumed:s0=t0=0; ds0=dt0=1.)
c*(1)*
      al = alpha - m1*r0
c*(2)*
      al = al + DMAX1(0.d0,-m1*dr0)+DMAX1(0.d0,-m2)+DMAX1(0.d0,-m3)
      tmp = DABS(m1)*dr0 + DABS(m2) + DABS(m3)
      n1 = DABS(m1)/tmp
      n2 = DABS(m2)/tmp
      n3 = DABS(m3)/tmp
      al = DMAX1(0.d0,DMIN1(1.d0,al/tmp))
c*(3)*
      al0 = DMIN1(al,1.d0-al)
c*(4)*
      b1 = DMIN1(n1*dr0,n2)
      b3 = DMAX1(n1*dr0,n2)
      b2 = n3
      if (b2 .LT. b1) then
         tmp = b1
         b1 = b2
         b2 = tmp
      else if (b2 .GT. b3) then
         tmp = b3
         b3 = b2
         b2 = tmp
      endif
      b12 = b1 + b2
      bm = DMIN1(b12,b3)
      pr = DMAX1(6.d0*b1*b2*b3,1.0d-50)
c*5*     
      if (al0 .LT. b1) then
         tmp = al0*al0*al0/pr
      else if (al0 .LT. b2) then
         tmp = 0.5d0*al0*(al0-b1)/(b2*b3) +  b1*b1*b1/pr
      else if (al0 .LT. bm) then
         tmp = (al0*al0*(3.d0*b12-al0) + b1*b1*(b1-3.d0*al0) +
     %        b2*b2*(b2-3.d0*al0))/pr
      else if (b12 .LT. b3) then
         tmp = (al0 - 0.5d0*bm)/b3
      else
         tmp = (al0*al0*(3.d0-2.d0*al0) + b1*b1*(b1-3.d0*al0) +
     %        b2*b2*(b2-3.d0*al0) + b3*b3*(b3-3.d0*al0))/pr
      endif
      
      if (al .LE. 0.5d0) then
         FL3D = tmp*dr0
      else
         FL3D = (1.d0-tmp)*dr0
      endif
c***  
      return
      end