c/* %W%      latest revision %G% %U% */
      subroutine DAMPS (x, ApB, AmB, CpD, CmD, E)
************************************************************************

c *** Special sub to avoid div by zero/noise in calculation of D
c *** Called instead of TPN.
c *** If x = cos(theta) > 0.3 then DAMPS simply passes everything
c *** over to TPN.  In the case that x is NOT less greater than
c *** 0.3 then DAMPS will still call on TPN to calculate the
c *** amplitudes.  However, DAMPS will then modify the values of
c *** D, C+D, and C-D before return to the calling routine.


      implicit real*8 (a-h, i, k, o-z)

              character*80 SCCSID

      complex*16  Ap, An, Bp, Bn, Cp, Cn, Dp, Dn, Ep, En
      complex*16  Dp1, Dp2, Dn1, Dn2

      dimension   ApB(4), AmB(4), CpD(4), CmD(4), E(4)
      dimension   CpD1(4), CpD2(4), CmD1(4), CmD2(4)

      common /amplit/   Ap, An, Bp, Bn, Cp, Cn, Dp, Dn, Ep, En


c>>>> FIRST EXECUTABLE STATEMENT <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

      if (abs(x) .gt. 0.3) go to 899
      x1 = -0.3
      x2 = +0.3

      call TPN (x1, ApB, AmB, CpD1, CmD1, E)

      Dp1 = Dp
      Dn1 = Dn

      call TPN (x2, ApB, AmB, CpD2, CmD2, E)

      Dp2 = Dp
      Dn2 = Dn

c *** Scale with proportional parts

 898  p2 = (x - x1)/(x2 - x1)
      p1 = (x2 - x)/(x2 - x1)

c *** Call TPN to reset other, non d, amplitudes

 899  call TPN (x, ApB, AmB, CpD, CmD, E)

      if (abs(x) .gt. 0.3) then
         return
      endif

      do 897 j = 1,4
         CpD(j) = p1 * CpD1(j) + p2 * CpD2(j)
         CmD(j) = p1 * CmD1(j) + p2 * CmD2(j)
 897  continue

      Dp = p1 * Dp1 + p2 * Dp2
      Dn = p1 * Dn1 + p2 * Dn2
      return
      end