c amatrixexport.f
c f77 amatrixexport.f -o amatrixexport
c
c see appendix C of Martin, Schwarz, and Mandy 1996, ApJ, 461, 265.
c
      common /af/ a(349,349)
      common /energies/ ecol(348)
c
      write(6,*) 
     $     ' unified export program containing all functionality'
      write(6,*) 
     $     ' referred to in appendix C of '
      write(6,*) 
     $     ' Martin, Schwarz, and Mandy 1996, ApJ, 461, 265'
      write(6,*) 
     $     ' (for collision rates see'
      write(6,*) 
     $     ' Martin and Mandy 1995, ApJ, 455, L89)'
      write(6,*) 
     $     '    reads in A values, tunneling rates, energies'
      write(6,*) 
     $     '    includes an interactive test'
      write(6,*) ' '
      write(6,*) 
     $     ' required datasets: amatrix.dat, tunnel.dat, ecol.dat'
      write(6,*) ' '
c
c read in the A values in s-1
c
      write(6,*) ' reading in the A values from amatrix.dat'
      write(6,*) ' '
      call amatrxread
c
c read in the energies for check and exercise below
c
      write(6,*) ' reading in the energies from ecol.dat'
      write(6,*) 
     $     '    the energies are in wavenumbers relative to'
      write(6,*) 
     $     '    the energy of two separated H atoms'
      write(6,*) ' '
      call enerin
c
      write(6,*) ' basic check of the A values for invalid entries'
      write(6,*) '    exercises the routine colvj which converts' 
      write(6,*) '    the column/row index to v and j'
      call amatrxcheck
      write(6,*) ' end of amatrxcheck'
      write(6,*) ' '
c
c read in the tunneling rates in s-1
c
      write(6,*) ' reading in the tunneling rates from tunnel.dat'
      write(6,*) ' '
      call amatrxtunn
c
c to see all the non-zero tunneling elements, run check program again
c      call amatrxcheck
c
c create the diagonal element
c
      write(6,*) ' creating the diagonal element'
      write(6,*) ' '
      call amatrxdiag
c
c now test the routine
c
      write(6,*) 
     $     ' interactive testing routine'
      write(6,*) 
     $     '    exercises the routine vjcol which converts'
      write(6,*) 
     $     '    input v and j to column/row index'
      write(6,*) 
     $     '    for the vector ecol and matrix a'
      write(6,*) 
     $     '    use 99,99 for the dissociated state'
      write(6,*) 
     $     '    to end, put in negative quantum number v'
      write(6,*) ' '
      write(6,*) 
     $     ' radiative transition is named; wavelength is in microns'
      write(6,*) ' tunneling rates are given'
      write(6,*) 
     $     ' for diagonal elements, the total exit rate is reported'
      write(6,*) ' '
      write(6,*) ' '
c
 1    continue
      write(6,*) ' initial state'
      read(5,*) iv0,ij0
      if(iv0.lt.0) go to 2
      call vjcol(iv0,ij0,icol0)
      if(icol0.lt.0) go to 1
      write(6,*) ' final state'
      read(5,*) ivf,ijf
      call vjcol(ivf,ijf,icolf)
      if(icolf.lt.0) go to 1
c
c handle those involving dissociated state separately
      if(icolf.eq.349) then
         write(6,*) ' final state is dissociation '
         if(icol0.eq.349) then
            write(6,*)  ' diagonal element'
            write(6,3) a(icolf,icol0)
        else
            write(6,*)  ' tunneling rate'
            write(6,3) a(icolf,icol0)
         endif
         go to 1
      endif
c
c other than the diagonal to check, these are zero
      if(icol0.eq.349) then
         write(6,*) ' initial state is dissociated '
         write(6,*)  
     $        ' no valid transition'
         write(6,3) a(icolf,icol0)
         go to 1
      endif
c         
      idi=ijf-ij0
      idif=iabs(idi)
      edif=ecol(icol0)-ecol(icolf)
c wavelength in microns
      if( (idif.eq.0 .or. idif.eq.2) .and. iv0.ge.ivf) then
         if(edif.le.0.) then
            if(icolf.ne.icol0) then
               write(6,*) ' energy inverted: no transition'
            else
               write(6,*) ' diagonal element: total exit'
            endif
         else
            if((ij0+ijf).eq.0) then
               write(6,*) ' no transition for j=jp=0'
            else
               alam=10000./(edif+1.e-35)
               if(idi.eq.-2) write(6,10) iv0,ivf,ijf,alam
               if(idi.eq.0) write(6,11) iv0,ivf,ijf,alam
               if(idi.eq.2) write(6,12) iv0,ivf,ijf,alam
 10            format(1x,i2,'-',i2,' S(',i2,')',f10.3)
 11            format(1x,i2,'-',i2,' Q(',i2,')',f10.3)
 12            format(1x,i2,'-',i2,' O(',i2,')',f10.3)
            endif
         endif
      else
         write(6,*) ' no transition'
      endif
      write(6,3) a(icolf,icol0)
 3    format(1x,1p8e9.2)
      go to 1
c
 2    continue
      write(6,*) ' goodbye'
c
c end of test
c
c
      stop
      end
c
      subroutine amatrxread
      common /af/ a(349,349)
      open(unit=1,status='old',file='amatrix.dat')
      do 400 i=1,349
         read(1,401) (a(i,k),k=1,349)
 400  continue
 401  format(1p8e9.2)
c
      close(unit=1)
      return
      end
c
      subroutine amatrxcheck
      common /af/ a(349,349)
      common /energies/ ecol(348)
      write(6,*) ' beginning amatrxcheck'
      iquasi=0
c      do 1 icol0=1,348
      do 1 icol0=1,349
         call colvj(icol0,iv0,ij0)
c         do 2 icolf=1,348
         do 2 icolf=1,349
            iflag=0
            call colvj(icolf,ivf,ijf)
            idi=ijf-ij0
            idif=iabs(idi)
            if(icol0.eq.349 .or. icolf.eq.349) then
               edif=1.
            else
               edif=ecol(icol0)-ecol(icolf)
            endif
            if(icol0.eq.349 .and. icolf.eq.349) edif=0.
c wavelength in microns
            alam=10000./(edif+1.e-35)
            if(idif.eq.0 .or. idif.eq.2 .and. iv0.ge.ivf) then
c possible transition
               if(edif.le.0.) then
                  if(a(icolf,icol0).ne.0.) then
                     iflag=1
                     if(icolf.ne.icol0) then
                        write(6,*) ' energy inverted: no transition'
                     endif
                  endif
               else
c allowed by energy difference
                  if((ij0+ijf).eq.0) then
                     if(a(icolf,icol0).ne.0.) then
                        iflag=1
                        write(6,*) ' no transition for j=jp=0'
                     endif
                  endif
c allowed transition
                  iflag=0
                  if(ecol(icol0).gt.0.) then
                     if(a(icolf,icol0).ne.0.) then
c                        iflag=1
c                        write(6,*) ' quasi bound non zero'
c put counter in to avoid diagnostic output; these are valid additions
c for states near bound states.
                        iquasi=iquasi+1
                     endif
                  endif
               endif
            else
               if(a(icolf,icol0).ne.0.) then
                  iflag=1
                  write(6,*) ' no transition'
               endif
            endif
            if(icolf.ne.icol0) then
               if(iflag.eq.1) write(6,*) icolf,icol0,iv0,ij0,ivf,ijf
               if(iflag.eq.1) write(6,3) a(icolf,icol0)
            endif
 3          format(1x,1p8e9.2)
 2       continue
 1    continue
      if(iquasi.ne.38) write(6,*) 
     $     ' incorrect number of A values for quasibound'
c
c end of check
c
      return
      end
c
      subroutine amatrxtunn
      common /af/ a(349,349)
      open(unit=2,status='old',file='tunnel.dat')
      read(2,402) (a(349,k),k=1,349)
 402  format(6(1x,1pd12.4))
c
      close(unit=2)
      return
      end
c
      subroutine amatrxdiag
      common /af/ a(349,349)
c
c generate diagonal element which describes depletion of a state j to
c anywhere (all i): best to do this on reading in for accuracy
c
      do 1400 j=1,349
c diagonal elements are all zero when read in
         sum=0.0
c sum along columns
         do 1402 i=1,349
            sum=sum-a(i,j)
 1402    continue
         a(j,j)=sum
 1400 continue
c
      return
      end
c
      subroutine vjcol(n,j,icol)
c  routine using general sindex to find the column index for populations
c
c  input
c  n      state vibrational quantum number 0, 14
c  j      state rotational quantum number 0, 38 or 1, 37
c  n=j=99, ind=1 for dissociation (used in other programs)
c
c  output
c  icol:  column index  1, 349
c     para:             1, 177
c     ortho:            178, 348 (171)
c     dissociation:     349
c
      ioff=177
      if(n.lt.99) then
         call sindex(n,j,icol)
         if(icol.lt.0) go to 1
         if(mod(j,2).eq.1) icol=icol+ioff
      else
         icol=349
      end if
 1    continue
      return
      end
c
      subroutine sindex(n,j,indvj)
c  this subroutine returns indvj [1-178] given 
c   n     vibrational quantum number
c   j     rotational quantum number
      dimension iv(39)
c      data iv/15,15,15,15,15,14,14,14,14,13,13,13,13,12,12,11,11,10,10,
      data iv/15,15,15,15,15,14,14,14,14,13,13,13,12,12,12,11,11,10,10,
     19,9,9,8,8,7,7,6,6,5,5,4,4,3,3,2,2,1,1,1/
      if( n.gt.14 .or. n.lt.0) go to 200
      if( j.gt.38 .or. j.lt.0) go to 200
      if( (n+1).gt.iv(j+1) ) go to 200
      istart=1
      if(mod(j,2).eq.1)istart=2
      icount=1
      istop=j
c avoid the uncertain convention of a single pass or not on negative first step
      if(istart.gt.istop)go to 20
      do 10 i=istart,istop,2
      icount=icount+iv(i)
   10 continue
   20 continue
      icount=icount+n
      indvj=icount
      return
 200  continue
      write(6,9998) n,j
 9998 format(1x,'n, j = ',2i5,' illegal v,j combination')
      indvj=-1
      return
      end
c
      subroutine colvj(icol,n,j)
C  this subroutine returns quantum numbers given the column index
c
c  input
c  icol:  column index  1, 349
c     para:             1, 177
c     ortho:            178, 348 (171)
c     dissociation:     349
c  output
C  n      state vibrational quantum number 0, 14
C  j      state rotational quantum number 0, 38 or 1, 37
c  n=j=99, ind=1 for dissociation
c
      dimension iv(39)
c  this array is the maximum of v+1 for j=0,38
c      data iv/15,15,15,15,15,14,14,14,14,13,13,13,13,12,12,11,11,10,10,
      data iv/15,15,15,15,15,14,14,14,14,13,13,13,12,12,12,11,11,10,10,
     19,9,9,8,8,7,7,6,6,5,5,4,4,3,3,2,2,1,1,1/
      idiss=349
      ioff=177
c  Compute v and j
      if(icol.ne.idiss) then
         if(icol.gt.ioff) then
c fortran indices
c ortho
            istart=2
            icheck=icol-ioff
         else
            icheck=icol
            istart=1
         endif
         do 100 i=istart,39,2
            icheck=icheck-iv(i)
            if(icheck.le.0)go to 101
 100     continue
         write(6,9999)icol
 9999    format(1x,'icol=',i3,' How did it get here?')
         stop
 101     continue
         j=i-1
         n=icheck+iv(i)-1
      else
c dissociation
         n=99
         j=99
      end if
      return
      end
c
      subroutine enerin
      common /energies/ ecol(348)
      open(unit=3,status='old',file='ecol.dat')
c energies in wavenumbers relative to dissociated H
            read(3,501) (ecol(k),k=1,348)
 501  format(8f10.2)
      close(unit=3)
      return
      end
c
