program overlap

c **********************************************************************
c *                                                                    *
c *                                                                    *
c *                                                                    *
c *                                                    Roger Yelle     *
c *                                                    21 Jan 1991     *
c *                                                                    *
c **********************************************************************

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

c       ---------------------------------------------------------
c            Arrays : IR band parameters
c       ---------------------------------------------------------
        dimension arad(12,3),pcol(12,3),dgn(12,3),rdmol(3),rmcol(3)
     &  , rmass(3), qrot(3), wave(12,3), wlortz(3)

c       ---------------------------------------------------------
c            Arrays : line listing
c       ---------------------------------------------------------
     &  , wav(nd4),str(nd4),wid(nd4),eng(nd4),nlines(4)

c       ---------------------------------------------------------
c          Arrays: absorption coefficients
c       ---------------------------------------------------------
     &  ,wt(100),tover(100),tline(100),tdopp(100)

c       ---------------------------------------------------------
c                           Data for Titan 
c       ---------------------------------------------------------
        data rmn2 / 28.d0 /
        data rdn2 / 3.d0 /
        data gm / 8.95d18 /
        data GMm / 4.16d-4 /
        data rplnt / 2575.d5 /
        data rh / 9.5d0 /
        data nu /3/
        data nlortz/300/
        data wlortz/5.2d-2,1.d-1,1.d-1/
        data nld/ 2.7d19/
        data tref/296.d0/

        write(unit=*,fmt=904)
        read (unit=*,fmt=*) tmp,prs,clm

c       ----------------------------------------------------------
c              Read in IR vibration transition parameters
c       ----------------------------------------------------------
        call bandread(wave,arad,pcol,dgn,rmass,rdmol,qrot)

c       ----------------------------------------------------------
c          Calculate reduced mass and collision diameters
c           for hydrocarbon-N2 collisions
c       ----------------------------------------------------------
        do im = 1,3 
           rdmol(im) = 1.d-8*half*(rdn2+rdmol(im))
           rmcol(im) = amu*rmn2*rmass(im)/(rmn2+rmass(im))
        end do

c       ----------------------------------------------------------
c             Read in line listing for fundamentals
c       ----------------------------------------------------------
        imol = 1
        call sread(1,wav,str,wid,eng,nlines,stot,wmean,wstart,wend)

        wstep=half*dsqrt(two*rkb*tmp/rmass(1)/amu)*wmean/cs

        write(unit=*,fmt=810) wstart,wend
 810    format(' wstart, wend = ',2f9.3)
        write(unit=*,fmt=811) 
        read (unit=*,fmt=*) wstart, wend, delw
 811    format(' Enter wstart, wend, delw > ',$)

        hckt = hck*(one/tmp - one/tref)
        wmin = wstart
        nb = 1
        do while(wmin .lt. wend)

           wmax = dmin1(wmin+delw,wend)
           wt(nb) = half*(wmin+wmax)
           wdop=dsqrt(two*rkb*tmp/rmass(1)/amu)*wt(nb)/cs
           lw1 = locate(wav,nlines,wmin)+1
           lw2 = locate(wav,nlines,wmax)

c          ------------------------------------------------------------
c                  Calculate transmission without overlap
c          ------------------------------------------------------------
           sumt = zero
           do lw = lw1, lw2 

              wlor = 1.d-6*prs*dsqrt(tmp/tref)*wid(lw)
              avgt = wlor/wdop
              wreach = dmax1(20.d0*wdop,nlortz*wlor)

              wlast = wav(lw)
              oplast = zero
              sfac = str(lw)*dexp(-hckt*eng(lw))/wdop/sqpi
              do while(wlast .lt. wav(lw)+wreach)

                 wnext = dmin1(wlast + wstep, wav(lw)+wreach)
                 wmid = half*(wnext+wlast)
                 wdel = wnext - wlast
                 xx = dabs(wmid-wav(lw))/wdop
                 call voigt(xx,avgt,vg)
                 opnext = vg*sfac
                 opmid = half*(opnext+oplast)
                 sumt=sumt+wdel*opmid*dexp(-opmid*clm)

                 wlast = wnext
                 oplast = opnext

              end do

           end do

           tline(nb) = nld*two*sumt/(wmax-wmin)

c          ------------------------------------------------------------
c                   Calculate transmission with overlap
c          -------------------------------------------------------------
           wlor = 1.d-6*prs*dsqrt(tmp/tref)*wlortz(1)
           wreach = dmax1(20.d0*wdop,nlortz*wlor)
           sumt = zero
           sumd = zero
           wlast = wmin 
           oplast = zero
           oplastk = zero
           do while(wlast .lt. wmax)

               wnext = dmin1(wlast + wstep, wmax)
               wmid = half*(wnext+wlast)
               wdel = wnext - wlast

               suml = zero
               sumk = zero
               do lw = lw1, lw2
                  dw = dabs(wmid-wav(lw))
                  sfac = str(lw)*dexp(-hckt*eng(lw))/wdop/sqpi
                  wlor = 1.d-6*prs*dsqrt(tmp/tref)*wid(lw)
                  avgt = wlor/wdop
                  if(dw .le. wreach) then
                     xw = dw / wdop
                     call voigt(xw,avgt,vg)
                     suml = suml + vg*sfac 
                     sumk = sumk + (dexp(-xw*xw))*sfac
                  end if
               end do
               opnext = suml
               opmid = half*(opnext+oplast)
               opnextk = sumk
               opmidk = half*(opnextk+oplastk)
               sumt=sumt+wdel*opmid*dexp(-opmid*clm)
               sumd=sumd+wdel*opmidk*dexp(-opmidk*clm)

               wlast = wnext
               oplast = opnext
               oplastk = opnextk

            end do

            tover(nb) = nld*sumt/(wmax-wmin)
            tdopp(nb) = nld*sumd/(wmax-wmin)

            write(unit=*,fmt=*) nb,tline(nb),tover(nb),tdopp(nb)

            wmin = wmax
            nb = nb + 1

        end do
        nbmax = nb - 1

c       ----------------------------------------------------------
c                        Write output to Disc
c       ----------------------------------------------------------
        write(unit=*,fmt=901)
        read (unit=*,fmt=900) fname
        open (unit=30, file=fname, status='new') 
           do nb = 1, nbmax
           write(unit=30,fmt='(1x,f10.4,2x,5e11.3)') 
     &      wt(nb),tover(nb),tline(nb),tdopp(nb)
     &      ,(tline(nb)-tover(nb))/tover(nb)
     &      ,(tdopp(nb)-tover(nb))/tover(nb)
           end do
        close(unit=30)

        stop

c       ----------------------------------------------------------
c                              Formats
c       ----------------------------------------------------------
  900   format(a)
  901   format(' Enter output file name > ',$)
  904   format(' Enter tmp, prs, clm > ',$)

        end