subroutine OnShl(enew, eold, Knew, Knew2, Dk2de)

c***********************************************************************
c     Onshl calc the on-shell momenta for every channel. 
c     The sign of the momenta (or which sheet of energies) is also
c          decided in this roution by looking the relative position 
c          of the current and previous momentum on the complex plane.
c     NOTE: Generally, Knew is complex because of complex physical mass
c           of Delta. To show 3D tmatrix struct, even the energy becomes
c	    complex which causes complex onshell Knew.
c***********************************************************************

      implicit none

      integer  nch, i
      double precision mtotal, mdiff, xold, yold, xnew, ynew, rek, imk
      double precision zero, two, three, four

      integer  Nifty, Isheet
      double precision mpion, Mass
      double complex enew, eold, Knew(2), Knew2(2), kold2(2), Dk2de(2)

      common /control/ Nifty(15)
      common /mass/ mpion, Mass(2)
      common /sheet/ Isheet(2)
      
      parameter (zero = 0.0d0, two = 2.0d0, three = 3.0d0, four = 4.0d0)
      
      nch = Nifty(2)
      if (eold .eq. 0.0d0) eold = mpion + Mass(1)

c     this if-statement serves as initialization for nifty(4)=3 
c     relativistic kimematics 
         do 10000 i = 1, nch
	    mtotal = mpion + Mass(i)
	    mdiff = mpion - Mass(i)
            Kold2(i) = (eold-mtotal)*(eold-mdiff)
     &                *(eold+mtotal)*(eold+mdiff)/(four*eold**2)
            Knew2(i) = (enew-mtotal)*(enew-mdiff)
     &                *(enew+mtotal)*(enew+mdiff)/(four*enew**2)
            Dk2de(i) = (enew**four - mtotal**two*mdiff**two)
     &                /(two*enew**three)
10000    continue

c the mometum k = sqrt(k**2), but the sign is undecided.
c we are try to decide the sign relative the the old k by looking where
c the straight line connecting Kold2 and Knew2 cross the 'cut'
c The 'cut' is choosen as from k=0 to +infinity alone positive real axis
c The real axis is belong to upper half plane
 
      do 10020 i = 1, nch
         xold = dble(Kold2(i))
         yold = imag(Kold2(i))
         xnew = dble(Knew2(i))
         ynew = imag(Knew2(i))
         if (yold*ynew .lt. zero) then
            if ((xold*ynew-xnew*yold)/(ynew-yold) .gt. zero)
     &         Isheet(i) = -Isheet(i)
         else if ((yold .eq. zero) .and. (xold .gt. zero)) then
            if (ynew .lt. zero) Isheet(i) = -Isheet(i)
         else if ((ynew .eq. zero) .and. (xnew .gt. zero)) then
            if (yold .lt. zero) Isheet(i) = -Isheet(i)
         endif
10020 continue

c compute channel momenta Knew(i) 
      do 10030 i = 1, nch
         Knew(i) = sqrt(Knew2(i))
         rek = dble(Knew(i))
         imk = imag(Knew(i))
         if (Isheet(i) .gt. 0) then
c force 0 <= arg(Knew(i) < Pi
            if (imk .lt. zero) then
               Knew(i) = -Knew(i)
            else if ((imk .eq. zero) .and. (rek .lt. zero)) then
               Knew(i) = -Knew(i)
            endif
         else
c force Pi <= arg(Knew(i) < 2Pi
            if (imk .gt. zero) then
               Knew(i) = -Knew(i)
            else if ((imk .eq. zero) .and. (rek .gt. zero)) then
               Knew(i) = -Knew(i)
            endif
         endif
10030 continue

      eold = enew

      return
      end