c/* %W%      latest revision %G% %U% */
      subroutine TNCM(kappa)
***********************************************************************

c *** Subroutine TNCM calcs averaged on-shell pN c.m. amplitudes
c ***
c *** kappa(or vector kap) = p-nucleon c.m. momentum
c ***
c *** Initialize if phase shifts or nucleus (Z, A) changes

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

c ********* tvm 2/24/92  change all dimensions of 100 to 160 ********
      dimension   nifty(20), Ecm(160), Tlab(160), kap(160), plab(160)
      dimension   d(40,160), eta(40,160), eps(8,160)
      dimension   t(40), tx(40), rewave(40,160), imwave(40,160)
      dimension   bn(40), bnuc(40), bnf(40), bnucf(40)
CCC	increase dimensionality to include other off-diagonal
CCC	ra(7)--> ra(14); r(7)-->r(14) : pjd 4/4/93
c      dimension   ta(7), tax(7), ra(7), rax(7)
      dimension   ta(7), tax(7), ra(14), rax(14)
CCC
      dimension   tawave(7,160), taxwve(7,160), retij(40), imtij(40)

      common /tjmats/   t, tx, ra, rax
      common /bandt/    bnuc, bn, bnucf, bnf, retij, imtij
      common /params/   hbarc, pi, Mp, Mn, nZ, nA, nes, nwaves
      common /switch/   nifty
      common /tape/     rewave, imwave, tawave, taxwve

      equivalence    (rewave(1,1), d(1,1)),  (imwave(1,1), eta(1,1))

c     data ndata/1/, nfold/0/

c *** Flip is more complicated, see Alexander's notes
c *** Nwave code (order for read in changed from Alex)
c ***
c *** I = 1
c ***
c ***  1) 1S0     2) 1D2     3) 1G4     4) 3P0     5) 3P1     6) 3P2
c ***  7) 3F2     8) 3F3     9) 3F4    10) 3H4    11) 3H5    12) 3H6
c ***
c ***
c *** I = 0
c ***
c *** 13) 1S0    14) 1P1    15) 1F3    16) 1H5    17) 3S1    18) 3D1
c *** 19) 3D2    20) 3D3    21) 3G3    22) 3G4    23) 3G5
c ***
c ***
c *** I = 1
c ***
c *** 24) 1I6    25) 3J6    26) 3J7
c ***
c ***
c *** I = 0
c ***
c *** 27) 3I5    28) 3I6    29) 3I7    30) 1J7    31) 3K7    32) 3K8
c ***

c *** After read in phases, the mixing parameters follow
c *** epsilon (eps) 1,2,3,.....

c *** N.B.  This is modified form from Alexander

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

      Mln = Mp
      nucln = 1
      ndata = 0
      IF (kappa .le. -100.0) ndata = 1
      kappa = abs(kappa)
      IF (ndata .eq. 0) go to 2001
      IF (nifty(8) .ne. 1) go to 2000

c *** Nifty(8) = 1, then read in phase shifts.
c *** Tmatrix (unfolded) for this nucleus.
c *** Amplitudes from input phase shifts and energy.
c ***
c *** Tlab in MeV, delta (d) read in in degrees, f calculated in MeV-1,
c *** eta = exp(-2 * Im(delta)), if read in as 0 then set to 1.
c ***

      open(1, file='tape1', status='unknown')
      nei = 1
      write(6,902) nes, nwaves
      rewind 1
      read(1,919) ne, nwaves

c *** Proton or neutron

      do 8011 ne = nei, nes
         read(1,9001) Tlab(ne),
     $                (eta(nwave,ne), d(nwave,ne), nwave=1,nwaves)
         read(1,9001) (eps(j,ne), j=1,8)

c *** Write out phases

         if( nifty(20) .eq. 1 .OR. nifty(20) .eq. 3 ) then
            write(6,9003) Tlab(ne),
     $                    (eta(nwave,ne), d(nwave,ne), nwave=1,nwaves)
            write(6,9003) Tlab(ne), (eps(j,ne), j=1,8)
         else if( ne .eq. 1 .OR. ne .eq. nes ) then
c                 write(6,9003) Tlab(ne),
c    $                    (eta(nwave,ne), d(nwave,ne), nwave=1,nwaves)
c                 write(6,9003) Tlab(ne), (eps(j,ne), j=1,8)
         endif

c *** Lab KE read in, convert to Ecm

         Ecm(ne) = sqrt(2. * Mln * (Tlab(ne) + 2. * Mln))
 8011 continue

      rewind 1
      do 548 nu = 1,7
         ta(nu) = 0.
         tax(nu) = 0.
 548  continue

c-------      Calculate kinematics

      do 1000 ne = 1,nes
         s = Ecm(ne) * Ecm(ne)
         kap(ne) = SQRT((s - (Mln + Mp) * (Mln + Mp))
     $                  * (s - (Mln - Mp) * (Mln - Mp))/4./s)
         plab(ne)= kap(ne) * SQRT(s)/Mln
         k = kap(ne)
         Epofk = SQRT( Mp*Mp + k*k )
         ENofk  = SQRT( Mln*Mln + k*k )
         ftot   = -(Epofk + ENofk)/4./pi/pi/Epofk/ENofk
         factor = ftot/k

c-----------  Compute alphas and convert deltas
c ***
c *** See LPT p. 311
c *** See also Stapp, Ypsilantis, Metropolis Phys Rev 105 #1,305/57
c ***
c *** Code for ta, tax p. 305 Stapp et al..
c ***          ta(1) = Re(alpha1/2i), ta(2) = Re(alpha2/2i) etc.
c ***          tax for imaginary parts of the same alphas.
c *** A factor of 1/2, not in the formula, included to agree
c *** with rewave imwave in TNCM then one factor of 1/2 not included
c *** in the formulas for t11, t1m1, t00, and tss.

         ds  = (d(18, ne) + d(17, ne)) * pi/180.
         fp  = (d( 7, ne) + d( 6, ne)) * pi/180.
         gd  = (d(21, ne) + d(20, ne)) * pi/180.
         hf  = (d(10, ne) + d( 9, ne)) * pi/180.
         ig  = (d(27, ne) + d(23, ne)) * pi/180.
         ajh = (d(25, ne) + d(12, ne)) * pi/180.
         aki = (d(31, ne) + d(29, ne)) * pi/180.

         do 544 nwave = 1,nwaves
            IF (eta(nwave, ne) .le. 0.) eta(nwave, ne) = 1.
  544    continue

         eta21 = SQRT( eta(18, ne) )
         eta01 = SQRT( eta(17, ne) )
         eta32 = SQRT( eta( 7, ne) )
         eta12 = SQRT( eta( 6, ne) )
         eta43 = SQRT( eta(21, ne) )
         eta23 = SQRT( eta(20, ne) )
         eta54 = SQRT( eta(10, ne) )
         eta34 = SQRT( eta( 9, ne) )
         eta65 = SQRT( eta(27, ne) )
         eta45 = SQRT( eta(23, ne) )
         eta76 = SQRT( eta(25, ne) )
         eta56 = SQRT( eta(12, ne) )
         eta87 = SQRT( eta(31, ne) )
         eta67 = SQRT( eta(29, ne) )
         et1s  = eta21 * eta01
         et2s  = eta32 * eta12
         et3s  = eta43 * eta23
         et4s  = eta54 * eta34
         et5s  = eta65 * eta45
         et6s  = eta76 * eta56
         et7s  = eta87 * eta67

c *** To make s = sin(2 * epsilon * pi/180)

         s1     = SIN( eps(1, ne) * pi/90. )
         s2     = SIN( eps(2, ne) * pi/90. )
         s3     = SIN( eps(3, ne) * pi/90. )
         s4     = SIN( eps(4, ne) * pi/90. )
         s5     = SIN( eps(5, ne) * pi/90. )
         s6     = SIN( eps(6, ne) * pi/90. )
         s7     = SIN( eps(7, ne) * pi/90. )
         x1     = s1 * et1s
         x2     = s2 * et2s
         x3     = s3 * et3s
         x4     = s4 * et4s
         x5     = s5 * et5s
         x6     = s6 * et6s
         x7     = s7 * et7s
         tax(1) = x1 * SIN(ds)
         tax(2) = x2 * SIN(fp)
         tax(3) = x3 * SIN(gd)
         tax(4) = x4 * SIN(hf)
         tax(5) = x5 * SIN(ig)
         tax(6) = x6 * SIN(ajh)
         tax(7) = x7 * SIN(aki)
         ta(1)  = x1 * COS(ds)
         ta(2)  = x2 * COS(fp)
         ta(3)  = x3 * COS(gd)
         ta(4)  = x4 * COS(hf)
         ta(5)  = x5 * COS(ig)
         ta(6)  = x6 * COS(ajh)
         ta(7)  = x7 * COS(aki)

c *** To make compatible with LPT convention for rewave, imwave

         do 545 ni = 1,7
            tawave(ni, ne) = 0.
            taxwve(ni, ne) = 0.
            tawave(ni, ne) = 0.5 * factor * ta(ni)
            taxwve(ni, ne) = 0.5 * factor * tax(ni)
 545     continue

         do 500 nwave = 1,nwaves
            etadum = eta(nwave, ne)
            ddum   = d(nwave, ne) * pi/180.
            IF (etadum .le. 0.) etadum = 1.
            rewave(nwave, ne) = 0.
            imwave(nwave, ne) = 0.
            IF (nwave .gt. nwaves) go to 500

c *** Nucleon case, mix partial waves, at least for MAW convention

            IF ((nwave .eq. 6) .OR. (nwave .eq. 7)) then
               etadum = etadum * COS( 2.0 * eps(2, ne) * pi/180. )
            endif
            IF ((nwave .eq. 9) .OR. (nwave .eq. 10)) then
               etadum = etadum * COS( 2.0 * eps(4, ne) * pi/180. )
            endif
            IF ((nwave .eq. 17) .OR. (nwave .eq. 18)) then
               etadum = etadum * COS( 2.0 * eps(1, ne) * pi/180. )
            endif
            IF ((nwave .eq. 20) .OR. (nwave .eq. 21)) then
               etadum = etadum * COS( 2.0 * eps(3, ne) * pi/180. )
            endif
            IF ((nwave .eq. 23) .OR. (nwave .eq. 27)) then
               etadum = etadum * COS( 2.0 * eps(5, ne) * pi/180. )
            endif
            IF ((nwave .eq. 12) .OR. (nwave .eq. 25)) then
               etadum = etadum * COS( 2.0 * eps(6, ne) * pi/180. )
            endif
            IF ((nwave .eq. 29) .OR. (nwave .eq. 31)) then
               etadum = etadum * COS( 2.0 * eps(7, ne) * pi/180. )
            endif

            rewave(nwave, ne) = 0.5 * etadum * SIN(2. * ddum) * factor
            imwave(nwave, ne) = 0.5 * (1. - etadum * COS(2. * ddum))
     $                              * factor
  500    continue

 1000 continue

c *** Store unfolded amplitudes on unit 3.

      rewind 3
      write(3) nes, nwaves
      write(3) (kap(ne),
     $         (rewave(nwave, ne), imwave(nwave, ne),
     $                            nwave = 1,nwaves), ne = 1,nes)
      write(3) (kap(ne), (tawave(ni, ne), taxwve(ni, ne), ni = 1,7),
     $                                                    ne = 1,nes)
      rewind 3
      ndata = 0
      write(6,919) nes, nwaves
c     write(6,940)

      if( nifty(20) .eq. 1 .OR. nifty(20) .eq. 3 ) then
         do 10 ne=1,nes
            write(6,941) kap(ne), Ecm(ne), plab(ne), Tlab(ne),
     $      (rewave(nwave, ne), imwave(nwave, ne), nwave = 1,nwaves),
     $      (tawave(ni, ne),    taxwve(ni,ne),        ni = 1,7)
 10      continue
      else
******************* tvm 8/26/92 phase shift printing c'ed out *******
c        do 11 ne = 1 , nes , nes-1
c           write(6,941) kap(ne), Ecm(ne), plab(ne), Tlab(ne),
c    $      (rewave(nwave, ne), imwave(nwave, ne), nwave = 1,nwaves),
c    $      (tawave(ni, ne),    taxwve(ni,ne),        ni = 1,7)
c11      continue
      endif

      garbage = 666.0
c     write(7,*) 'Booga, Booga, Booga!  Ooongowa!'
c-----Amplitudes all stored, interpolate on kappa, using 2 point lagrange
c *** Amplitudes all stored on tape 3, read back

 2000 continue
      ndata = 0
      read(3) nstore, nwaves
      write(6,923) nstore, nwaves

c *** Change do loop so UBC wont mess up the read

      read(3) (kap(ne), (rewave(nwave, ne), imwave(nwave, ne),
     $                   nwave = 1,nwaves), ne = 1,nstore)
      read(3) (kap(ne), (tawave(ni, ne),    taxwve(ni, ne),
     $                      ni = 1,7),      ne = 1,nstore)
      rewind 3
      IF (nes .le. nstore) go to 2001
      write(6,921) nstore
      nes = nstore
 2001 continue
      x = kappa
c     IF (kappa .gt. kap(nes)) x = kap(nes)
****************************************************************
      IF (kappa .lt. kap(1))   x = kap(1)


c *** Don't want scattering length for nucleon projectile.

c *** Here t and tx refer to interpolated values of rewave and imwave

c ***    tvm 2/24/92  replace all 100 's with 160's in next 4 lines ****
 2100 call LAGRNG (x, kap,   t, rewave, nes, nwaves, 2, 160, 40)
      call LAGRNG (x, kap,  tx, imwave, nes, nwaves, 2, 160, 40)
      call LAGRNG (x, kap,  ra, tawave, nes, 7, 2, 160, 7)
      call LAGRNG (x, kap, rax, taxwve, nes, 7, 2, 160, 7)

      t(13) = 0.
      tx(13) = 0.

CCC	For on-shell, the off-diagonal components are equal, therefore
CCC	set ra(i) = ra(i+7) etc. (pjd 4/4/93)
	do 2002 ip=1,7
	ra(ip+7) = ra(ip)
2002	rax(ip+7) = rax(ip)
CCC

c *** Special case H, I, J, and K NN-partial waves set = 0.

      IF (nifty(7) .eq. 1) go to 8006

      do 2424 nwave = 1,nwaves
         IF  ( nwave .eq. 10
     $               .OR.
     $         nwave .eq. 11
     $               .OR.
     $         nwave .eq. 12
     $               .OR.
     $         nwave .eq. 16
     $               .OR.
     $         nwave .ge. 24 )
     $   then
            t(nwave)  = 0.0
            tx(nwave) = 0.0
         endif
 2424 continue

      do 2525 nica = 4,7
         ra(nica)  = 0.0
         rax(nica) = 0.0
CCC	As above (+7)
	ra(nica+7) = 0.0
	rax(nica+7) = 0.0
 2525 continue

 8006 call MIXUPN (bn, bnf, t, tx)

      do 410 n = 1,32
         retij(n) = t(n)
         imtij(n) = tx(n)
 410  continue

      return

c-----formats
 902  format('0','# of energies, and # of deltas for each Ecm =',
     $        2i10,//)
 919  format(2i5)
 921  format(' ******nes > nstore=', i5,'nes set = store *****,',/)
 923  format('0', ' Tmatrices (TNCM)',
     $            ' calculated from Arndt SM83',
     $            ' phases nstore(nes) nwaves =', 2i5/)
 940  format('0','kap  Ecm  plab  Tlab  first line after rewave',
     $           ' imwave',/
     $           ' ', 13x, '1S0', 23x, '1D2', 23x, '1G4', 23x, '3P0',
     $           23x, '3P1',/
     $           ' ', 13x, '3P2', 23x, '3F2', 23x, '3F3', 23x, '3F4',
     $           23x, '3H4',/
     $           ' ', 13x, '3H5', 23x, '3H6', 23x, '1S0', 23x, '1P1',
     $           23x, '1F3',/
     $           ' ', 13x, '1H5', 23x, '3S1', 23x, '3D1', 23x, '3D2',
     $           23x, '3D3',/
     $           ' ', 13x, '3G3', 23x, '3G4', 23x, '3G5', 23x, '1I6',
     $           23x, '3J6',/
     $           ' ', 13x, '3J7', 23x, '3I5', 23x, '3I6', 23x, '3I7',
     $           23x, '1J7',/
     $           ' ', 13x, '3K7', 23x, '3K8')
 941  format('0', 4(f10.3),/, 100( 10e13.5,/))
 9001 format(8f10.6)
 9003 format(1x,f10.3/120(8f13.5/))
      end

