program hcnheat 

        implicit double precision (o-z, a-h)
        include "parameters.dat"
        character fname*64

        dimension wav(200),str(200),wid(200),eng(200),tmp(28),ht(28)
c       ----------------------------------------------------------
c             Read in line listing for fundamentals
c       ----------------------------------------------------------
        ib=9
        call sread(ib, wav, str, wid, eng, nlines,stot,wav0
     &       ,wmin,wmax)

c        open(unit=20, file='/home/asgard/yelle/lines/hcna.dat')
c           read(unit=20, fmt=*) nlines
c           stot = 0.d0
c           do n = 1, nlines
c              read(unit=20, fmt='(f10.6,e10.3,f5.3,f10.3)') 
c     &         wav(n), str(n), wid(n), eng(n)
c              stot = stot + str(n)
c           end do
c        close(unit=20) 

c          -----------------------------------------------------
c                     Loop through lines (debug )
c          -----------------------------------------------------
           do lw = 1, nlines
              write(unit=*,fmt='(i5,f10.6,1x,e10.3,1x,f11.6)') 
     &             lw,wav(lw),str(lw),eng(lw)
           end do

c       --------------------------------------------------------
c                     Loop through temperatures
c       --------------------------------------------------------
        ndat = 28
        do n = 1, ndat
           tmp(n) = 30.d0 + 10.*(n-1) 
           write (*,*) n, tmp(n)
           hckt = hck*(one/tmp(n)-one/296.d0)
           qrot = 296.d0/tmp(n)

c          -----------------------------------------------------
c                     Loop through lines
c          -----------------------------------------------------
           sum = zero
           do lw = 1, nlines

c              write(unit=*,fmt='(i5,f10.6,1x,e10.3,1x,f11.6)') 
c     &         lw,wav(lw),str(lw),eng(lw)

              hckw = hck*wav(lw)
              hckwt0 = hckw/296.d0
              hckwt = hckw/tmp(n)
              strn = qrot*str(lw)*dexp(-hckt*eng(lw))
     &         *(one-dexp(-hckwt))/(one-dexp(-hckwt0))
              bp = bplanck(wav(lw),tmp(n))
              sum = sum + bp*strn

           end do
           ht(n) = four*pi*sum
           write (*,*) n, tmp(n), ht(n)
        end do

c       ----------------------------------------------------------
c                        Write output to Disc
c       ----------------------------------------------------------
        write(unit=*, fmt=905) 
        read (unit=*, fmt=900) fname
        open (unit=30, file=fname, status='new') 

           do n = 1, ndat
              write(unit=30, fmt=906) tmp(n),ht(n)
           end do

        close(unit=30)

        stop

c       ----------------------------------------------------------
c                              Formats
c       ----------------------------------------------------------
  900   format(a)
  905   format(' Enter filename for output > ',$)
  906   format(1x,0pf7.1,1x,1p9e11.3)

        end


        function bplanck(wav, temp)
        implicit double precision (o-z, a-h)
        parameter (one = 1.d0, c1 = 1.192d-5, c2 = 1.439565d0)
        c3 = dexp(- ((c2 * wav) / temp))
        bplanck = ((c1 * (wav ** 3)) * c3) / (one - c3)
        return 
        end