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