SUBROUTINE SUSDKM( p1, p2, p3, ps, iway, ierr )
C --------------------------------------------------------------------------
C   Purpose: Perform three-body decay of particle. The phase space
C            integral is re-done every 100 events, and is carried
C            out by a call to SUSBCH.
C
C     Notes: p1,2,3 = output 4 momenta.
C            ps     = 4 momenta of decay CMS frame in LAB.
C            iway   = decay type
C            ierr   < 0 skip event. 
C
C    Author: JMB 10-Mar-91
C    Last change 19/2/92
C --------------------------------------------------------------------------
C
      IMPLICIT NONE
C
      INCLUDE 'SUSANA.INC'
C
      INTEGER ierr, i, iway, nstore, itype
      REAL p1(4), p2(4), p3(4), ps(4), ECMS
      REAL mop1, mop2, mop3, ey1, ey2, ey3, rtheta, rphi
      REAL randum, rn32
      REAL E12(2), temp,  phase, pherr 
      REAL lowe(2), highe(2)
      DOUBLE PRECISION beta, alph
      DOUBLE PRECISION theta(3), phi(3)
C
C
C -- ! Start in CMS frame, positron along z-axis ! --
C
      ecms = msus(8)
C     Identify decay class.
      IF (iway.LE.2) THEN 
        itype = 1
      ELSE
        itype = 2
      ENDIF
C
C -----> Pick out energies from store.
C
C -   Generate more energies if necessary
      IF (e_used(itype).GE.100) THEN
        IF ((nevent-event_done).GE.100) THEN
          nstore = 100
        ELSE
          nstore = nevent - event_done
        ENDIF
        IF (nstore.GT.0) THEN
          CALL SUSBCH(1, nstore, ierr)
        ELSE
          ierr = 1
          WRITE(isustr(1),*) 'ABNORMAL END!'
          GOTO 8900
        ENDIF
      ENDIF
C
      e_used(itype) = e_used(itype) + 1
C     Load energies
      IF (itype.EQ.1) THEN
        Ey1 = e_e( 2*e_used(itype) - 1 )
        Ey2 = e_e( 2*e_used(itype) )
      ELSE
        Ey1 = e_nu( 2*e_used(itype) - 1 )
        Ey2 = e_nu( 2*e_used(itype) )
      ENDIF
C
C --  Ey1 & E2 are now fixed. Calculate E3.
      Ey3 = ECMS - Ey1 - Ey2
C
C --  Calculate momenta of particles.
      IF ( ((Ey1**2-m1**2).LT.0.0) .OR. ((Ey2**2-m2**2).LT.0.0) .OR. 
     &     ((Ey3**2-m3**2).LT.0.0) ) THEN
        IF (isustr(3).GT.0) THEN
          WRITE(isustr(1),*) 'SUSDKM:Neg.Sqrt:Skipping an event'
        ENDIF
        ierr = -5
        GOTO 8900
      ENDIF
      mop1 = sqrt(Ey1**2 - m1**2)
      mop2 = sqrt(Ey2**2 - m2**2)
      mop3 = sqrt(Ey3**2 - m3**2)
C
      IF (mop1.LE.0.0) mop1 = small
      IF (mop2.LE.0.0) mop2 = small
      IF (mop3.LE.0.0) mop3 = small
      temp = (mop2**2-mop3**2-mop1**2)/(2.0*mop1*mop3) 
      IF ( abs( temp ).GT.1.0) 
     &   THEN
        IF (isustr(3).GT.0) THEN
          write(isustr(1),*) 
     &   'SUSDKM6:Unphysical momenta:skipping event'
        ENDIF
        IERR = -6
        GOTO 8900
      ENDIF
C
C --  Calculate the angles of the tracks.
      beta  = acos( temp )
      IF (abs((mop1/mop2 )*sin(beta)).GT.1.0) THEN
        IF (isustr(3).GT.0) THEN
          write(isustr(1),*) 
     &   'SUSDKM:Unphysical momenta - skipping event'
        ENDIF
        ierr = -2
        GOTO 8900
      ENDIF        
      alph = asin( -(mop1/mop2 )*sin(beta))
      temp = (-mop3 - mop1*temp) / mop2
      IF (temp.GT.0.0) THEN
C --    alpha is OK.
      ELSE
        alph = -pi - alph
      ENDIF
      alph = twopi + alph
C
C --> Thus far, the z-axis was defined by the direction of E3.
C     Now rotate all tracks by a random theta, phi.
C
      rphi   = rn32(randum)*pi
      rtheta = rn32(randum)*twopi
C
      theta(1) = beta  + rtheta
      theta(2) = alph  + rtheta
      theta(3) = rtheta
      IF (beta.GT.twopi) beta = beta - twopi
      IF (alph.GT.twopi) alph = alph - twopi
C
C
 2    CONTINUE
C --> Now the z-axis is the ZEUS z axis. 
C
C --  Lorentz boost by velocity of photino along the photino direction.
C     Go to cartesians...
      p1(1) = mop1 * sin(theta(1)) * abs(cos(rphi))
      p1(2) = mop1 * sin(theta(1)) * abs(sin(rphi))
      p1(3) = mop1 * cos(theta(1))
      p1(4) = Ey1
C
      p2(1) = mop2  * sin(theta(2)) * abs(cos(rphi))
      p2(2) = mop2  * sin(theta(2)) * abs(sin(rphi))
      p2(3) = mop2  * cos(theta(2))
      p2(4) = Ey2
C
      p3(1) = mop3 * sin(theta(3)) * abs(cos(rphi))
      p3(2) = mop3 * sin(theta(3)) * abs(sin(rphi))
      p3(3) = mop3 * cos(theta(3))
      p3(4) = Ey3
C
      CALL LORENB( ECMS, ps, p1, p1 )
      CALL LORENB( ECMS, ps, p2, p2 )
      CALL LORENB( ECMS, ps, p3, p3 )
C
 8900 CONTINUE
      RETURN
      END