subroutine rybicki(tmin, dtau, u, st, ntau, xu, wu, nu
     & , opp, wpp, oppsum, npp, b1, b2, q)

c     --------------------------------------------------------
c          rybicki calculates the mean band intensity 
c     --------------------------------------------------------
      implicit double precision (o-z, a-h)
      logical itst
      include "parameters.dat"

c     --------------------------------------------------------
c                        Input Arrays
c     --------------------------------------------------------
      dimension dtau(1), u(nd1,1), st(nd1,1), xu(1), wu(1)
     & , opp(nd1, 1), wpp(1),  b1(1), b2(1), oppsum(1)

c     --------------------------------------------------------
c                        Output Arrays
c     --------------------------------------------------------
      dimension q(1)

c     --------------------------------------------------------
c                        Internal Arrays
c     --------------------------------------------------------
      dimension s(nd1), e(nd1, nd1), ta(nd1), tb(nd1), tc(nd1)
     & , w(nd1), tu(nd1, nd1), gf(nd1), dt(nd1)
      integer indx(nd1)

c     --------------------------------------------------------
c                      Initialize Arrays
c     --------------------------------------------------------
      do k = 1, ntau
         q(k) = zero
         do m = 1, ntau
            e(m,k) = zero
         end do
         e(k,k) = one
      end do

c     ---------------------------------------------------------
c                  Loop Through Frequencies 
c     ---------------------------------------------------------
      do j = 1, npp

            bf1 = b1(j)
            bf2 = b2(j)
            wgf = wpp(j)
            do k = 1, ntau
               gf(k) = opp(k,j)
            end do

c           -----------------------------------------------
c                  Monochromatic optical depth
c           -----------------------------------------------
            do k = 1, ntau - 1
               dt(k) = half*(gf(k)+gf(k+1))*dtau(k)
c               dt(k) = dmax1(half*(gf(k)+gf(k+1))*dtau(k),tmin)
            end do

c           -----------------------------------------------
c                 find optically thin level
c           -----------------------------------------------
            tmax = zero 
            do k = 1, ntau -1
               tmax = tmax + dt(k)
            end do

            k = 1
            tot = dt(k)
            do while((tot .lt. tmin).and.(k.lt.ntau))
               k = k + 1
               tot = tot + dt(k)
            end do
            ntmn = max(k,1)

            if ((tmax .lt. tmin).or.(ntmn .eq. ntau)) then

               do k = 1, ntau
                  q(k)=q(k)+half*wgf*gf(k)*(bf1+half*bf2)/oppsum(k)
               end do

            else

               do i = 1, nu

c                 -------------------------------------
c                    Calculate W matrix (Diagonal)
c                 -------------------------------------
                  do k = 1, ntau
                     w(k) = wu(i)*gf(k)*wgf/oppsum(k)
                  end do 

c                 -------------------------------------
c                    Calculate T matrix (Tridiagonal)
c                 -------------------------------------
                  do k = 1, ntmn-1
                     ta(k) = zero
                     tb(k) = one/tmin/tmin
                     tc(k) = -one/tmin/tmin
                  end do
                  call tmatrix(xu(i), dt(ntmn), ntau-ntmn+1
     &             , ta(ntmn), tb(ntmn), tc(ntmn))

c                 -------------------------------------
c                        LU decomposition of T
c                 -------------------------------------
c                  call ludcmpt(ta,tb,tc,ntau)

c                 -------------------------------------
c                     Calculate U matrix (Diagonal)
c                 -------------------------------------
                  do k = 1, ntau
                     do m = 1, ntau
                        tu(m,k) = zero
                     end do
                  end do

                  do k = ntmn, ntau
                     tu(k,k) = u(k,j)
                  end do

c                 -------------------------------------
c                    Calculate T(-1)*U store in TU
c                 -------------------------------------
c                  do k = 1, ntau
c                     call lubksbt(ta, tb, tc, tu(1,k), ntau)
c                  end do
                  do k = 1, ntau
                     call tridag(ta,tb,tc,tu(1,k),ntau,ierr)
                     if(ierr .eq. 1) then
                         write(unit=*,fmt=*) ' Error in rybicki'
                         do m = 1, ntau
                            write(unit=*,fmt=*) ta(m),tb(m),tc(m)
                         end do
                         stop
                     end if
                  end do

c                 ------------------------------------
c                   Calculate W*T(-1)*U add into E
c                 ------------------------------------
                  itst = .false.
                  do k = 1, ntau
                     aamax = zero
                     do m = 1, ntau
                        e(m,k) = e(m,k) - w(m) * tu(m,k)
                        if(dabs(e(m,k)) .gt. aamax) aamax = e(m,k)
                     end do
                     if(aamax .eq. zero) itst = .true.
                  end do

                  if(itst) then 
                  write(unit=*,fmt='(18h j,ntmn,tmax,tmin=,2i6,2e11.3)') 
     &             j,ntmn,tmax,tmin
                  do k = 1, ntau
                     write(unit=*,fmt='(i6,e11.3,3e16.8)') k,dt(k)
     &                ,ta(k),tb(k),tc(k)
                  end do
                  stop 100
                  end if
c                 -------------------------------------
c                    Calculate S vector (Rybicki's K)
c                 -------------------------------------
                  do k = 1, ntmn-1
                     s(k) = zero
                  end do
                  do k = ntmn, ntau-1
                     s(k) = st(k,j)
                  end do
                  s(ntau) = st(ntau,j) + (two*xu(i)/dt(ntau-1))
     &             *(bf1 + xu(i)*bf2)

c                 -------------------------------------
c                     Calculate T(-1)*S store in S
c                 -------------------------------------
c                  call lubksbt(ta, tb, tc, s, ntau)
                  call tridag(ta,tb,tc,s,ntau,ierr)
                  if(ierr .eq. 1) then
                      write(unit=*,fmt=*) ' Error in rybicki'
                      do k = 1, ntau
                         write(unit=*,fmt=*) ta(k),tb(k),tc(k)
                      end do
                      stop
                  end if

c                 -----------------------------------
c                    Calculate W*T(-1)*S add into Q
c                 ------------------------------------
                  do k = 1, ntau
                     q(k) = q(k) + w(k) * s(k)
                  end do

               end do 

            end if 

         end do
c     ----------------------------------------------------------
c                   Solve for mean line intensity
c     ----------------------------------------------------------
      call ludcmp(e, ntau, nd1, indx, dinx, ierr)
      if(ierr .eq. 1) then
         write(unit=*,fmt=*) ' Error in rybicki '
         do n = 1, ntau
            write(unit=*,fmt='(1x,10e11.3)') (e(i,j),j=1,ntau)
         end do
         stop 300
      end if
      call lubksb(e, ntau, nd1, indx, q)

      return 
      end

c ********************************************************************

      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
      dtp = dt(k)
      tb(k) = one + two * xu / dtp + two * u2 / dtp / dtp
      tc(k) = - two * u2 / dtp / dtp 

c     -------------------------------------
c                      Middle
c     -------------------------------------
      do k = 2, ntau - 1
        dtm = dtp
        dtp = dt(k)
        if(dtp .gt. tmin) then
           dta = half * (dtm + dtp)
           tfm = u2 / dta / dtm
           tfp = u2 / dta / dtp
           ta(k) = - tfm
           tb(k) = one + tfm + tfp
           tc(k) = - tfp
        else
           ta(k) = zero
           tb(k) = one/tmin/tmin
           tc(k) = -tb(k)
         end if
      end do

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

      return 
      end

c **********************************************************************

      subroutine ludcmpt(a, b, c, nmax)
      implicit double precision (o-z, a-h)
      dimension a(1), b(1), c(1)
      do n = nmax - 1, 1, -1
         if(dabs(b(n+1)) .lt. 1.d-30) then
            write(unit=*,fmt=*) n,' b = 0 in LUDCMPT '
            stop 111
         end if
         if(dabs(b(n+1)) .gt. 1.d+30) then
            write(unit=*,fmt=*) n+1,' b = inf in LUDCMPT '
            do m = nmax, 1, -1
               write(unit=*,fmt=*) m,a(m), b(m), c(m) 
            end do
            stop 112
         end if
         c(n) = c(n) / b(n + 1)
         b(n) = b(n) - c(n) * a(n + 1)
      end do
      return 
      end

      subroutine lubksbt(a, b, c, d, nmax)
      implicit double precision (o-z, a-h)
      dimension a(1), b(1), c(1), d(1), z(1000)
      z(nmax) = d(nmax)
      do n = nmax - 1, 1, -1
         z(n) = d(n) - (c(n) * z(n + 1))
      end do
      if(dabs(b(1)) .lt. 1.d-30) then
         write(unit=*,fmt=*) n,' b = 0 in LUBKSBT '
         stop 222 
      end if
      d(1) = z(1) / b(1)
      do n = 2, nmax, 1
         if(dabs(b(n)) .lt. 1.d-30) then
            write(unit=*,fmt=*) n,' b = 0 in LUBKSBT '
            stop 333 
         end if
         d(n) = (z(n) - (a(n) * d(n - 1))) / b(n)
      end do
      return 
      end

      subroutine tridag(a,b,c,d,n,ierr)
      implicit double precision (a-h,o-z)
      dimension a(1),b(1),c(1),d(1),g(1000)

      ierr = 0
      if(b(1) .eq. 0.d0) then 
        ierr = 1
        return
      end if
      bet = b(1)
      d(1) = d(1)/bet
      do j = 2, n
         g(j) = c(j-1)/bet
         bet = b(j)-a(j)*g(j)
         if(bet .eq. 0.d0) then
            ierr = 1
            return
         end if 
         d(j) = (d(j)-a(j)*d(j-1))/bet
      end do
       
      do j = n-1, 1, -1
         d(j) = d(j)-g(j+1)*d(j+1)
      end do

      return
      end