! To extract a HCN/HNC  linelist from the Einstein A, coeffients 
! of Harris et al. 2002, ApJ, 578,657 and energy levels of
! Harris et al. MNRAS (Submitted).
!
! INPUT FILES:
!
! fort.99              An Einstein A coefficent file from Harris et al. 2002.
! enlvs-combined.txt   The energy level file from Harris et al. 2005.
!
! OUTPUT FILE:
!
! linelist.txt         The extracted linelist of the format outlined in
!                      Harris et al. 2005.
!
! The user my wish to chnage the following variables
! 
! T:     The temperature for which absolute intensity and the partition function
!        are computed
! simin: The minimum intensity (cm mol^{-1}) to be output. 
! numin: Minimum frequncy output (cm^{-1})
! numax: Maximum frequncy output (cm^{-1})
! zai  : if true output an ab initio linelist, if false output the linelist
!        which uses lab determined energy levels where possible.

! iout : Selects the format of the output linelist.
!        1:  output the minimum data necisarry for an opacity calculation
!            with Einstein A coeffients (nu, E", J", J', A).
!        2:  outputs linelist in format given in Harris et al. (2005), plus the
!            integrated absorption intesity at T (cm molecule^{-1}).
!        3:  outputs data in compressed Einstien A coeffient file 
!            (index1, index2, A) 
!        4:  outputs data in compressed Einstien A coeffient file, with frequency
!            (index1, index2, A, nu)
!        5:  Outputs linelist in format given in Harris et al. (2005), plus the 
!            dimensionless weighted oscillator strength (gf).
!        6:  output the minimum data necisarry for an opacity calculation
!            with weighted oscillatior strengths (nu, E", gf).
!        ANY OTHER VALUE: Outputs linelist in format given in Harris et al. (2005).


! zmin : if true output the minimum data necisarry for an opacity calculation
!        (nu, E", J", J', A).
! zSi  : if .true. and zmin is false, this outputs the integrated absorption
!        intesity at T (cm molecule^{-1}) as well as the Einstein A coefficent.
!
! Author G. J. Harris     August 2005.
      implicit none
      integer in, nen,out
      parameter (in=10,out=11,nen=168110)
! storeage space for energy level data.
      real(8), allocatable, dimension(:) :: eai,ex
      integer, allocatable, dimension(:) :: JJ,Pp,nn,v1,v2,ll,v3,iso
! data
      real(8) const,k,gfk
! user changable variables
      real(8) T, simin, numin,numax
      integer iout
      logical zai 
! local variables
      real(8) A, Si,nu, El,eu,e1,e2,beta,expfact,Q,intS,gf
      integer ind1,ind2,iu,il,count
      character(2) at

      data const /1.3270914d-12/
      data gfk /1.499197d0/
      data k /0.6950386d0/

!----------------------------
! user changeable variables
      T=3000.0d0
      simin=1.0d-26
      numin=2000.0d0
      numax=3800.0d0
      zai=.false.
      iout=5
! ---------------------------

      beta=1.0d0/(k*T)

      allocate(eai(nen),ex(nen),JJ(nen),Pp(nen),iso(nen))
      allocate(nn(nen),v1(nen),v2(nen),ll(nen),v3(nen))

! read energy levels
      call read_energy(eai,ex,JJ,Pp,nn,iso,v1,v2,ll,v3)
! compute partition function
      call pfcalc(eai,ex,jj,T,Q)
      write(6,*) "T = ", T, "Q = ", Q
! read transition strength file, line by line.
      open(unit=in,form='formatted',file='fort.99')
      open(unit=out,form='formatted',file='linelist.txt')
      ints=0.0d0
      count=0

10    read(in,1001,end=92) ind2, ind1, A

       if(ex(ind1) .ge. 0.0d0 .and. ex(ind2) .ge. 0.0d0 .and. &
          .not. zai) then
        nu=ex(ind2)-ex(ind1)
        e1=ex(ind1)
        e2=ex(ind2)
        at='lb'
       else
        nu=eai(ind2)-eai(ind1)
        e1=eai(ind1)
        e2=eai(ind2)
        at='ai'
       endif
       if (nu .eq. 0.0d0) goto 10

       if(nu .lt. 0.0d0) then
        il=ind2
        iu=ind1
        el=e2
        eu=e1
        nu=-nu
       else
        il=ind1
        iu=ind2
        el=e1
        eu=e2
       endif

       if(nu .lt. numin .or. nu .gt. numax) goto 10

       expfact=exp(-El*beta)*(1.0d0-exp(-nu*beta))
       SI=(const/(Q*(nu**2)))*expfact*A*dble((2*JJ(iu))+1)

       if(Si .ge. simin) then
        count=count+1
        ints=ints+si
        if(iout .eq. 1) then
         write(out,1004) nu,el,JJ(il),JJ(iu),A
        elseif(iout .eq. 2) then
         write(out,1002) nu,el,JJ(il),pp(il),nn(il),JJ(iu),pp(iu),&
                     nn(iu),A,Si,il,iu,&
                     iso(il),v1(il),v2(IL),ll(il),v3(il),&
                     iso(iu),v1(iu),v2(Iu),ll(iu),v3(iu), at
        elseif(iout .eq. 3) then
         write(out,1001) iu, il, A
        elseif(iout .eq. 4) then
         write(out,1005) iu, il, A, nu
        elseif(iout .eq. 5) then
         gf=(gfk*dble((2*JJ(iu))+1)*A) / (nu*nu)
         write(out,1002) nu,el,JJ(il),pp(il),nn(il),JJ(iu),pp(iu),&
                     nn(iu),A,gf,il,iu,&
                     iso(il),v1(il),v2(IL),ll(il),v3(il),&
                     iso(iu),v1(iu),v2(Iu),ll(iu),v3(iu), at
        elseif(iout .eq. 6) then
         gf=(gfk*dble((2*JJ(iu))+1)*A) / (nu*nu)
         write(out,1006) nu,el,gf
        else
         write(out,1003) nu,el,JJ(il),pp(il),nn(il),JJ(iu),&
                     pp(iu),nn(iu),A,il,iu,&
                     iso(il),v1(il),v2(IL),ll(il),v3(il),&
                     iso(iu),v1(iu),v2(Iu),ll(iu),v3(iu), at
        endif
       endif

      goto 10
92    close(in)
      close(out)
      write(6,*) "Integrated intensity ", ints
      write(6,*) "number of lines output ", count

      deallocate(eai,ex,JJ,Pp,nn,iso,v1,v2,ll,v3)
 1001 format(2i7,1p,d14.6)
 1002 format(f12.6,f13.6,i3,i2,i5,i3,i2,i5,1p,2e10.3,0p,2i7,1x,5i3,1x,5i3,1x,a2)
 1003 format(f12.6,f13.6,i3,i2,i5,i3,i2,i5,1p,e10.3,0p,2i7,1x,5i3,1x,5i3,1x,a2)
 1004 format(f12.6,f13.6,2i3,1p,e10.3)
 1005 format(2i7,1p,d14.6,0p,f13.6)
 1006 format(f12.6,f13.6,e10.3)
     end

! -----------------------------------------------------------------

      subroutine read_energy(eai,ex,JJ,Pp,nn,iso,v1,v2,ll,v3)
      ! the primary function of this subroutine is to read and store the 
      ! energy level data.
      implicit none
      integer in, nen
      parameter (in=10,nen=168110)
! arguments
      real(8) Eai(*), Ex(*)
      integer JJ(*),Pp(*),nn(*),iso(*),v1(*),v2(*),ll(*),v3(*)
! local variables
      integer it, jt,pt,nt,isot,v1t,v2t,v3t,lt
      real(8) eait,ext
      integer i

      open(unit=in,form='formatted',file='enlvs-combined.txt')
      do 10 i=1, nen
       read(in,1001) it,Jt,Pt,nt,eait,isot,v1t,v2t,lt,v3t,ext
!       write(6,1001) it,Jt,Pt,nt,eait,isot,v1t,v2t,lt,v3t,ext
       JJ(it)=Jt
       nn(it)=nt
       pp(it)=pt
       iso(it)=isot
       v1(it)=v1t
       v2(it)=v2t
       v3(it)=v3t
       ll(it)=lt
       eai(it)=eait
       ex(it)=ext
       if(it .gt. 1 .and. ex(it) .lt. 1.0d-10) ex(it)=-1.0d0
!       write(6,1001) it,Jt,Pt,nt,eait,isot,v1t,v2t,lt,v3t,ex(it)
10    continue
      close(in)

      

      return
1001  format(i6,i3,i2,i5,f13.6,5i3,f13.6,1p,e11.2,0p,2x,a1)
      end

! ---------------------------------------------

     subroutine pfcalc(eai,ex,jj,T,Q)
     implicit none
     integer nen
     parameter (nen=168110)
! arguments
     real(8), dimension(*) :: eai,ex
     integer, dimension(*) :: jj
     real(8) T,Q
! data
     real(8) k
! local variables
     real(8) beta, entmp
     integer i
     data k /0.6950333d0/

     beta=1.0d0/(k*T)
     Q=0.0d0

     do 10 i=1, nen

      if(ex(i) .lt. 0.0d0) then
       entmp=eai(i)
      else
       entmp=ex(i)
      endif

      Q=Q+(dble((2*jj(i))+1)*exp(-entmp*beta))

10   continue

     return
     end

! ---------------------------------------------------







