subroutine tmatrix(xu, dt, ntau, ta, tb, tc)

c     ------------------------------------------------
c          TMATRIX calculates the tridaigonal transfer
c          matrix for the Feautier form of the R-T eq.
c          Second order boundary conditions from
c          Auer(1967) Ap J,vol 150,p l53
c     ------------------------------------------------

      implicit double precision (o-z, a-h)
      parameter (zero = 0.d0, half = 0.5d0, one = 1.d0, two = 2.d0)
      dimension dt(1),  ta(1), tb(1), tc(1)

c     -------------------------------------
c                  Upper Boundary
c     -------------------------------------
      u2 = xu * xu
      k = 1
      tb(k) = one + two * xu / dt(1) + two * u2 / dt(1) / dt(1)
      tc(k) = - two * u2 / dt(1) / dt(1) 

c     -------------------------------------
c                      Middle
c     -------------------------------------
      do k = 2, ntau - 1
        dta = half * (dt(k-1) + dt(k))
        tfm = u2 / dta / dt(k-1)
        tfp = u2 / dta / dt(k)
        ta(k) = - tfm
        tb(k) = one + tfm + tfp
        tc(k) = - tfp
      end do

c     -------------------------------------
c                  Lower Boundary
c     -------------------------------------
      k = ntau
      dtm = dtp
      ta(k) = - two * u2 / dt(k-1) / dt(k-1)
      tb(k) = one + two * xu / dt(k-1) + two * u2/ dt(k-1) / dt(k-1)

      return 
      end