program fluo

c----------------------------------------------------------------------
c  version 002    2 May, 1999
c----------------------------------------------------------------------
c            written by      Daniel Haskel 
c                            Department of Physics
c                            Box 351560
c                            University of Washington
c                            Seattle, WA   USA 98195
c                      
c                 phone      (206) 543-0435 
c                 e-mail     haskel@phys.washington.edu
c----------------------------------------------------------------------
c      description of the code:
c 
c  Fluo normalizes, aligns, and corrects xmu scans for self
c  absorption effects in fluorescence data. The normalization is done
c  by subtracting a pre-edge line then extrapolating a post-edge line to
c  e0. The e0 intercept is used for edge-step normalization and alignment.  
c  The first derivative is also used in the alignment procedure. The
c  self absorption correction is evaluated using tabuted cross sections in
c  McMaster tables.
c----------------------------------------------------------------------
c       sample input file:
c 
c   e0 = 37442
c   central = ba   edge=k
c   atmwei = true
c   angin=60    angout=30
c   formin=uwexafs formout=uwexafs      
c   out =bafile
c   atoms
c	ba    1
c	fe    12
c	o     19
c    files
c         bafeo_amor.xmu 1   * ba edge amorphous 36 K
c         bafeo_amor.xmu 3   * ba edge * amorphous 300 K
c  !	  febao_crys300.dat  * fe k-edge crystalline sample

c----------------------------------------------------------------------
c       version history
c     001      pre-pre-pre version. No alignment. 
c     002      Mcfluo incorporated to calculate Xsections mutf, mub, muc
c              from atom list with relative amounts, central atom + edge.
c              Energy alignment of several files incorporated.
c     5-21-99  If fluo-corrections calculated, added another edge-step
c     normalization after correction is done
c----------------------------------------------------------------------

      parameter(nptx=2**11, nfilx=50, ndocx=25, iou=10,
     $            nwdx=5)
      parameter(epsi=1.e-5, epsi2=2.e-3, zero=0.e0, pi= 3.141592654 )
      character*5 vnum
      parameter(vnum='0.02')
      
      dimension     energy(nptx), xmuraw(nptx), xmuout(nptx)
      dimension     foo1(nptx), foo2(nptx), foo3(nptx)
      dimension     nkey(nfilx)
      dimension     der(nptx)
      logical       eefind, stfind, vaxflg, lfluo, langle, le0fin
      character*3   cnum, test
      character*10  filtyp, frmin, frmout, skey(nfilx), skeyo, tmpe0
      character*78  outfil, file(nfilx), line, of, string,
     $              commnt(nfilx), out2
      character*100 doc(ndocx)
      real          mutf, mub, muc
      real          maxder, xmuint
      external      xmuint

      common /idata/ nxmu
      save /idata/
      common /rdata/ e0, energy, xmuout
      save /rdata/

      character*78 comm

 400  format(a)
 410  format(1x,f13.6,3x,f13.6)
 420  format(75('='))
 430  format(i3)
 432  format(i3.3)
 435  format(f9.2)
 440  format('Fluo --',a5,'Correcting for self absorption.'16x,
     $     'by Daniel Haskel')

      vaxflg   = .false.
c VAX USERS:  change the line above to .true.
c VAX VAX VAX VAX VAX VAX VAX VAX VAX VAX VAX VAX VAX VAX VAX VAX VAX 
      test = 'abc'

      write(line,420)
      call messag(line)
      write(string,440)vnum
      call messag(string)
      call messag(line)

      call noinit(nkey, predg1, predg2, enor1, enor2,
     $            foo1, foo2, foo3,
     $            of, file, skey, commnt, angin, angout,
     $            mutf, mub, muc)

      call normin(nfile, nkey, file, e0, 
     $            enor1, enor2, predg1, predg2,
     $            frmin, frmout, of, skey, commnt, 
     $            mutf, mub, muc, angin, angout,
     $            le0fin, lfluo, langle)


      do 50 i=1,nfile

c       --- read the data
        filtyp = ' '
        nxmu   = nptx
        ndoc   = ndocx
        call inpdat(filtyp, frmin, file(i), vaxflg, skey(i), nkey(i), 
     $       ndoc, doc, nxmu, energy, xmuraw, foo1, foo2, foo3 )
        call case(test,filtyp)
        call triml(frmin)
        if ((frmin(1:2).eq.'uw').and.(filtyp.ne.'xmu')) then
            call messag(' *** ERROR:  uwexafs input file must be '//
     $                  'of type XMU.')
            call messag('             Skipping data set.')
            goto 50
        endif


        if (frmout.eq.' ')   frmout = frmin



c       --- run time message about input data file

        if (lfluo) then
            comm = ' Normalizing and correcting data from "'
        elseif (i.eq.1) then
            comm = ' Normalizing data from  "'
        elseif (i.ne.1) then
            comm = ' Normalizing and aligning data from "'
        endif    

        ii = istrln(file(i))
        ic = istrln(comm)
        if (frmin(1:2).eq.'uw') then
            write(cnum,430)nkey(i)
            call messag(comm(:ic)//file(i)(:ii)//'", nkey = '//cnum)
        else
            call messag(comm(:ic)//file(i)(:ii)//'"')
        endif

c       --- remove pre-edge and normalize.  use step as initial guess
c       --- in all other data
        step   = 0
        stfind = .true.
        slopre = 0
        bpre   = 0


        if (i.eq.1) then
             if ((le0fin).and.(e0.gt.zero)) then
              eefind=.false.
             else
              eefind=.true.
             endif
   
             call prenew(eefind, e0, predg1, predg2, enor1, enor2,
     $              nxmu, energy, xmuraw, step, stfind, slopre,
     $              bpre, xmuout)
                if (eefind) then
                   write(tmpe0, 435) e0
ccc                   write(tmpe0, '(f7.2)') e0
                  call messag(' e0 not provided; Found by fluo to be'//
     $                         ' e0 = '//tmpe0)
                endif   
                e0ref = e0
c                print*,'e0ref=', e0ref
                if (step.lt.epsi) step = 1.e0
                do 20 j=1,nxmu
                   xmuout(j) = xmuout(j) / step
 20             continue 
                h=1.e0
                do 21 j=1,nxmu
                   der(j)= dfridr(xmuint,energy(j),h,err)
 21             continue
                maxder=0.e0
                do 22 j=1,nxmu
                   maxder=max(der(j),maxder)
 22             continue
                inmax=1
                do 23 j=1,nxmu
                   if (abs((der(j)-maxder)).lt.epsi) inmax=j  
 23             continue
                do 24 j=inmax,nxmu
c                   if (abs((der(j)-zero)).lt.epsi2) then
                   if (der(j).lt.zero) then
                      inzero=j
                      goto 25
                   endif
 24             continue
 25             continue
c                print*,'inmax=',inmax
c                print*,'inzero=',inzero
                derref=energy(inzero)
                emader=energy(inmax)
c                print*,'derref=',derref
c                print*,'emader=',emader
        else
            eefind = .true.
            call prenew(eefind, e0, predg1, predg2, enor1, enor2,
     $              nxmu, energy, xmuraw, step, stfind, slopre,
     $              bpre, xmuout)
c                print*, e0
                e0shft = e0 - e0ref
c                print*,'e0shft=', e0shft
               if (step.lt.epsi) step = 1.e0
                do 30 mm=1,nxmu
                   xmuout(mm) = xmuout(mm) / step
                   energy(mm) = energy (mm) - e0shft
 30             continue
                h=1.e0
                do 31 mm=1,nxmu
                   der(mm)=dfridr(xmuint,energy(mm),h,err)
 31             continue
                maxder=0.e0
                do 32 mm=1,nxmu
                   maxder=max(der(mm),maxder)
 32             continue
                inmax=1
                do 33 mm=1,nxmu
                   if (abs((der(mm)-maxder)).lt.epsi) inmax=mm
 33             continue
                do 34 mm=inmax,nxmu
c                   if (abs((der(mm)-zero)).lt.epsi2) then
                   if (der(mm).lt.zero) then
                      inzero=mm
                      goto 35
                   endif
 34             continue
 35             continue
c                print*,'inmax=',inmax
c                print*,'inzero=',inzero
                emader=energy(inmax)
c                print*,'emader=',emader
                dersh= energy(inzero) - derref
c                print*,'dersh=', dersh
                do 36 mm=1,nxmu
                   energy(mm)=energy(mm) - dersh
 36             continue

        endif        

        if (lfluo) then       
           if (langle) then
              angin = angin * (pi)/(180.0)
              angout = angout * (pi)/(180.0)
           else
              angin = 45.0 *(pi)/(180.0)
              angout = 45.0 *(pi)/(180.0)
           endif   
       
          if (angout.le.epsi)  angout = epsi      
            
          mutf = mutf* (sin(angin))/(sin(angout))
          alpha = ((mutf/muc) + (mub/muc))
          betha = ( (mutf/muc) + (mub/muc) + 1 )

         do 70 j=1,nxmu  
         xmuout(j) = xmuout(j)*(alpha)/(betha - xmuout(j))
 70      continue 
         
         do 71 j=1,nxmu
            xmuraw(j) = xmuout(j)
 71      continue 
         
         eefind = .true.
         call prenew(eefind, e0, predg1, predg2, enor1, enor2,
     $        nxmu, energy, xmuraw, step, stfind, slopre,
     $        bpre, xmuout)
         
         
         if (step.lt.epsi) step = 1.e0
         do 72 mm=1,nxmu
            xmuout(mm) = xmuout(mm) / step
 72      continue

        endif
 
        if (of.ne.' ') then
            outfil = of
            out2   = of
            if (frmout.eq.'ascii') then
                write(cnum,432)i
                jj=index(outfil,'.')
                if (jj .eq. 0) then
                   im=istrln(outfil)
                   outfil = outfil(1:im)//'.'//cnum
                   out2   = outfil(1:im)//'_der.'//cnum
                else
                   outfil = outfil(1:jj)//cnum
                   out2   = outfil(1:jj-1)//'_der.'//cnum
                endif
            endif    
        else
            if (frmout.eq.'ascii') then
                write(cnum,432)i
                ind = index(file(i), '.')
                ist = istrln(file(i))
                if (ind.ne.0) then
                    outfil = file(i)(:ind-1)//'.'//cnum
                    out2   = file(i)(:ind-1)//'_der.'//cnum
                else
                    outfil = file(i)(:ist)//'.'//cnum
                    out2   = file(i)(:ist)//'_der.'//cnum
                endif
            else
               ind = index(file(i), '.')
               ist = istrln(file(i))
               if (ind.ne.0) then
                  outfil = file(i)(1:ind-1)//'.nor'
                  out2   = file(i)(1:ind-1)//'_der'//'.nor'
               else
                  outfil = file(i)(1:ist)//'.nor'
                  out2   = file(i)(1:ist)//'_der'//'.nor'
               endif
            endif
        endif

        filtyp = 'xmu'
        nkeyo  = 0
        skeyo  = ' '
        iexist = 0
        call addocn(ndoc, doc, commnt(i), step, slopre, bpre)
        call outdat(filtyp, frmout, outfil, vaxflg, skeyo, nkeyo,
     $    ndoc, doc, nxmu, energy, xmuout, foo1, foo2, foo3, iexist)

        call addocd(ndoc, doc, commnt(i), maxder, emader)
        call outdat(filtyp, frmout, out2, vaxflg, skeyo, nkeyo,
     $    ndoc, doc, nxmu, energy, der, foo1, foo2, foo3, iexist)


        ii = istrln(outfil)
        if (frmout(1:2).eq.'uw') then
           call messag('  >>> Wrote uwexafs output to "'//
     $          outfil(:ii)//'"')
        else
           call messag('  >>> Wrote ascii output to "'//
     $          outfil(:ii)//'"')
        endif
        
 50   continue 
      call messag(' fluo is done. Have a nice day ')
      call messag(line)

      stop
c end main program fluo
      end


      subroutine noinit(nkey, predg1, predg2, enor1, enor2,
     $            foo1, foo2, foo3,
     $            of, file, skey, commnt, angin, angout, mutf,
     $            mub, muc)

      parameter(nptx=2**11, nfilx=50)

      character*10  skey(nfilx)
      character*78  file(nfilx), of, commnt(nfilx)
      integer       nkey(nfilx)
      real          mutf, mub, muc
      dimension energy(nptx), xmuout(nptx)
      dimension foo1(nptx), foo2(nptx), foo3(nptx)
      common /idata/ nxmu
      common /rdata/ e0, energy, xmuout
      parameter(zero=0.e0, nulval=-99 999)

c --------- integers ----------------------------------------------
      do 10 i=1,nfilx
        nkey(i) = 0
 10   continue 

c --------- reals -------------------------------------------------
      e0     = real(nulval-1)
      predg1 = -200.e0
      predg2 = -50.e0
      enor1  = 100.e0
      enor2  = 500.e0
      angin = 45.0
      angout = 45.0
      mutf = 10000.0
      mub = 10000.0
      muc = 10000.0

      do 20 i=1,nptx
        energy(i) = zero
        xmuout(i) = zero
        foo1(i)   = zero
        foo2(i)   = zero
        foo3(i)   = zero
 20   continue 

c --------- logicals ----------------------------------------------
c      eefind = .false.
c      lalign = .true.
c      lderiv = .false.

c --------- characters --------------------------------------------
      of = ' '
      do 510 i=1,nfilx
        file(i)   = ' '
        skey(i)   = ' '
        commnt(i) = ' '
 510  continue 

      return 
c  end subroutine noinit
      end 

      subroutine normin(nfile, nkey, file, e0, 
     $            enor1, enor2, predg1, predg2,
     $            frmin, frmout, of, skey, commnt,
     $            mutf, mub, muc, angin, angout,
     $            le0fin, lfluo, langle)

c--------------------------------------------------------------
c  copyright 1996 university of washington  Dani Haskel and Bruce Ravel
c--------------------------------------------------------------
      implicit real(a-h,o-z)
      implicit integer(i-n)
c      implicit double precision (a-h,o-z)
c---------------------------------------------------------------------   
c  this parses the lines of the command file looking for keywords then
c  reads in the value for that keyword from the next word.  
c---------------------------------------------------------------------

      parameter(nwdx=20, nfilx=50)
      parameter(zero=0.e0, nulval=-99999, imel=20)
      dimension     perc(imel)
      dimension     nkey(nfilx)
      character*2   test, cent, edge, elem(imel)
      character*10  frmin, frmout, skey(nfilx)
      character*20  fname, errmes*60
      character*78  file(nfilx), string, messg, of, words(nwdx),
     $              commnt(nfilx), temp
      logical       inpnul, there, le0fin, langle, lfluo, latoms, 
     $              latwei
      real          mutf, muc, mub

 1410 format(bn,f10.0)
 1440 format(a)

c---------------------------------------------------------------------
c  initialize some things used only in this routine
      iline  = 0
      mel    = 0
      inpnul = .true.
      test   = 'ab'
      le0fin = .false.
      langle = .false.
      lfluo = .false.
      latoms = .false.
      latwei = .false.
      frmin = ' '
      frmout = ' '

c  the value of the variable test must be in the same case as the 
c  keyword names in the long block of elseif's

c---------------------------------------------------------------------
c  open fluo.inp
      fname = 'fluo.inp'
      call upper(fname)
      inquire(file=fname,exist=there)
      if (.not.there) then
          call lower(fname)
          inquire(file=fname,exist=there)
      endif
      if (.not.there) then
          call messag(' ')
          call messag('* * * WARNING * * *')
          call messag('"fluo.inp" not found.  fluo stopping.')
          call messag(' ')
          stop
      endif
      open(unit=1,file=fname,status='old')

c---------------------------------------------------------------------
c  begin reading the input file, words cleared each time
101   continue
      nwds = nwdx
      do 110 iw=1,nwds
        words(iw)=' '
110   continue
      read (1,1440,end=191)string
      call untab(string)
      call triml(string)
c                                           - denotes end of job
      if  (string(1:1).eq.'-') goto 191
c                                           skip a comment line
      if  ((string(1:1).eq.'!').or.(string(1:1).eq.'*')
     $ .or.(string(1:1).eq.'%')) goto 101
c                                           skip a blank line
      if  (string.eq.' ') goto 101
c                                           begin reading line, i counts words
      i=1
      call bwords(string,nwds,words)
      inpnul = .false.
      iline  = iline+1
      do 120 ib = 1,nwds
         call triml(words(ib))
120   continue



c************************ input file parsing ***************************
c  read a word, identify it, assign the value from the following word(s)
c  increment i and come back.  i points to position in string, when i
c  exceeds nwds go read a new line.
 130  continue
      call case(test,words(i))
c                                           skip a blank line
      if     (words(i).eq.' ') then
         goto 101
c                                           ignore everything after !,*,%
      elseif ((words(i)(1:1).eq.'!').or.(words(i)(1:1).eq.'*').or.
     $          (words(i)(1:1).eq.'%')) then
          goto 101
c                                           finish early
      elseif (words(i)(1:3).eq.'end') then
          goto 191

c                                           central atom
      elseif (words(i)(1:4).eq.'cent') then
          cent=words(i+1)(1:2)
          i=i+2
c                                           edge of interest
      elseif (words(i).eq.'edge') then
          edge=words(i+1)(1:2)
          i=i+2
c                                           atomic weight flag
      elseif (words(i).eq.'atmwei') then
          if ((words(i+1)(1:1).eq.'t').or.(words(i+1)(1:1).eq.'T'))
     $        then
             latwei = .true.
             i=i+2
          endif   
c                                           e0
      elseif ((words(i).eq.'e0').or.(words(i).eq.'ee')) then
          call getrea(words(i), words(i+1), e0)
          le0fin = .true.
          i=i+2
c                                           pre-edge boundaries
      elseif (words(i).eq.'pre1') then
          call getrea(words(i), words(i+1), predg1)
          i=i+2
      elseif (words(i).eq.'pre2') then
          call getrea(words(i), words(i+1), predg2)
          i=i+2
c                                           normalization boundaries
      elseif (words(i).eq.'nor1') then
          call getrea(words(i), words(i+1), enor1)
          i=i+2
      elseif (words(i).eq.'nor2') then
          call getrea(words(i), words(i+1), enor2)
          i=i+2
      elseif (words(i).eq.'mutf') then
          call getrea(words(i), words(i+1), mutf)
          lfluo = .true.
          i=i+2
      elseif (words(i).eq.'mub') then
          call getrea(words(i), words(i+1), mub)
          i=i+2
      elseif (words(i).eq.'muc') then
          call getrea(words(i), words(i+1), muc)
          i=i+2
      elseif (words(i).eq.'angin') then
          call getrea(words(i), words(i+1), angin)
          langle = .true.
          i=i+2
      elseif (words(i).eq.'angout') then
          call getrea(words(i), words(i+1), angout)
          i=i+2

c                                           format
      elseif (words(i).eq.'format') then
          frmin  = words(i+1)
          frmout = words(i+1)
          call triml(frmin)
          call case(test,frmin)
          call triml(frmout)
          call case(test,frmout)
          i=i+2
c                                           formin
      elseif (words(i).eq.'formin') then
          frmin = words(i+1)
          call triml(frmin)
          call case(test,frmin)
          i=i+2
c                                           formout
      elseif (words(i).eq.'formout') then
          frmout = words(i+1)
          call triml(frmout)
          call case(test,frmout)
          i=i+2
c                                           outfile
      elseif (words(i)(1:3).eq.'out') then
          of = words(i+1)
          i=i+2

c                                              read in atoms list
      elseif (words(i)(1:4).eq.'atom') then
         mel    = 0
         do 528 ix=1,imel
               elem(ix) = ' '
               perc(ix) = zero
 528     continue
 529     read (1,1440,end=191)string
         call untab(string)
         call triml(string)
         if  (string(1:1).eq.'-') goto 191
c     skip a comment line
         if  ((string(1:1).eq.'!').or.(string(1:1).eq.'*')
     $        .or.(string(1:1).eq.'%')) goto 529
c     skip a blank line
         if  (string.eq.' ') goto 529
         nwds = nwdx
         call bwords(string,nwds,words)
         do 720 ib = 1,nwds
            call triml(words(ib))
 720     continue
         if (words(1)(1:3).eq.'fil') then 
            nfile=0
            goto 140
         else
            lfluo = .true.
            latoms = .true.
            mel=mel+1
            elem(mel)=words(1)(1:2)
            call getrea('atomic concentration',words(2),perc(mel))
c            print*, elem(mel), perc(mel)
            goto 529
         endif   

c                                               read in data filenames
      elseif ((words(i)(1:3).eq.'fil').or.(words(i)(1:3).eq.'dat'))
     $                then
          nfile = 0
 140      continue
          read(1,1440,end=191)string
          call untab(string)
          call triml(string)
          nnw = nwdx
          do 150 ii=1,nwdx
            words(ii)=' '
 150      continue
          call bwords(string,nnw,words)
          do 170 n=1,nnw
             call triml(words(n))
 170      continue
          if ((words(1)(:1).eq.'!').or.(words(1)(:1).eq.' ').or.
     $         (words(1)(:1).eq.'%').or.(words(1)(:1).eq.'*')) then
              goto 140
          elseif ((words(1)(1:3).eq.'---').or.(words(1).eq.'end')) then
              goto 191
          elseif ((words(1).eq.'e0').or.(words(1)(1:4).eq.'mess').or.
     $            (words(1).eq.'nofit').or.(words(1)(1:4).eq.'noal').or.
     $            (words(1).eq.'pre1').or.(words(1)(1:4).eq.'norm').or.
     $            (words(1).eq.'pre2').or.(words(1)(1:3).eq.'out').or.
     $            (words(1).eq.'nor1').or.(words(1)(1:3).eq.'fil').or.
     $            (words(1).eq.'nor2').or.(words(1).eq.'data').or.
     $            (words(1).eq.'emax').or.(words(1).eq.'format').or.
     $            (words(1).eq.'ee').or.
     $            (words(1).eq.'formin').or.(words(1).eq.'formout'))
     $                    then
              goto 130
          else
              nfile = nfile+1
              file(nfile) = words(1)
              call triml(file(nfile))

c             print*, ' frmin = ', frmin(1:5)
              nf = 2
              if (frmin(1:2).eq.'uw') then 
                  nf = 3
                  if (istrln(words(2)).le.3) then
                      call getint('nkey',words(2),nkey(nfile))
                  else
                      skey(nfile) = words(2)
                  endif
              endif

c will require user supplied file comment to be in same line after '*'

              if (frmin.eq.' ') then 
                 if ((words(2)(1:1).ne.'*').and.
     $               (words(2)(1:1).ne.' '))  then
                   nf = 3
                   if (istrln(words(2)).le.3) then
                      call getint('nkey',words(2),nkey(nfile))
                   else
                      skey(nfile) = words(2)
                   endif
                 endif
              endif


c              print*, ' nfile, nkey, skey = ', nkey(nfile), skey(nfile)
c             --- individual file comments
              if ((nnw.ge.nf).and.(words(nf)(1:1).ne.'!').and.
     $                    (words(nf)(1:1).ne.'%').and.
     $                    (words(nf)(1:1).ne.'#').and.
     $                    (words(nf)(1:1).eq.'*')) then
                  commnt(nfile) = string
                  do 160 k=1,nf-1
                    kk   = istrln(words(k))
                    temp = commnt(nfile)(kk+1:)
                    commnt(nfile) = temp
                    call triml(commnt(nfile))
 160              continue 
              endif
              goto 140
          endif
      else
          iunk  = istrln(words(i))
          messg = words(i)(:iunk)//' is an unknown keyword.'
          call messag(messg)
          i=i+2
      endif
c     if read entire line then read next line else read next word in line
      if (i.gt.nwds) goto 101
      goto 130

c     done reading lines
191   continue
      if (inpnul) then
         errmes = ' Your input file has no useful information '
          goto 666
      endif

c    calculate McMaster cross sections         
      if (latoms) then 
         print*, 'Atom list given ...'
         call mcfluo (cent, edge, mel, elem, perc,
     $        latwei, mutf, mub, muc)         
c         print*, mutf, mub, muc
      endif

      if (nfile.eq.0) then
          errmes = 'You specified no input data.'
          goto 666
      endif

300   continue
      return

 666  continue 
      ii = istrln(errmes)
      call messag('* * * ERROR!  '//errmes(:ii))
      call messag('      fluo is stopping.')
      stop

c end subroutine normin
      end

       subroutine prenew(eefind, ee, predg1, predg2, enor1, enor2,
     $     nxmu, energy, xmuraw, step, stfind, slopre, bpre, xmuout)
c
c     pre-edge subtraction and normalization of exafs data.
c
c inputs:
c   eefind   logical flag for picking e0 at the maximum derivative
c   ee       e0,  energy origin of data
c   predg1   region for picking pre-edge line
c   predg2   region for picking pre-edge line
c   enor1    region for picking normalization
c   enor2    region for picking normalization
c   nxmu     length of array energy, xmuraw, and xmuout
c   energy   array of energy points
c   xmuraw   array of raw absorption points
c   step     edge step for normalization (found if zero on input)
c   stfind   logical flag to tell whether to find edge step or not
c outputs:
c   ee       e0,  energy origin of data
c   predg1   region used for picking pre-edge line
c   predg2   region used for picking pre-edge line
c   enor1    region used for picking normalization
c   enor2    region used for picking normalization
c   step     edge step for normalization
c   slopre   slope of pre-edge line
c   bpre     intercept of pre-edge line (value at energy = 0)
c   xmuout   array of unnormalized absorption points
c            after pre-edge line subtracted
c
c   requires subroutine polyft (and so also function nofx)
c
c   copyright 1992  university of washington :          matt newville
c   Modified by Dani , march 96. Changes include checking pre1 and pre2 
c   and nor1, nor2 to be in data range. Also e0 has to be in data range.
c   Criteria changed to at least 6 consecutive increasing points.
c-----------------------------------------------------------------
      integer     ninc
      parameter   (ninc = 6)
      logical     eefind, stfind, inc(ninc)
      dimension   energy(nxmu), xmuraw(nxmu), xmuout(nxmu), anorm(2)
      real        zero, tiny, one, onepls
      parameter (zero = 0.0, tiny = 1.e-10, one = 1.0, onepls = 1.01)
c
c   if ee was not specified, it is found in the data from
c   these conditions:
c       1) ee is the third increasing point in a row.
c       2) ee has the largest derivative up to that point.
c   these conditions will almost always pick a point on the
c   absorption step, and will fail less often than simply
c   taking the point with maximum derivative.
c
      if (nxmu.le.5)  return

      if ((.not. eefind).and. ((ee.lt.energy(1)).or.
     $     (ee.gt.energy(nxmu)))) then
ccccc      if (ee.lt.energy(1)) then
         eefind =.true.
         print*, ' The value of e0 given is outside data range'
         print*, ' fluo is finding its own value'
      endif

      if (eefind) then
         do 100 j = 1, ninc
            inc(j) = .false.
 100     continue
         dxde  = zero
         demx  = zero
         ntry  = max(2, int(nxmu/2)) + 3
         do 150 i = 2, ntry
            deltax  = xmuraw(i) - xmuraw(i-1)
            deltae  = max(tiny, (energy(i) - energy(i-1)))
            dxde    = deltax/deltae
            inc(1)  = (dxde.gt.zero).and.(deltae.gt.tiny)
            if ((dxde.gt.demx).and.inc(6).and.inc(5).and.inc(4)
     $           .and.inc(3).and.inc(2).and.inc(1)) then
               ee   = energy(i)
               demx = dxde * onepls
            end if
            if (inc(6).and.(.not.inc(5)).and.(.not.inc(4))
     $           .and.(.not.inc(3)).and.(.not.inc(2))
     $           .and.(.not.inc(1))
     $                .and.(demx.gt.zero)) go to 200
            do 130 j  = ninc, 2, -1
               inc(j) = inc(j - 1)
 130        continue
 150     continue
      end if
 200  continue
cc      print*, ' pre-edge: loop stopped at = ', i 
c  at this point, ee is known. so we go on
c  to do the pre-edge  and normalization
c
c  get pre-edge line (slope and intercept) by linear fit
      if ( (predg1.eq.zero).and.(predg2.eq.zero) )  then
         predg2 = -50
         predg1 = -200
      end if
      pre1  = ee + predg1
      pre2  = ee + predg2
      if (pre1.lt.energy(1))   pre1 = energy(2)
      if (pre2.lt.energy(1))   pre2 = ( ee + pre1 ) /2
      call polyft(pre1, pre2, energy, xmuraw, nxmu, 2, anorm)
      bpre   = anorm(1)
      slopre = anorm(2)
c
c  apply e0 shift, do pre-edge subtraction, and calculate
c  the energy relative to edge for easier normalization
      do 300 i = 1, nxmu
         xmuout(i)  = xmuraw(i) - bpre - slopre*energy(i)
         energy(i)  = energy(i) - ee
 300  continue
      
c  normalization : make pre-edge 0.0 and post-edge 1.0
c    if step size wasn't given, get it by extracting to ee a
c    line that best fits the data on the range (ee+enor1,ee+enor2)
      if (stfind) then
         if ( (enor1.eq.zero).and.(enor2.eq.zero) )   then
            enor1 = 100
            enor2 = 300
         end if
         if (enor2.gt.energy(nxmu))   enor2 = energy(nxmu)
         if (enor1.gt.energy(nxmu))   enor1 = ( ee + enor2 ) /2
         call polyft(enor1, enor2, energy, xmuout, nxmu, 2, anorm)
         step   = anorm(1)
      end if
      if (abs(step).lt.tiny) step = one
c  put energy array back to original values, 
c  recover from bad pre-edge fits
      do 500 i = 1, nxmu
         energy(i)  = energy(i) + ee
c     for "post-edge" pre-edge fit:
         if (predg1.gt.zero) then
            xmuout(i) = one + xmuout(i)
         end if
 500  continue
      return
c end subroutine prenew
      end

      subroutine addocn(ndoc, doc, commnt, step, slopre, bpre)
      
      parameter(ndocx=25)
      character*78  commnt
      character*100 doc(ndocx)

c     If file comment is given, need to make room for 2 new lines; if not 
c     just need one document line. Copy line i to line i+2 (or i+1). The 
c     first line of doc in not altered.

      nnd = 2
      if (commnt.ne.' ') then
        do 10 i=min(ndoc, ndocx-nnd), nnd, -1
         doc(i+2) = doc(i)
 10     continue 

       write(doc(nnd),400)commnt(2:)
 400   format('fluo: ',a)

       write(doc(1+nnd),410)step, slopre, bpre
 410   format('fluo: STEP = ',f7.4', pre-edge line: mux = ',
     $            f9.7' *E +',f9.5)

      else

      do 20 i=min(ndoc, ndocx-nnd), nnd, -1
        doc(i+1) = doc(i)
 20     continue
        write(doc(nnd),410)step, slopre, bpre
      endif

      return 
c  end subroutine addoc      
      end




c----------------------------------------------------------------------
c          input/output routines for data files
c               for the uwxafs programs
c
c   input/output routines for data files for the uwxafs programs
c
c   copyright 1992  university of washington, seattle, washington
c   written by      matthew newville
c                   department of physics, fm-15
c                   university of washington
c                   seattle, wa   usa 98195
c   phone           (206) 543-0435
c   e-mail          newville@u.washington.edu
c
c  these routines are the basic input/output routines for getting
c  numerical and document data from files into the uwxafs programs.
c  there are currently two data formats supported:
c
c 1. 'uw' :  a binary file format known as the uwxafs file handling
c            routines. this is very efficient way to store data, and
c            can store several (191) data sets in a single file. the
c            drawback is that the files are not extremely portable.
c
c 2. 'asc':  these are column files in a format that is fairly easy
c            for anything to deal with. the files have several lines
c            of documents. if the first character of the document is
c            '#' this character will be removed. after the documents
c            is a line with minus signs for characters(3:6), then an
c            ignored line (for column labels), and then the data. up
c            to five columns are used. the expected order is:
c                  x, real(y), imag(y), ampl(y), phase(y).
c            if any column representing y is zero, the appropriate
c            value will be calculated and returned. the files in this
c            format hold only one data set, and use more memory than
c            the uwxafs files, but are portable and convenient.
c
c  other file types can be added without too much difficulty.
c  the routines listed here are:
c      inpdat : retrieve data and documents from a file
c      inpcol : retrieve data and documents from an ascii file
c      inpuwx : retrieve data and documents from a uwxafs file
c      outdat : write data and documents to a file
c      outcol : write data and documents to an ascii file
c      outuwx : write data and documents to a uwxafs file
c
c  note: the fortran input/output unit number 11 is used for all
c        unit numbers in these routines. conflicts between these
c        routines will not happen, but conflicts may arise if
c        unit = 11 indicates an open file in a calling subprogram.
c----------------------------------------------------------------------
       subroutine inpdat(filtyp, format, filnam, vax, skey, nkey, 
     $       ndoc, doc, ndata, xdata, yreal, yimag, yampl, yphas )
c
c   copyright 1992  university of washington :          matt newville
c
c    retrieve data and documents from a file acording
c    to the format specified by 'format'.
c inputs:
c   filtyp    file type to open. if may be ' '
c   format    file format (uwxafs, ascii, column)
c   filnam    file name
c   vax       logical flag for being on a vax machine (binary file)
c   skey      symbolic key for record in uwxafs file
c   ndata     maximum number of elements in data arrays
c   nkey      numeric key for record in uwxafs file
c   ndoc      maximum number of document lines to get
c               note:   ndoc cannot be less than or equal to zero!
c outputs:
c   skey      symbolic key of record in uwxafs file
c   ndoc      number of document lines returned
c   doc       array of document lines
c   ndata     number of elements in data arrays
c   xdata     array of x values of data
c   yreal     array of real part of y data values
c   yimag     array of imaginary part of y data values
c   yampl     array of amplitude part of y data values
c   yphas     array of phase part of y data values
c---------------------------------------------------------------------
       character*(*)  filtyp, format, skey, filnam, doc(*)
       character*10   type, symkey, form, formin, errmsg*80
       dimension      xdata(*), yreal(*), yimag(*)
       dimension      yampl(*), yphas(*)
       logical        vax 
       integer        irecl, ndatmx, ndocmx
       data  irecl, ndocmx, ndatmx  / 512 , 19, 4096/
c---------------------------------------------------------------------
c some initializations
       if (vax) irecl = 128   
       type = filtyp
       call triml(type)
       call upper(type)
       symkey = skey
       call triml(symkey)
       call upper(symkey)
c       
c determine format of the input file
       formin = format
       call triml(formin)
       call smcase(formin, 'a')
       call testrf(filnam, irecl, form, ier)
       call smcase(form, 'a')
       if (ier.eq.-1) then 
          call messag('  inpdat error: file not found  ')
       elseif (ier.eq.-2) then 
          call messag('  inpdat error: unknown file format = '//formin)
       elseif (ier.eq.-3) then 
          call messag('  inpdat error: poorly formatted ascii data?  ')
       elseif (ier.eq.-4) then 
          call messag('  inpdat error: no data in ascii format file? ')
       end if
       if (ier.ne.0) then 
          errmsg =    '    for file ' // filnam
          ilen   = istrln(errmsg)
          call messag( errmsg(1:ilen) )
          stop
       endif
       if ((formin.ne.' ').and.(formin(1:2).ne.form(1:2))) then
          call messag('  inpdat warning: the requested format was'//
     $         ' incorrect!')
          call messag('  form    = '//form(1:5)  )
          call messag('  formin  = '//formin(1:5)  )
       end if
c  now call the appropriate routine to get the data,
c  according to the format.
       ndata = max(1, min(ndata, ndatmx) )
       ndoc  = max(1, min(ndoc , ndocmx) )
       if (form(1:2).eq.'uw') then
          ndoc = ndocmx
          call inpuwx(type, filnam, skey, nkey, irecl, ndoc, doc,
     $         ndata, xdata, yreal, yimag, yampl, yphas )
       elseif ((form(1:2).eq.'co').or.(form(1:2).eq.'as')) then
          call inpcol(filnam, ndoc, doc,
     $         ndata, xdata, yreal, yimag, yampl, yphas )
          skey   = 'ascii'
          call upper(skey)
       else
          call messag('  inpdat error: unknown file format = '// form)
          ilen   = min(54, max(1, istrln(filnam)))
          errmsg = '                for file ' // filnam(1:ilen)
          call messag( errmsg(1:ilen+26) )
          stop
       end if
       filtyp = type
       format = form
c     
       return
c end subroutine inpdat
       end
       subroutine inpcol(filnam, ndoc, doc, ndata,
     $                   xdata, yreal, yimag, yampl, yphas)
c
c   copyright 1992  university of washington :       matt newville
c
c   open and get all information from a column file. document
c   lines are read until a line of '----', then a label line is
c   skipped and the column data are read in.  the data is read 
c   and stored in the following order:
c                xdata  yreal  yimag  yampl  yphas
c inputs:
c   filnam    file name containing data
c   ndoc      maximum number of document lines to get
c   ndata     maximum number of elements in data arrays
c outputs:
c   ndoc      number of document lines returned
c   doc       array of document lines
c   ndata     number of elements in data arrays
c   xdata     array of x values of data
c   yreal     array of real part of y data values
c   yimag     array of imaginary part of y data values
c   yampl     array of amplitude part of y data values
c   yphas     array of phase part of y data values
c---------------------------------------------------------------------
       integer   ilen , istrln, j, i, mxword, ndoc, ndata, iounit
       integer   iexist, ierr, nwords, idoc, id
       real      zero , small
       parameter( zero = 0.e0 , small = 1.e-6, mxword = 5)
       real      xdata(*), yreal(*), yimag(*), yampl(*), yphas(*)
       real      xinp(mxword)
       logical   isdat
       character*(*) filnam, doc(*)
       character*30  words(mxword), line*100, status*10, file*60
       external      istrln, isdat
c--------------------------------------------------------------------- 
 10    format(a)
       file = filnam
       ilen = istrln(file)
       if (ilen.le.0)  then
           call messag( ' inpcol:  no file name given')
           stop
       end if
c  initialize buffers
       do 80 j = 1, ndoc
          doc(j) = ' '
  80   continue
       do 100 i = 1, mxword
          words(i) = '0. '
          xinp(i)  = zero
 100   continue
       do 120 j = 1, ndata
          xdata(j) = zero
          yreal(j) = zero
          yimag(j) = zero
          yampl(j) = zero
          yphas(j) = zero
 120   continue
c  open data file
      iounit = 7
      status ='old'
      call openfl(iounit, filnam, status, iexist, ierr)
      if ((iexist.lt.0).or.(ierr.lt.0)) go to 900
c
c  get documents from header: up to ndoc
c       read file header, save as document lines,
c       remove leading '#' and '%' both of which are
c       known to be extraneous comment characters.
       nwords = 5
       idoc = 0
       id   = 1
 200   continue
          read(iounit, 10, end = 950, err = 960) line
          call triml (line)
c  if line is '----', read one more line, go read numerical data
          if (line(3:6) .eq. '----')  then
             read(iounit, 10, end = 950, err = 960) line
             goto 400
          end if
c  remove leading '#' or '%' from line 
          if ( (line(1:1).eq.'#').or.(line(1:1).eq.'%') ) then
             line(1:1) = ' '
             call triml(line)
c  if the line is all numbers, then this is data!
          elseif (isdat(line)) then
             goto 410
          end if
c  save line in doc if there's room
          if ((idoc .lt. ndoc) .and. (istrln(line).gt.0) ) then
             idoc = idoc + 1
             doc(idoc) = line
          endif
          goto 200
c     
c  read numerical data
 400   continue
          nwords = 5
          read(iounit, 10, end = 600, err = 980) line
 410      continue 
          call untab(line)
          call bwords(line,nwords,words)
          if (nwords.le.1) goto 600 
          do 450 i = 1, nwords
c >>> changed for phit
              call sgetvl(words(i), xinp(i), ierr) 
c <<<
              if (ierr.ne.0) goto 600
 450      continue
          xdata(id) = xinp(1)
          yreal(id) = xinp(2)
          yimag(id) = xinp(3)
          yampl(id) = xinp(4)
          yphas(id) = xinp(5)
          if (id.ge.ndata) go to 610
          if ( id.ge.2 ) then
             if (abs(xdata(id) - xdata(id-1)).ge.small) id = id + 1
          else
             id = id + 1
          end if
          goto 400
 600   continue
       id    = id - 1
       if (id.lt.1) go to 950
 610   continue
       ndata = id 
       if (idoc.le.0) then
          ndoc =  1
          doc(1) = 'inpdat: no document line found'
       end if
c  make sure that all columns are filled:
c   if yampl and yphas are both zero, compute them from yreal, yimag
c   if yreal and yimag are both zero, compute them from yampl, yphas
       do 800 i = 1, ndata
          if ( ( (yampl(i).eq.zero).and.(yphas(i).eq.zero) ) .and.
     $         ( (yreal(i).ne.zero).or. (yimag(i).ne.zero) ) ) then
            yampl(i) = sqrt( yreal(i)**2 + yimag(i)**2 )
            yphas(i) = atan2( yimag(i), yreal(i) )
             if (i.gt.1) call pijump( yphas(i), yphas(i-1) )
 
          elseif ( (yreal(i).eq.zero).and.(yimag(i).eq.zero)
     $        .and.(yampl(i).ne.zero)   ) then
            yreal(i) = yampl(i) * cos ( yphas(i) )
            yimag(i) = yampl(i) * sin ( yphas(i) )
 
          end if
 800   continue
c  close data file and return
       close(iounit)
       return
c error handling
c  open file - error
 900   continue
         call messag(' inpcol: error opening file '//file(1:ilen) )
         go to 990
c  end or error at reading documents 
 950   continue
 960   continue
         call messag( ' inpcol: error reading file '//file(1:ilen) )
         call messag('         during reading of documents.') 
         go to 990
c  error at reading numerical data
 980   continue
         call messag( ' inpcol: error reading file '//file(1:ilen) )
         call messag('         during reading of numerical data.')

 990     continue 
         close(iounit)
         stop
c end error handling
c end subroutine inpcol
       end
       subroutine inpuwx(ftypin, filein, skey, nkey, irecl, ndoc, 
     $           documt, ndata, xdata, yreal, yimag, yampl, yphas )
c
c   copyright 1992  university of washington :          matt newville
c
c     open and get all information from a uwxafs file
c
c inputs:
c   ftypin   file type to open, checked for compatibility, may be ' '
c   filein   file name containing data
c   skey     symbolic key for record in data file (only one of these)
c   nkey     numeric key for record in data file  (two is needed    )
c   ndoc     maximum number of document lines to get
c outputs:
c   skey      symbolic key of record in uwxafs file
c   ndoc      number of document lines returned
c   docu      array of document lines
c   ndata     number of elements in data arrays
c   xdata     array of x values of data
c   yreal     array of real part of y data values
c   yimag     array of imaginary part of y data values
c   yampl     array of amplitude part of y data values
c   yphas     array of phase part of y data values
c
c notes:
c  1   the full 'noabort' error checking is done for the calls to
c      the uwxafs routines, which means that marginally useful
c      error messages will be given when one of the uwxafs
c      filehandling routines dies.
c
c  2    currently, the following file types are supported:
c           xmu,  chi,  rsp,  env,  rspep, rip
c
c  3    uwxafs file handling routines only do single precision.
c       this routine can be made implicit double precision if the
c       array buffer is maintained as single precision:
c           implicit double precision(a-h,o-z)
c           real           buffer(maxpts)
c---------------------------------------------------------------------
       parameter( maxpts = 2048, zero = 0. )
       character*(*)  ftypin, skey, filein, documt(*)
       character*10   type, ftype, safefl*8, abrtfl*8
       character*70   filnam, messg
       dimension      xdata(*), yreal(*), yimag(*)
       dimension      yampl(*), yphas(*)
       real           buffer(maxpts)
c---------------------------------------------------------------------
c initialize
 10    format(a)
 20    format(2x,2a)
 30    format(2x,a,i3)
       safefl = ' '
       abrtfl = 'noabort'
       call upper(abrtfl)
 
       ftype = ftypin
       filnam= filein
 
       call upper(skey)
       call triml(skey)
       call triml(ftype)
       call triml(filnam)
       ilen = max(1, istrln(filnam))
c note: uwxafs requires ftype to be upper case.
       call upper (ftype)
        do 100 i = 1,ndata
            xdata(i)  = zero
            yreal(i)  = zero
            yimag(i)  = zero
            yampl(i)  = zero
            yphas(i)  = zero
100    continue
       do 110 i = 1, maxpts
            buffer(i) = zero
110    continue
c  call uwxafs file handling routines:
c : open data file
       iounit = 11
       call openrf(iounit, filnam, abrtfl, safefl, ftype, irecl, ier)
       if (ier.ne.0) then
                messg = 'inpuwx: error opening file '
           call messag(messg//filnam(:ilen))
                write (messg, '(9x,a,i4)') 'openrf error code ',ier
           call messag(messg)
           stop
       end if
c : check file type
       call gftype(iounit, type, ier)
       if (ier.ne.0) then
                messg = 'inpuwx: error getting file type for '
           call messag(messg//filnam(:ilen))
                write (messg, '(9x,a,i4)') 'gftype error code ',ier
           call messag(messg)
           stop
       end if
       call upper(type)
 
       if (ftype.eq.' ') then
           ftype = type
       elseif (ftype.ne.type) then
                messg = 'inpuwx: incorrect file type for '
           call messag(messg//filnam(:ilen))
                messg = '     file type for this file is '
           call messag(messg//type)
                messg = '     file type requested was '
           call messag(messg//ftype)
           stop
       endif
       ftypin = ftype
 
c : find out how many records there are in the file
       call gnie (iounit, nie, ier)
       if (nie.le.0) then
               messg = 'inpuwx:  no data records in '
          call messag(messg//filnam(:ilen) )
          stop
       end if
c : get skey if it wasn't given as input
       if (skey.eq.' ') then
           call gskey(iounit, nkey, skey, ier)
           if (ier.ne.0) then
                  messg = 'inpuwx: error getting skey for '
             call messag(messg//filnam(:ilen))
                  write (messg, '(9x,a,i4)') 'gskey error code ',ier
             call messag(messg)
             stop
           end if
           if (skey.eq.' ') then
             write (messg, '(1x,2a,i4)') 'inpuwx: found no skey ',
     $                                  'for nkey =',nkey
             call messag(messg)
             call messag('        in file = '//filnam(:ilen))
             stop
           end if
       end if
 
c : get nkey if it wasn't given as input
       if (nkey.eq.0) then
           call gnkey(iounit, skey, nkey, ier)
           if (ier.ne.0) then
                  messg = 'inpuwx: error getting nkey for '
             call messag(messg//filnam(:ilen))
                  write (messg, '(9x,a,i4)') 'gnkey error code ',ier
             call messag(messg)
             stop
           end if
       end if
c 
c : get documents : up to ndoc
c   first check how many document lines there are
       call gdlen(iounit, nkey, ndocln, ier)
       if (ier.ne.0) then
               messg = 'inpuwx: error getting document length for '
          call messag(messg//filnam(:ilen))
               write (messg, '(9x,a,i4)') 'gdlen error code ',ier
          call messag(messg)
          stop
       end if
       if (ndoc.gt.ndocln) ndoc = ndocln
c   then get the documents
       call getdoc(iounit, documt, ndoc, skey, nkey, ndsent, ier)
       if (ier.eq.6) then
               messg = 'inpuwx error: reading file '
          call messag(messg//filnam(:ilen) )
               messg = '  no skey or nkey given to specify record, '
          call messag(messg)
               messg = '  or an incorrect skey or nkey given '
          call messag(messg)
          stop
       elseif (ier.ne.0) then
               messg = 'inpuwx: error getting documents for '
          call messag(messg//filnam(:ilen))
               write (messg, '(9x,a,i4)') 'getdoc error code ',ier
          call messag(messg)
          stop
       end if
       ndoc = ndsent
 
c : get data
       call getrec(iounit, buffer, maxpts, skey, nkey, nbuff, ier)
       if (ier.ne.0) then
               messg = 'inpuwx: error getting data for '
          call messag(messg//filnam(:ilen))
               write (messg, '(9x,a,i4)') 'getrec error code ',ier
          call messag(messg)
          stop
       end if
 
c : close file
       call closrf(iounit,ier)
       if (ier.ne.0) then
               messg = 'inpuwx: error closing data file '
          call messag(messg//filnam(:ilen))
               write (messg, '(9x,a,i4)') 'closrf error code ',ier
          call messag(messg)
          stop
       end if
c-----------------------------------------------------------------
c finished with uwxafs routines, so now sort the data into
c xdata, re(y), imag(y), ampl(y), phase(y) according to file type
c
c convert ftype to the case of this routine.
c   'case' controls the the case of this routine
       call smcase (ftype, 'case')
c- xmu: nbuff energy, then nbuff y-values
       if (ftype.eq.'xmu') then
            ndata   = nbuff/2
            do 400 i = 1, ndata
               xdata(i) = buffer(i)
               yreal(i) = buffer(ndata + i)
               yampl(i) = yreal(i)
400         continue
c-  chi: xmin, deltax, chi(kmin + i*deltak)
       elseif (ftype.eq.'chi') then
            ndata   = nbuff - 2
            do 500 i = 1, ndata
               xdata(i) = buffer(1) + (i-1)*buffer(2)
               yreal(i) = buffer(2 + i)
               yampl(i) = yreal(i)
500         continue
c-  env,rspep: kmin, deltak, phase, amplitude pairs (kmin + i*deltak)
       elseif ( (ftype.eq.'env').or.(ftype.eq.'rspep')  ) then
            ndata   = (nbuff - 1) / 2
            do 600 i = 1, ndata
               xdata(i) = buffer(1) +(i-1)*buffer(2)
               yphas(i) = buffer(2*i+1)
               yampl(i) = buffer(2*i+2)
               yreal(i) = yampl(i) * cos ( yphas(i) )
               yimag(i) = yampl(i) * sin ( yphas(i) )
600         continue
c  rsp, rip: kmin, deltak, real, imaginary pairs (kmin + i*deltak)
       elseif ( (ftype.eq.'rsp').or.(ftype.eq.'rip')  ) then
            ndata   = (nbuff - 1) / 2
            do 700 i = 1, ndata
               xdata(i) = buffer(1) +(i-1)*buffer(2)
               yreal(i) = buffer(2*i+1)
               yimag(i) = buffer(2*i+2)
               yampl(i) = sqrt( yreal(i)**2 + yimag(i)**2 )
               yphas(i) = atan2( yimag(i), yreal(i) )
                 if (i.gt.1) call pijump( yphas(i), yphas(i-1) )
700         continue
       else
                  messg = 'inpuwx: unrecognized file type for '
             call messag(messg//filnam(:ilen))
                  messg = '        file type for this file is '
             call messag(messg//ftype)
             stop
       end if
       return
c end subroutine inpuwx
       end
       subroutine outdat(filtyp, format, filnam, vax, skey, nkey,
     $    ndoc, doc, ndata, xdata, yreal, yimag, yampl, yphas, iexist)
c
c   copyright 1992  university of washington :          matt newville
c
c    write data and documents to a file acording to the
c    format specified
c inputs:
c   filtyp    file type to open, may be ' '.
c   format    file format (uwxafs, ascii, column)
c   filnam    file name
c   vax       logical flag for being on a vax machine (binary file)
c   ndoc      number of document lines returned
c   doc       array of document lines
c   ndata     number of elements in data arrays
c   xdata     array of x values of data
c   yreal     array of real part of y data values
c   yimag     array of imaginary part of y data values
c   yampl     array of amplitude part of y data values
c   yphas     array of phase part of y data values
c   iexist    flag for whether to write redundant data to uwxafs file
c             iexist = 1 : do not write redundant data
c             iexist = 0 : do write redundant data
c outputs:
c   skey      symbolic key for record in uwxafs file
c   nkey      numeric key for record in uwxafs file
c   ndoc      number of document lines written
c
c---------------------------------------------------------------------
       character*(*)  filtyp, format, filnam, skey, doc(*)
       character*30   type, form
       dimension      xdata(*), yreal(*), yimag(*)
       dimension      yampl(*), yphas(*)
       logical        vax 
       integer        irecl
c---------------------------------------------------------------------
       irecl  = 512
       if (vax) irecl = 128   
c
       idoc = ndoc
       skey = ' '
       form = format
       type = filtyp
       call upper(type)
       call triml(type)
       call triml(form)
c convert form to the case of this routine.
c   'case' controls the the case of this routine
       call smcase (form, 'case')
c
       if (form(1:2).eq.'uw') then
             if ( (idoc.le.0).or.(idoc.gt.19) )  idoc = 19
             call outuwx(type, filnam, skey, nkey, irecl, idoc, doc,
     $           ndata, xdata, yreal, yimag, yampl, yphas, iexist)
       elseif ( (form(1:3).eq.'col').or.(form(1:3).eq.'asc') ) then
              call outcol(type, filnam, idoc, doc,
     $           ndata, xdata, yreal, yimag, yampl, yphas)
              skey = 'ascii'
             call upper(skey)
       else
         call messag('outdat: unknown file format = '//form)
         stop
       end if
c
       return
c end subroutine outdat
       end
       subroutine outcol(filtyp, filnam, ndoc, doc, ndata,
     $                  xdata, yreal, yimag, yampl, yphas)
c
c   copyright 1992  university of washington :          matt newville
c
c  open and write all information to a column file. document lines are
c  written, followed by a line of '----', then a label line, and then
c  the data are written.  the file type tells what to use for the label
c  and how many columns to write. it may be left blank.
c
c inputs:
c   filtyp    file type to write (may be ' ' : used for label only)
c   filnam    file name to write (' ' and '*' mean write to unit 6)
c   ndoc      maximum number of document lines to write
c   documt    array of document lines
c   ndata     number of elements in data arrays
c   xdata     array of x values of data
c   yreal     array of real part of y data values
c   yimag     array of imaginary part of y data values
c   yampl     array of amplitude part of y data values
c   yphas     array of phase part of y data values
c outputs:
c   ndoc      number of document lines written
c
c---------------------------------------------------------------------
       parameter(zero = 0.000000)
       character*(*)  filtyp, filnam, doc(*)
       character*80   filout, errmsg, type*10, status*10
       dimension      xdata(*), yreal(*), yimag(*)
       dimension      yampl(*), yphas(*)
c---------------------------------------------------------------------
       type   = filtyp
       call triml(type)
c convert type to the case of this routine.
       call smcase(type, 'a')
       filout = filnam
       call triml(filout)
       if (ndata.le.0) ndata = 2
 25    format(2a)
 26    format(a,5x,a)
c open data file
c     if file name is ' ' or '*', write to standard output (unit 6)
       iounit = 6
       if ((filout.ne.' ').and.(filout.ne.'*')) then
          iounit = 0
          status ='unknown'
          call openfl(iounit, filout, status, iexist, ierr)
          if ((ierr.lt.0).or.(iexist.lt.0)) go to 990
       endif
c     
c     write documents to header
       jdoc  = 0
       mxl   = 76
       mxlp1 = mxl + 1
       do 200, idoc = 1, ndoc
          call triml(doc(idoc))
          ilen = istrln(doc(idoc))
          if (ilen.ge.1) then
             jdoc = jdoc + 1
             if (ilen.gt.mxl) then
                write(iounit, 25) '# ', doc(idoc)(1:mxl)
                write(iounit, 26) '# ', doc(idoc)(mxlp1:ilen)
             else
                write(iounit, 25) '# ', doc(idoc)(1:ilen)
             end if
          end if
 200   continue
       ndoc = jdoc
       write(iounit, 25) '#----------------------------------',
     $                   '-----------------------------------'
c  write line of minus signs and the label for columns
       if (type.eq.'xmu') then
          write(iounit, 25) '#     energy          xmu       ', ' '
c >>> changed for phit
       elseif (type.eq.'2col') then
          write(iounit, 25) '#     x               y(x)       ', ' '
c <<<
       elseif (type.eq.'chi') then
          write(iounit, 25) '#     k              chi(k)     ', ' '
       elseif (type.eq.'env') then
          write(iounit, 25) '#     k          real(chi(k))   ',
     $         'imag(chi(k))   ampl(chi(k))   phase(chi(k))'
       elseif ( (type.eq.'rsp').or.(type.eq.'rspep') ) then
          write(iounit, 25) '#     r          real(chi(r))   ',
     $         'imag(chi(r))   ampl(chi(r))   phase(chi(r))'
       else
          write(iounit, 25) '#     x          real(y(x))     ',
     $         'imag(y(x))     ampl(y(x))     phase(y(x))  '
       end if
c
c  make sure that all of re(y), im(y), amp(y), and phase(y) are known
       do 300 i = 1, ndata
          if ( ( (yampl(i).eq.zero).and.(yphas(i).eq.zero) ) .and.
     $         ( (yreal(i).ne.zero).or. (yimag(i).ne.zero) ) ) then
             yampl(i) = sqrt( yreal(i)**2 + yimag(i)**2 )
             yphas(i) = atan2( yimag(i), yreal(i) )
             if (i.gt.1) call pijump( yphas(i), yphas(i-1) )
          elseif ( (yreal(i).eq.zero).and.(yimag(i).eq.zero)
     $                               .and.(yampl(i).ne.zero) ) then
             yreal(i) = yampl(i) * cos ( yphas(i) )
             yimag(i) = yampl(i) * sin ( yphas(i) )
          end if
 300   continue
c
c  write out data : some file types only write out a few columns
 380   format(2x,e13.7,3x,e13.7)
 390   format(2x,e13.7,2x,e13.7,2x,e13.7,2x,e13.7,2x,e13.7)
       if ( (type.eq.'xmu').or.(type.eq.'chi').or.
     $             (type.eq.'2col') ) then
          do 400 i = 1, ndata
             write(iounit, 380) xdata(i), yreal(i)
 400      continue
       else
          do 450 i = 1, ndata
             write(iounit, 390) xdata(i), yreal(i), yimag(i),
     $                          yampl(i), yphas(i)
 450      continue
       end if
c
c  close data file and return
       close(iounit)
       return
 990   continue
       ilen   = max(1, istrln(filnam))
       errmsg = 'outcol: error opening file '//filnam(:ilen)
       imsg   = istrln(errmsg)
       call messag(errmsg(:imsg))
       stop
c end subroutine outcol
       end

       subroutine outuwx(ftypin, filein, skey, nkey, irecl, ndoc, doc,
     $           ndata, xdata, yreal, yimag, yampl, yphas, iexist)
c
c     write out data and documents to a uwxafs file
c
c inputs:
c   ftypin    file type to write to, may be ' ' if filnam exists.
c   filein    file name to write to
c   skey      symbolic key of record in uwxafs file
c   ndoc      number of document lines returned
c   doc       array of document lines
c   ndata     number of elements in data arrays
c   xdata     array of x values of data
c   yreal     array of real part of y data values
c   yimag     array of imaginary part of y data values
c   yampl     array of amplitude part of y data values
c   yphas     array of phase part of y data values
c   iexist    flag for whether to write redundant data to file
c               iexist = 1 : do not write redundant data
c               iexist = 0 : do write redundant data
c
c   copyright 1992  university of washington :          matt newville
c-----------------------------------------------------------------------
       parameter(maxpts = 2048, maxdoc = 19, zero = 0.)
       character*(*)  filein, ftypin, doc(*), skey
       character*10   skyout, ftype, type, filnam*70, messg*70
       character*100  docout(maxdoc), abrtfl*8, safefl*8
       dimension      xdata(*), yreal(*), yimag(*)
       dimension      yampl(*), yphas(*)
       real           buffer(maxpts)
c-----------------------------------------------------------------------
c initialize
 10    format(a)
       safefl = ' '
       abrtfl = 'noabort'
       call upper(abrtfl)
       skyout = ' '
       type   = ' '
       filnam = filein
       call triml(filnam)
       ilen   = max(1, istrln(filnam))
       ftype  = ftypin
       call upper(ftype)
c >>> new for phit
       if (ftype.eq.'2COL') ftype = 'XMU'
c <<<
       do 60 i = 1, maxdoc
          docout(i) = ' '
 60    continue
c output documents
       idoc = 0
 80    continue
          idoc  = idoc + 1
          if ((idoc.ge.maxdoc).or.(idoc.gt.ndoc)) then
             idoc  = idoc - 1
             go to 100
          end if
          docout(idoc) = doc(idoc)
          call triml(docout(idoc))
          go to 80
100    continue
ccccc       ndoc = idoc
c  open data file to check file type
       iounit = 11
       call openrf(iounit, filnam, abrtfl, safefl, ftype, irecl, ier)
       if (ier.ne.0) then
               messg = 'outuwx: error opening file '//filnam(:ilen)
               imsg  = max(1, istrln(messg))
          call messag(messg(:imsg))
               write(messg, '(9x,a,i3)' ) 'openrf error code ',ier
          call messag(messg)
          stop
       end if
 
c  check file type
       call gftype(iounit, type, ier)
       call upper(type)
c  if file type was not given, close and the re-open data file
c           with file type just found, so we can write to file
       if (ftype.eq.' ')  then
           ftype = type
           call closrf(iounit,ier)
           call openrf(iounit, filnam, abrtfl, safefl, ftype,irecl,ier)
c  if file type was given but it was wrong, stop
       elseif (ftype.ne.type) then
                 messg = 'outuwx: incorrect file type for '
            call messag(messg//filnam(:ilen))
                 messg = '        file type for this file is '
            call messag(messg//type)
                 messg = '        file type requested was '
            call messag(messg//ftype)
            stop
       endif
c
c  make sure that all of re(y), im(y), amp(y), and phase(y) are known
       do 300 i = 1, ndata
          if ( ( (yampl(i).eq.zero).and.(yphas(i).eq.zero) ) .and.
     $         ( (yreal(i).ne.zero).or. (yimag(i).ne.zero) ) ) then
            yampl(i) = sqrt( yreal(i)**2 + yimag(i)**2 )
            yphas(i) = atan2( yimag(i), yreal(i) )
             if (i.gt.1) call pijump( yphas(i), yphas(i-1) )
 
          elseif ( (yreal(i).eq.zero).and.(yimag(i).eq.zero)
     $        .and.(yampl(i).ne.zero)   ) then
            yreal(i) = yampl(i) * cos ( yphas(i) )
            yimag(i) = yampl(i) * sin ( yphas(i) )
 
          end if
300    continue
c
c  put data into a single buffer according to data type
c convert ftype to the case of this routine.
c   'case' controls the the case of this routine
       call smcase(ftype, 'case')
c  usually buffer(1) and buffer(2) are xdata(1) and xdata(2) -xdata(1)
       buffer(1) = xdata(1)
       buffer(2) = xdata(2) - xdata(1)
c   xmu: nbuff energy, then nbuff y-values
       if (ftype.eq.'xmu') then
            nbuff    = 2*ndata
            do 400 i = 1, ndata
               buffer(i)         = xdata(i)
               buffer(ndata + i) = yreal(i)
400         continue
c   chi: kmin, deltak, chi(kmin + i*deltak)
       elseif (ftype.eq.'chi') then
            nbuff     = ndata + 2
            do 500 i  = 1, ndata
               buffer(2 + i)     = yreal(i)
500         continue
c   env: kmin, deltak, phase, amplitude pairs (kmin + i*deltak)
       elseif ( (ftype.eq.'env').or.(ftype.eq.'rspep') ) then
            nbuff     = 2* (ndata + 1)
            do 600 i  = 1, ndata
               buffer(2*i+1)     = yphas(i)
               buffer(2*i+2)     = yampl(i)
600         continue
c   rsp: kmin, deltak, real, imaginary pairs (kmin + i*deltak)
       elseif ( (ftype.eq.'rsp').or.(ftype.eq.'rip') ) then
            nbuff     = 2* (ndata + 1)
            do 700 i = 1, ndata
               buffer(2*i+1)     = yreal(i)
               buffer(2*i+2)     = yimag(i)
700         continue
c   other data types not yet supported
       else
            call messag('outuwx: not able to decipher ftype ='//ftype)
            stop
       end if
c 
c  generate skyout for data with hash
       call hash(buffer, nbuff, docout, idoc, skyout)
 
c  check if this record is already in the file,
c    and decide whether or not to write data and
c    documentation for the record to the file
 
       call gnkey(iounit, skyout, nkey, ier)
       if ( (iexist.eq.1).and.(nkey.ne.0) ) then
          skey = ' '
       else
          call putrec(iounit, buffer, nbuff, skyout, 0, ier)
          call putdoc(iounit, docout, idoc,  skyout, 0, ier)
          skey = skyout
       end if
 
       ftypin = ftype
c  close file and leave
       call closrf(iounit, ierr)
       return
c end subroutine outuwx
       end
      subroutine getrea(keywrd,string,value)
      
      character*(*) keywrd, string
      character*72  messg
      integer       j, k, id, ie
      real          value
      logical       isnum

 400  format(bn,f13.0)
 410  format(bn,e19.13)
 420  format(bn,d19.13)

      if (isnum(string)) then
          call lower(string)
          ie = index(string, 'e')
          id = index(string, 'd')
          if ((id.eq.0).and.(ie.eq.0)) then
              read(string,400)value
          elseif (ie.ne.0) then
              read(string,410)value
          elseif (id.ne.0) then
              read(string,420)value
          endif
      else
          j = istrln(string)
          k = istrln(keywrd)
          messg = 'Error reading "'//string(:j)//'" as the value for "'
     $                //keywrd(:k)//'"'
          call messag(messg)
          messg = '"'//string(:j)//'" must be a real number.'
          call messag(messg)
          stop
      endif

      return
c  end subroutine getrea
      end

c**********************************************************************

      subroutine getint(keywrd,string,ivalue)
      
      character*(*) keywrd, string
      character*72  messg
      integer       ivalue, j, k

 400  format(bn,i10)

      read(string,400,iostat=ierr)ivalue
      if (ierr.ne.0) then
          j = istrln(string)
          k = istrln(keywrd)
          messg = 'Error reading '//string(:j)//' as the value for '//
     $                keywrd(:k)
          call messag(messg)
          messg = string(:j)//' must be an integer.'
          call messag(messg)
          stop
      endif

      return
c  end subroutine getint
      end

c**********************************************************************

      subroutine getlog(keywrd,string,lvalue)
      
      character*(*) keywrd, string
      character*72  test*2
      logical       lvalue

      test   = 'ab'
      lvalue = .false.
      call triml(string)
      call case(test,string)
      if ( (string(:1).eq.'t') .or. (string(:1).eq.'y')
     $                         .or. (string(:2).eq.'on') )
     $            lvalue=.true. 

      return
c  end subroutine getlog
      end

c**********************************************************************

      subroutine gettit(keywrd,string,title,ntit)
      
      character*(*) keywrd, string, title
      character*72  messg, toss
      integer       ntit, i, j

      ntit  = ntit + 1
      toss  = string
      call case(keywrd,toss)
      i     = index(toss, keywrd)
      j     = istrln(keywrd)
      title = string(i+j+1:60)
      call triml(title)
      if ( (title(:1).eq.'=') .or. (title(:1).eq.',') ) then
          toss  = title(2:)
          title = toss
          call triml(title)
      endif
      messg = ' title > '//title
      call messag(messg)
      
      return
c  end subroutine gettit
      end

c       subroutine case(test,word)
c c  returns *word* in the same case as *test*
c c  note that this is just the reverse of smcase !
c       character*(*) test, word
c       call smcase (word, test)
c       return
c c  end subroutine case
c       end


c       function  nofx(x,array,npts)
c 
c       implicit integer(i-n)
c       implicit real(a-h,o-z)
c c      implicit double precision(a-h,o-z)
c 
c c
c c   function nofx
c c
c c   purpose
c c     given a value x and an array of values, find the index
c c     corresponding to the array element closest to x
c c
c c   usage
c c     n = nofx(x,array,npts)
c c
c c   parameters
c c     x     - a given value
c c     array - array of values, assumed to be stored in order of
c c             increasing value
c c     npts  - number of elements in array
c c
c c   subroutines and function subprograms required
c c     none
c c
c c   written  8/11/81 by j.m. tranquada
c c
c       dimension  array(npts)
c       imin = 1
c       imax = npts
c       inc = ( imax - imin ) / 2
c    10 continue
c       it  = imin + inc
c       xit = array(it)
c       if ( x .lt. xit ) then
c          imax = it
c       else if ( x .gt. xit ) then
c          imin = it
c       else
c          nofx = it
c          return
c       endif
c       inc = ( imax - imin ) / 2
c       if ( inc .gt. 0 ) go to 10
c       xave = ( array(imin) + array(imin+1) ) / 2.
c       if ( x .lt. xave ) then
c          nofx = imin
c       else
c          nofx = imin + 1
c       endif
c       return
c c end function nofx
c       end
c 

       subroutine openfl(iunit, file, status, iexist, ierr)
c  
c  open a file, 
c   if unit <= 0, the first unused unit number greater than 7 will 
c                be assigned.
c   if status = 'old', the existence of the file is checked.
c   if the file does not exist iexist is set to -1
c   if the file does exist, iexist = iunit.
c   if any errors are encountered, ierr is set to -1.
c
c   note: iunit, iexist, and ierr may be overwritten by this routine
       character*(*)  file, status, stat*10
       integer        iunit, iexist, ierr, iulow
       logical        opend, exist
       data  iulow /7/
c
c make sure there is a unit number
       ierr   = -3
       iexist = 0
       if (iunit.le.0) then
          iunit  = iulow
 10       continue
          inquire(unit=iunit, opened = opend)
          if (opend) then
             if (iunit.gt.20) return
             iunit = iunit + 1
             go to 10
          endif
       end if
c
c if status = 'old', check that the file name exists
       ierr = -2
       stat =  status                          
       call smcase(stat,'a')
       if (stat.eq.'old') then
          iexist = -1
          inquire(file=file, exist = exist)
          if (.not.exist) return
          iexist = iunit
       end if
c 
c open the file
       ierr = -1
       open(unit=iunit, file=file, status=status, err=100)
       ierr = 0
 100   continue
       return
c end  subroutine openfl
       end

      subroutine pijump (ph, old)
c
c     removes jumps of 2*pi in phases
c     ph = current value of phase (may be modified on output, but
c          only by multiples of 2*pi)
c     old = previous value of phase
 
ccc         implicit double precision (a-h, o-z)
      parameter (pi = 3.14159 26535 89793 23846 26433)
      parameter (twopi = 2 * pi)
      dimension xph(3)
 
      xph(1) = ph - old
      jump =  (abs(xph(1))+ pi) / twopi
      xph(2) = xph(1) - jump*twopi
      xph(3) = xph(1) + jump*twopi
 
      xphmin = min (abs(xph(1)), abs(xph(2)), abs(xph(3)))
      do 10  i = 1, 3
         if (abs (xphmin - abs(xph(i))) .le. 0.01)  isave = i
 10   continue
 
      ph = old + xph(isave)
 
      return
c end subroutine pijump
      end



c%%%        integer function nbrstr(string)
c%%%  c
c%%%  c  find a number in a string
c%%%  c  given a string that is known to begin with a digit or sign.
c%%%  c  return the position of the end of the number.
c%%%  c  nbrstr : position of end of number
c%%%        integer   istrln, i, ilen, iback
c%%%        character*(*)  string 
c%%%        character*1    digit*10, plus, minus, d, e, decml, s, sp
c%%%        logical     lexp, ldecml
c%%%        data digit  /'1234567890'/
c%%%        data minus,plus,d,e,decml /'-','+','d','e','.'/
c%%%  c------
c%%%        ldecml = .false.
c%%%        lexp   = .false.
c%%%        ilen   = istrln(string)
c%%%        nbrstr = ilen
c%%%        if (ilen.gt.1) then
c%%%           iback  = 1
c%%%  c find end of number :  digits are always ok.
c%%%  c stop at second d, e, decml, or sign that's not preceded by (d,e)
c%%%           do 200 i = 2, ilen
c%%%              sp = string(i-1:i-1)
c%%%              s  = string(i:i)
c%%%              if (index(digit,s).eq.0) then 
c%%%                 if ( ( (s.ne.plus).and.(s.ne.minus).and.(s.ne.d)
c%%%       $                 .and.(s.ne.e).and.(s.ne.decml))
c%%%       $          .or.( lexp.and.((s.eq.d).or.(s.eq.e)))
c%%%       $          .or.( ldecml.and.(s.eq.decml)) 
c%%%       $          .or.( ( (s.eq.plus).or.(s.eq.minus)).and.
c%%%       $                (sp.ne.d).and.(sp.ne.e)) )     go to 210
c%%%                 lexp   = lexp.or.(s.eq.d).or.(s.eq.e)
c%%%                 ldecml = ldecml.or.(s.eq.decml)
c%%%              end if
c%%%   200     continue 
c%%%           iback = 0
c%%%   210     continue 
c%%%           nbrstr = i - 1 - iback 
c%%%        end if
c%%%        return
c%%%  c  end function nbrstr
c%%%        end
c%%%        subroutine unblnk (string)
c%%%  c
c%%%  c remove blanks from a string
c%%%        integer        i, ilen, j
c%%%        character*(*)  string, str*256
c%%%        ilen = min(256, max(1, istrln(string)))
c%%%        j   = 0
c%%%        str = ' '
c%%%        do 10 i = 1, ilen
c%%%           if (string(i:i).ne.' ') then
c%%%              j = j+1
c%%%              str(j:j) = string(i:i)
c%%%           end if
c%%%   10   continue 
c%%%        string = ' '
c%%%        string = str(1:j)
c%%%        return
c%%%  c end subroutine unblnk
c%%%        end

       subroutine sgetvl(word, xreal, ierr)
c  reads real number "xreal" from input character string "word".
c  returns xreal = 0. and (ierr.ne.0)  if word cannot be a number
      character*(*) word
      real          xreal
      integer       ierr
      logical       number, isnum
      external      isnum
c
      xreal   = 0.
      ierr    = -999
      number  = isnum(word)
      if (number) then
          ierr = 0
          read ( word, '(bn,f30.0)', iostat = ierr)  xreal
      end if    
      return
c end subroutine sgetvl
      end
c---------------------------------------------------------------- 
c%%%        subroutine uneql (string)
c%%%  c replace commas or equal signs with blanks 
c%%%        integer        i, ilen
c%%%        character*(*)  string, equal*1
c%%%        data   equal / '=' /
c%%%        ilen = max(1, istrln(string))
c%%%        do 10 i=1,ilen
c%%%          if (string(i:i).eq.equal) string(i:i) = ' '
c%%%   10   continue 
c%%%        return
c%%%  c end subroutine uneql
c%%%        end
c---------------------------------------------------------------- 
      integer function nxtunt(iunit)

c  this function returns the value of the next unopened unit number equal
c  to or larger than iunit.  it will return neither unit numbers 0, 5, 
c  or 6 nor a negative unit number 

      logical open

      iun = iunit
      if (iun.le.0) iun = 1
 10   continue 
      if ((iun.eq.5).or.(iun.eq.6)) then
          iun = 7
          goto 10
      endif
      inquire (unit=iun, opened=open)
      if (open) then
          iun = iun + 1
          goto 10
      endif
      
      nxtunt = iun
      return 
c  end integer function nxtunt
      end

       logical function isdat(string)
c  tests if string contains numerical data 
c    returns true if the first (up to eight) words in string can
c    all be numbers. requires at least two words, and tests only 
c    the first eight columns
       integer nwords, mwords, i 
       parameter (mwords = 8)
       character*(30)  string*(*), words(mwords), line*(256)
       logical isnum
       external isnum
c       
       isdat = .false.
       do 10 i = 1, mwords 
          words(i) = 'no'
 10    continue
c
       nwords = mwords
       line   = string
       call triml(line)
       call untab(line)
       call bwords(line, nwords, words)
       if (nwords.ge.2) then
          isdat = .true.
          do 50 i = 1, nwords
             if (.not. ( isnum( words(i) ) ) ) isdat = .false.
 50       continue
       end if
c
       return
       end


       subroutine testrf(flnam, irecl, flform, ier)
c
c   test whether a data file can be interpreted as  uwxafs binary 
c   data file or  ascii column data file.
c
c   uwxafs binary files use direct access binary files 
c   with word size irecl, which is a machine dependent parameter
c
c ier = -1 : file not found
c ier = -2 : broken uwxafs file?
c ier = -3 : not uwxafs file, but can't find data. 
c ier = -4 : looks like ascii, saw line  of minus signs, 
c             but 2nd following line doesn't have data
c
c   copright 1994 university of washington   matt newville
c -----------------------------------------------------
      integer   i, irecl, iunit
      parameter (mwords = 3)
      character*(*) flnam, flform, line*80
      integer*2    indx(4)
      logical    exist, opend, isdat, prevdt, lisdat
      external  isdat
c -----------------------------------------------------
      flform = 'none'
      ier    = -1
      iunit  = 7
 10   continue
      inquire(unit=iunit, opened = opend)
      if (opend) then 
         if (iunit.gt.20) return
         iunit = iunit + 1   
         go to 10
      endif
      inquire(file = flnam, exist = exist)
      if (.not.exist) return
      ier    = -2
c -----------------------------------------------------
c try reading file as a uwxafs binary file
c     which have patriotic magic numbers embedded in them
      indx(3) = 0
      indx(4) = 0
      open(iunit, file= flnam, recl = irecl, err = 20, 
     $      access = 'direct', status = 'old' )
 20   continue
      read(iunit, rec=1, err = 25) (indx(i), i=1,4)
 25   continue
      if ((indx(3).eq.1776).and.(indx(4).eq.704)) then
         flform = 'uwxafs'
         ier  = 0
         go to 900
      end if
c -----------------------------------------------------
c try to read file as ascii data file
      close(iunit)
      open(iunit, file=flnam, status='old')
      prevdt = .false.
 200  continue
         ier  = -3
         read(iunit, '(a)', end = 900, err = 900) line
         call triml (line)
         if (line(3:6) .eq. '----') then 
            ier = -4
            read(iunit, '(a)', end = 900, err = 900) line
            read(iunit, '(a)', end = 900, err = 900) line
            lisdat = isdat(line)
            if (lisdat ) then
               flform = 'ascii'
               ier = 0
            end if
            go to 900
         end if
c if two lines in a row have all words being numbers, it is  data
         lisdat = isdat(line)
         if (lisdat.and.prevdt)  then 
            flform = 'ascii'
            ier = 0
            go to 900
         end if
         prevdt = lisdat
         go to 200
c---------------------
 900  continue
      close(iunit)
      return
c end subroutine testrf
      end
      subroutine mcfluo (cent, edge, mel, elem, perc,
     $           latwei, mutf, mub, muc)         

c   This subroutine will calculate the three total cross sections 
c   mutf, mub and muc needed for calculating self absorption corrections
c   to fluorescence data. 
c   input variables:
c      cent - central atom symbol
c      edge - edge of interest symbol: k, l1, l2, l3
c      mel  - length of array with atoms symbols and relative amounts
c      elem - array of characters of length mel 
c      perc - array of reals of length mel
c      latwei - logical flag = .true. if relative amounts of elements given 
c               in units of atomic weigth. The default is at. %
c   output variables:
c      mutf - Total cross section at the fluorescent energy
c      mub - background absorption xsection 
c      muc - central atom xsection of edge of interest

      parameter(zero=0.e0, imel=20)
      parameter( weimin = 1.e-5, alpmin = 0.001, ediff = 0.01 )
      dimension    perc(imel), abxsec(imel), blxsec(imel), 
     $             flxsec(imel), xsec(10), energy(9),
     $             weigat(imel)   
      dimension    izel(imel)
      character*2  elem(imel), test, cent, edge
      character*1  unit
      character*78 errom1
      real         mutf, mub, muc, percen, ecent, eabove, ebelow
      integer      incen, izcent
      logical      latwei, erf

c     Initialize some variables used in this subroutine
      print*, ' Mcfluo calculating corrections ...'
      test ='ab'
      do 5 i=1, imel
         abxsec(i)=zero
         blxsec(i)=zero
         flxsec(i)=zero
         weigat(i)=zero
         izel(i)=0
 5    continue
      incen=0

      izcent=is2z(cent)
      if (izcent.eq.0)  goto 80
      do 10 i= 1,mel
         call case (test, elem(i))
         izel(i)=is2z(elem(i))
         if (izel(i).eq.izcent) then
            incen=i
            percen=perc(i)
         endif   
         if (izel(i).eq.0) goto 81
 10   continue

      if (incen.eq.0) goto 79

      call smcase(edge, test)
      if ((edge(1:2) .ne. 'k ').and.(edge(1:2) .ne. 'l1')
     $     .and.(edge(1:2).ne. 'l2').and.(edge(1:2).ne. 'l3'))
     $        goto 78

      ecent = s2e(cent, edge)
      eabove = ecent + ediff 
      ebelow = ecent - ediff

c     initialize mucal variables
      erf = .false.

      do 11 i=1,9
         energy(i)=zero
         xsec(i)  =zero
 11   continue
      xsec(10)=zero

      do 30 i=1, mel
      unit = 'b'
c above
         call mucal(eabove, elem(i), izel(i), unit, xsec, energy, 
     $              erf, ier)
         abxsec(i)=xsec(4)
         weigat(i)=xsec(7)
c below 
         call mucal(ebelow, elem(i), izel(i), unit, xsec, energy, 
     $              erf, ier)
         blxsec(i)=xsec(4)
         if ((i.eq.incen).and.(edge(1:1).eq.'k')) alphfl = energy(6)
         if ((i.eq.incen).and.(edge(1:1).eq.'l')) alphfl = energy(8)
 30   continue
      if (alphfl.le.alpmin)  go to 90
c      print*, alphfl
      do 40 i=1, mel
         unit = 'b'
         call mucal (alphfl, elem(i), izel(i), unit, xsec, energy, 
     $               erf, ier)
         flxsec(i)=xsec(4)
 40   continue

c     If percent is given in atomic weight units divide by atomic weight
      
      if (latwei) then
         do 41 i=1, mel
           weigat(i) = max( weimin,  weigat(i) )
           perc(i)=perc(i)/(weigat(i))
 41      continue
      endif

c     Calculate mub, the total "background" xsection. 
c     All atoms at eabove except central atom at ebelow
      
c     calculate mutf, total xsection at fluorescence energy
      mub  = zero
      mutf = zero
      do 50 i=1, mel
         mub  = mub  + abxsec(i)*perc(i)
         mutf = mutf + flxsec(i)*perc(i)
 50   continue

c     correct mub for having taken central atom above edge

      mub= mub - abxsec(incen)*perc(incen) + blxsec(incen)*perc(incen)

c     calculate muc, edge of interest xsection

      muc =  (abxsec(incen)-blxsec(incen))*perc(incen)

 70   continue
      return

 78   continue
      errom1 = '** The edge of interest has to be either '//
     $         ' K, L1, L2 or L3 **'
      goto 82

 79   continue
      errom1 = '** Your central atom must be one of the '//
     $         ' elements in the atom list **'
      goto 82
      
 80   continue
      errom1 = '** You did not specify a central atom or'//
     $         ' its symbol is unrecognized **'
      goto 82
 81   continue
      errom1 = '** The atom symbol "'//elem(i)//'" is not an'//
     $         ' element **'
 82   continue
      mi = istrln(errom1)
      call messag (errom1(1:mi))
      call messag ('*** Fluo is stopping ***')
      stop
 90   continue
      call messag ('** no fluorescence energy found **')
      call messag ('*** Fluo is stopping ***')
      stop
   
c     end subroutine mcfluo
      end

c****************************************************************
      subroutine mucal(en,sym,iz,unit,xsec,energy,erf,ier)

c---------------------------------------------------------------------
c standardized, made double precision, stylized and adapted for atom
c by b ravel, september 1992
c---------------------------------------------------------------------
c      this is a routine to calculate x-sections using data from 
c      the may 1969 edition of mcmaster. 
c      written by pathikrit bandyopadhyay at the university of 
c      notre dame.
c
c      this program has data for all the elements from z=1 to 94
c      with the following exceptions:
c      84  po \
c      85  at |
c      87  fr |
c      88  ra  > mcmaster does not publish data for these elements.
c      89  ac |
c      91  pa |
c      93  np /
c
c****** please correct data errors as they are found.  thanx. ******
c
c      input:  
c        en:   energy at which to calculate the x-section, en<0 means 
c              to skip calculation of absorption
c   **   symb: name of material  
c        unit: units to be used. 'c' for cm**2/gm, 'b' for barns/atom   
c   **   iz:   z number of element
c        erf:  logical print flag (true = print run-time error messages)
c   **either or both supplied on input, both returned on output
c
c      output:
c        xsec(1)   = photoelectric x-section
c        xsec(2)   = coherent x-section
c        xsec(3)   = incoherent x-section
c        xsec(4)   = total x-section
c        xsec(5)   = conversion factor
c        xsec(6)   = absorption coefficient
c        xsec(7)   = atomic weight
c        xsec(8)   = density
c        xsec(9)   = l2-edge jump
c        xsec(10)  = l3-edge jump
c        energy(1) = k-edge energy
c        energy(2) = l1-edge energy
c        energy(3) = l2-edge energy
c        energy(4) = l3-edge energy
c        energy(5) = m-edge energy
c        energy(6) = k-alpha1
c        energy(7) = k-beta1
c        energy(8) = l-alpha1
c        energy(9) = l-beta1
c        ier       = error code
c
c      error codes:
c        ier=1: energy input is zero
c        ier=2: name does not match z
c        ier=3: no documentation for given element (z<94)
c        ier=4: no documentation for given element (z>94)
c        ier=5: l-edge calculation may be wrong for z<30 as mcmaster
c               uses l1 only.
c        ier=6: energy at the middle of edge
c        ier=7: no name or z supplied
c
c--------------------------------------------------------------------
      implicit integer(i-n)
      implicit real(a-h,o-z)
c      implicit double precision(a-h,o-z)

      parameter(nelem=94, eps=.001)

c      double precision ka(nelem),kb(nelem),la(nelem),lb(nelem)
c      double precision l2(nelem),l3(nelem),lj3(nelem),lj2
      real             ka(nelem),kb(nelem),la(nelem),lb(nelem)
      real             l2(nelem),l3(nelem),lj3(nelem),lj2
      logical          erf
      character*1      unit
      character*2      name(nelem),symb,sym,test
      dimension        ek(nelem),el(nelem),em(nelem),den(nelem),
     $                 atwt(nelem),cf(nelem)
      dimension        ak(0:3,nelem),al(0:3,nelem),am(0:3,nelem),
     $                 an(0:3,nelem)
      dimension        cih(0:3,nelem),coh(0:3,nelem)
      dimension        xsec(10),energy(9)

c###the following is the symbol, edge energies, mcmaster
c   coefficients, flouescent energies, density, conversion 
c   constant, and atomic weight of each supported element 
c   up to plutonium

      data name( 1),  ek( 1),el( 1)/
     $  'h ',  0.140000e-01,  0.000000e+00/
      data ak(0, 1),ak(1, 1),ak(2, 1),ak(3, 1)/
     $   0.244964e+01, -0.334953e+01, -0.471370e-01,  0.709962e-02/
      data al(0, 1),al(1, 1),al(2, 1),al(3, 1)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data am(0, 1),am(1, 1),am(2, 1),am(3, 1)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data an(0, 1),an(1, 1),an(2, 1),an(3, 1)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den( 1),cf( 1),atwt(1)/1.008, 0.167400e+01,1.008/
      data coh(0, 1),coh(1, 1),coh(2, 1),coh(3, 1)/
     $  -0.119075e+00, -0.937086e+00, -0.200538e+00,  0.106587e-01/
      data cih(0, 1),cih(1, 1),cih(2, 1),cih(3, 1)/
     $  -0.215772e+01,  0.132685e+01, -0.305620e+00,  0.185025e-01/
      data name( 2),  ek( 2),el( 2)/
     $  'he',  0.250000e-01,  0.000000e+00/
      data ak(0, 2),ak(1, 2),ak(2, 2),ak(3, 2)/
     $   0.606488e+01, -0.329055e+01, -0.107256e+00,  0.144465e-01/
      data al(0, 2),al(1, 2),al(2, 2),al(3, 2)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data am(0, 2),am(1, 2),am(2, 2),am(3, 2)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data an(0, 2),an(1, 2),an(2, 2),an(3, 2)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den( 2),cf( 2),atwt(2)/1.785e-04, 0.664700e+01,4.003/
      data coh(0, 2),coh(1, 2),coh(2, 2),coh(3, 2)/
     $   0.104768e+01, -0.851805e-01, -0.403527e+00, 0.269398e-01/
      data cih(0, 2),cih(1, 2),cih(2, 2),cih(3, 2)/
     $  -0.256357e+01,  0.202536e+01, -0.448710e+00,  0.279691e-01/
      data name( 3),  ek( 3),el( 3)/
     $  'li',  0.550000e-01,  0.000000e+00/
      data ak(0, 3),ak(1, 3),ak(2, 3),ak(3, 3)/
     $   0.775370e+01, -0.281801e+01, -0.241378e+00,  0.262542e-01/
      data al(0, 3),al(1, 3),al(2, 3),al(3, 3)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data am(0, 3),am(1, 3),am(2, 3),am(3, 3)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data an(0, 3),an(1, 3),an(2, 3),an(3, 3)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den( 3),cf( 3),atwt(3)/0.534, 0.115200e+02,6.940/
      data coh(0, 3),coh(1, 3),coh(2, 3),coh(3, 3)/
     $   0.134366e+01,  0.181557e+00, -0.423981e+00,  0.266190e-01/
      data cih(0, 3),cih(1, 3),cih(2, 3),cih(3, 3)/
     $  -0.108740e+01,  0.103368e+01, -0.190377e+00,  0.779955e-02/
      data name( 4),  ek( 4),el( 4)/
     $  'be',  0.112000e+00,  0.000000e+00/
      data ak(0, 4),ak(1, 4),ak(2, 4),ak(3, 4)/
     $   0.904511e+01, -0.283487e+01, -0.210021e+00,  0.229526e-01/
      data al(0, 4),al(1, 4),al(2, 4),al(3, 4)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data am(0, 4),am(1, 4),am(2, 4),am(3, 4)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data an(0, 4),an(1, 4),an(2, 4),an(3, 4)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den( 4),cf( 4),atwt(4)/1.848, 0.149600e+02,9.012/
      data coh(0, 4),coh(1, 4),coh(2, 4),coh(3, 4)/
     $   0.200860e+01, -0.461920e-01, -0.337018e+00,  0.186939e-01/
      data cih(0, 4),cih(1, 4),cih(2, 4),cih(3, 4)/
     $  -0.690079e+00,  0.946448e+00, -0.171142e+00,  0.651413e-02/
      data name( 5),  ek( 5),el( 5)/
     $  'b ',  0.188000e+00,  0.000000e+00/
      data ak(0, 5),ak(1, 5),ak(2, 5),ak(3, 5)/
     $   0.995057e+01, -0.274174e+01, -0.215138e+00,  0.227845e-01/
      data al(0, 5),al(1, 5),al(2, 5),al(3, 5)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data am(0, 5),am(1, 5),am(2, 5),am(3, 5)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data an(0, 5),an(1, 5),an(2, 5),an(3, 5)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den( 5),cf( 5),atwt(5)/2.340, 0.179500e+02,10.811/
      data coh(0, 5),coh(1, 5),coh(2, 5),coh(3, 5)/
     $   0.262862e+01, -0.207916e+00, -0.286283e+00,  0.144966e-01/
      data cih(0, 5),cih(1, 5),cih(2, 5),cih(3, 5)/
     $  -0.791177e+00,  0.121611e+01, -0.239087e+00,  0.117686e-01/
      data name( 6),  ek( 6),el( 6)/
     $  'c ',  0.284000e+00,  0.000000e+00/
      data ak(0, 6),ak(1, 6),ak(2, 6),ak(3, 6)/
     $   0.106879e+02, -0.271400e+01, -0.200530e+00,  0.207248e-01/
      data al(0, 6),al(1, 6),al(2, 6),al(3, 6)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data am(0, 6),am(1, 6),am(2, 6),am(3, 6)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data an(0, 6),an(1, 6),an(2, 6),an(3, 6)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den( 6),cf( 6),atwt(6)/2.250, 0.199400e+02,12.01/
      data coh(0, 6),coh(1, 6),coh(2, 6),coh(3, 6)/
     $   0.310861e+01, -0.260580e+00, -0.271974e+00,  0.135181e-01/
      data cih(0, 6),cih(1, 6),cih(2, 6),cih(3, 6)/
     $  -0.982878e+00,  0.146693e+01, -0.293743e+00,  0.156005e-01/
      data name( 7),  ek( 7),el( 7)/
     $  'n ',  0.402000e+00,  0.000000e+00/
      data ak(0, 7),ak(1, 7),ak(2, 7),ak(3, 7)/
     $   0.112765e+02, -0.265400e+01, -0.200445e+00,  0.200765e-01/
      data al(0, 7),al(1, 7),al(2, 7),al(3, 7)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den( 7),cf( 7),atwt(7)/1.250e-03, 0.232600e+02,14.008/
      data coh(0, 7),coh(1, 7),coh(2, 7),coh(3, 7)/
     $   0.347760e+01, -0.215762e+00, -0.288874e+00,  0.115131e-01/
      data cih(0, 7),cih(1, 7),cih(2, 7),cih(3, 7)/
     $  -0.123693e+01,  0.174510e+01, -0.354660e+00,  0.198705e-01/
      data name( 8),  ek( 8),el( 8)/
     $  'o ',  0.537000e+00,  0.000000e+00/
      data ak(0, 8),ak(1, 8),ak(2, 8),ak(3, 8)/
     $   0.117130e+02, -0.257229e+01, -0.205893e+00,  0.199244e-01/
      data al(0, 8),al(1, 8),al(2, 8),al(3, 8)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den( 8),cf( 8),atwt(8)/1.429e-03, 0.265700e+02,16.0/
      data coh(0, 8),coh(1, 8),coh(2, 8),coh(3, 8)/
     $   0.377239e+01, -0.148539e+00, -0.307124e+00,  0.167303e-01/
      data cih(0, 8),cih(1, 8),cih(2, 8),cih(3, 8)/
     $  -0.173679e+01,  0.217686e+01, -0.449050e+00,  0.264733e-01/
      data name( 9),  ek( 9),el( 9)/
     $  'f ',  0.686000e+00,  0.000000e+00/
      data ak(0, 9),ak(1, 9),ak(2, 9),ak(3, 9)/
     $   0.120963e+02, -0.244148e+01, -0.234461e+00,  0.219954e-01/
      data al(0, 9),al(1, 9),al(2, 9),al(3, 9)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den( 9),cf( 9),atwt(9)/1.108, 0.315500e+02,19.0/
      data coh(0, 9),coh(1, 9),coh(2, 9),coh(3, 9)/
     $   0.400716e+01, -0.560908e-01, -0.332017e+00,  0.187934e-01/
      data cih(0, 9),cih(1, 9),cih(2, 9),cih(3, 9)/
     $  -0.187570e+01,  0.232016e+01, -0.475412e+00,  0.280680e-01/
      data name(10),  ek(10),el(10)/
     $  'ne',  0.867000e+00,  0.000000e+00/
      data ak(0,10),ak(1,10),ak(2,10),ak(3,10)/
     $   0.124485e+02, -0.245819e+01, -0.212591e+00,  0.196489e-01/
      data al(0,10),al(1,10),al(2,10),al(3,10)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den(10),cf(10),atwt(10)/9.000e-04, 0.335100e+02,20.183/
      data coh(0,10),coh(1,10),coh(2,10),coh(3,10)/
     $   0.420151e+01,  0.416247e-01, -0.356754e+00,  0.207585e-01/
      data cih(0,10),cih(1,10),cih(2,10),cih(3,10)/
     $  -0.175510e+01,  0.224226e+01, -0.447640e+00,  0.255801e-01/
      data name(11),  ek(11),el(11)/
     $  'na',  0.107200e+01,  0.000000e+00/
      data ak(0,11),ak(1,11),ak(2,11),ak(3,11)/
     $   0.126777e+02, -0.224521e+01, -0.274873e+00,  0.250270e-01/
      data al(0,11),al(1,11),al(2,11),al(3,11)/
     $   0.102355e+02, -0.255905e+01, -0.119524e+00,  0.000000e+00/
      data den(11),cf(11),atwt(11)/0.970, 0.381900e+02,22.997/
      data coh(0,11),coh(1,11),coh(2,11),coh(3,11)/
     $   0.426374e+01,  0.134662e+00, -0.370080e+00,  0.214467e-01/
      data cih(0,11),cih(1,11),cih(2,11),cih(3,11)/
     $  -0.967717e+00,  0.161794e+01, -0.287191e+00,  0.131526e-01/
      data ka(11),kb(11),la(11),lb(11)/
     $   0.104100e+01,  0.106700e+01,  0.000000e+00,  0.000000e+00/
      data name(12),  ek(12),el(12)/
     $  'mg',  0.130500e+01,  0.630000e-01/
      data ak(0,12),ak(1,12),ak(2,12),ak(3,12)/
     $   0.128793e+02, -0.212574e+01, -0.299392e+00,  0.267643e-01/
      data al(0,12),al(1,12),al(2,12),al(3,12)/
     $   0.105973e+02, -0.289818e+01,  0.234506e+00,  0.000000e+00/
      data den(12),cf(12),atwt(12)/1.740, 0.403800e+02,24.320/
      data coh(0,12),coh(1,12),coh(2,12),coh(3,12)/
     $   0.439404e+01,  0.137858e+00, -0.359540e+00,  0.202380e-01/
      data cih(0,12),cih(1,12),cih(2,12),cih(3,12)/
     $  -0.571611e+00,  0.135498e+01, -0.222491e+00,  0.830141e-02/
      data ka(12),kb(12),la(12),lb(12)/
     $   0.125400e+01,  0.129700e+01,  0.000000e+00,  0.000000e+00/
      data name(13),  ek(13),el(13)/
     $  'al',  0.156000e+01,  0.870000e-01/
      data ak(0,13),ak(1,13),ak(2,13),ak(3,13)/
     $   0.131738e+02, -0.218203e+01, -0.258960e+00,  0.222840e-01/
      data al(0,13),al(1,13),al(2,13),al(3,13)/
     $   0.108711e+02, -0.277860e+01,  0.175853e+00,  0.000000e+00/
      data den(13),cf(13),atwt(13)/2.720, 0.447800e+02,26.970/
      data coh(0,13),coh(1,13),coh(2,13),coh(3,13)/
     $   0.415995e+01,  0.140549e+00, -0.352441e+00,  0.193692e-01/
      data cih(0,13),cih(1,13),cih(2,13),cih(3,13)/
     $  -0.439322e+00,  0.130867e+01, -0.211648e+00,  0.754210e-02/
      data ka(13),kb(13),la(13),lb(13)/
     $   0.148700e+01,  0.155300e+01,  0.000000e+00,  0.000000e+00/
      data name(14),  ek(14),el(14)/
     $  'si',  0.183900e+01,  0.118000e+00/
      data ak(0,14),ak(1,14),ak(2,14),ak(3,14)/
     $   0.132682e+02, -0.198174e+01, -0.316950e+00,  0.273928e-01/
      data al(0,14),al(1,14),al(2,14),al(3,14)/
     $   0.112237e+02, -0.273694e+01,  0.127557e+00,  0.000000e+00/
      data den(14),cf(14),atwt(14)/2.330, 0.466300e+02,28.086/
      data coh(0,14),coh(1,14),coh(2,14),coh(3,14)/
     $   0.464678e+01,  0.162780e+00, -0.358563e+00,  0.196926e-01/
      data cih(0,14),cih(1,14),cih(2,14),cih(3,14)/
     $  -0.414971e+00,  0.134868e+01, -0.222315e+00,  0.841959e-02/
      data ka(14),kb(14),la(14),lb(14)/
     $   0.174000e+01,  0.183200e+01,  0.000000e+00,  0.000000e+00/
      data name(15),  ek(15),el(15)/
     $  'p ',  0.214900e+01,  0.153000e+00/
      data ak(0,15),ak(1,15),ak(2,15),ak(3,15)/
     $   0.133735e+02, -0.186342e+01, -0.339440e+00,  0.288858e-01/
      data al(0,15),al(1,15),al(2,15),al(3,15)/
     $   0.115508e+02, -0.292200e+01,  0.254262e+00,  0.000000e+00/
      data den(15),cf(15),atwt(15)/1.820, 0.514300e+02,30.975/
      data coh(0,15),coh(1,15),coh(2,15),coh(3,15)/
     $   0.478525e+01,  0.168708e+00, -0.360383e+00,  0.197155e-01/
      data cih(0,15),cih(1,15),cih(2,15),cih(3,15)/
     $  -0.476903e+00,  0.146032e+01, -0.251331e+00,  0.107202e-01/
      data ka(15),kb(15),la(15),lb(15)/
     $   0.201500e+01,  0.213600e+01,  0.000000e+00,  0.000000e+00/
      data name(16),  ek(16),el(16)/
     $  's ',  0.247200e+01,  0.193000e+00/
      data ak(0,16),ak(1,16),ak(2,16),ak(3,16)/
     $   0.137394e+02, -0.204786e+01, -0.273259e+00,  0.229976e-01/
      data al(0,16),al(1,16),al(2,16),al(3,16)/
     $   0.118181e+02, -0.264618e+01, -0.968049e-01,  0.000000e+00/
      data am(0,16),am(1,16),am(2,16),am(3,16)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data an(0,16),an(1,16),an(2,16),an(3,16)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den(16),cf(16)/2.000, 0.532400e+02/
      data em(16),atwt(16)/ 0.170000e-01,32.066/
      data coh(0,16),coh(1,16),coh(2,16),coh(3,16)/
     $   0.492707e+01,  0.165746e+00, -0.359424e+00,  0.195505e-01/
      data cih(0,16),cih(1,16),cih(2,16),cih(3,16)/
     $  -0.656419e+00,  0.165408e+01, -0.298623e+00,  0.142979e-01/
      data ka(16),kb(16),la(16),lb(16)/
     $   0.230800e+01,  0.246400e+01,  0.000000e+00,  0.000000e+00/
      data name(17),  ek(17),el(17)/
     $  'cl',  0.282200e+01,  0.238000e+00/
      data ak(0,17),ak(1,17),ak(2,17),ak(3,17)/
     $   0.136188e+02, -0.171937e+01, -0.354154e+00,  0.290841e-01/
      data al(0,17),al(1,17),al(2,17),al(3,17)/
     $   0.120031e+02, -0.241694e+01, -0.240897e+00,  0.000000e+00/
      data am(0,17),am(1,17),am(2,17),am(3,17)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data an(0,17),an(1,17),an(2,17),an(3,17)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den(17),cf(17)/1.560, 0.588700e+02/
      data em(17),atwt(17)/ 0.170000e-01,35.457/
      data coh(0,17),coh(1,17),coh(2,17),coh(3,17)/
     $   0.507222e+01,  0.149127e+00, -0.352858e+00,  0.189439e-01/
      data cih(0,17),cih(1,17),cih(2,17),cih(3,17)/
     $  -0.718627e+00,  0.174294e+01, -0.319429e+00,  0.158429e-01/
      data ka(17),kb(17),la(17),lb(17)/
     $   0.262200e+01,  0.281500e+01,  0.000000e+00,  0.000000e+00/
      data name(18),  ek(18),el(18)/
     $  'ar',  0.320200e+01,  0.287000e+00/
      data ak(0,18),ak(1,18),ak(2,18),ak(3,18)/
     $   0.139491e+02, -0.182276e+01, -0.328827e+00,  0.274382e-01/
      data al(0,18),al(1,18),al(2,18),al(3,18)/
     $   0.122960e+02, -0.263279e+01, -0.736600e-01,  0.000000e+00/
      data am(0,18),am(1,18),am(2,18),am(3,18)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data an(0,18),an(1,18),an(2,18),an(3,18)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den(18),cf(18)/1.784e-03, 0.663200e+02/
      data em(18),atwt(18)/ 0.270000e-01,39.944/
      data coh(0,18),coh(1,18),coh(2,18),coh(3,18)/
     $   0.521079e+01,  0.135618e+00, -0.347214e+00,  0.184333e-01/
      data cih(0,18),cih(1,18),cih(2,18),cih(3,18)/
     $  -0.682105e+00,  0.174279e+01, -0.317646e+00,  0.156467e-01/
      data ka(18),kb(18),la(18),lb(18)/
     $   0.295700e+01,  0.319200e+01,  0.000000e+00,  0.000000e+00/
      data name(19),  ek(19),el(19)/
     $  'k ',  0.360700e+01,  0.341000e+00/
      data ak(0,19),ak(1,19),ak(2,19),ak(3,19)/
     $   0.137976e+02, -0.154015e+01, -0.394528e+00,  0.323561e-01/
      data al(0,19),al(1,19),al(2,19),al(3,19)/
     $   0.124878e+02, -0.253656e+01, -0.104892e+00,  0.000000e+00/
      data am(0,19),am(1,19),am(2,19),am(3,19)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data an(0,19),an(1,19),an(2,19),an(3,19)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den(19),cf(19)/0.862, 0.649300e+02/
      data em(19),atwt(19)/ 0.340000e-01,39.102/
      data coh(0,19),coh(1,19),coh(2,19),coh(3,19)/
     $   0.525587e+01,  0.188040e+00, -0.359623e+00,  0.193085e-01/
      data cih(0,19),cih(1,19),cih(2,19),cih(3,19)/
     $  -0.344007e+00,  0.149236e+01, -0.254135e+00,  0.107684e-01/
      data ka(19),kb(19),la(19),lb(19)/
     $   0.331300e+01,  0.358900e+01,  0.000000e+00,  0.000000e+00/
      data name(20),  ek(20),el(20)/
     $  'ca',  0.403800e+01,  0.400000e+00/
      data ak(0,20),ak(1,20),ak(2,20),ak(3,20)/
     $   0.142950e+02, -0.188644e+01, -0.283647e+00,  0.226263e-01/
      data al(0,20),al(1,20),al(2,20),al(3,20)/
     $   0.127044e+02, -0.255011e+01, -0.943195e-01,  0.000000e+00/
      data am(0,20),am(1,20),am(2,20),am(3,20)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data an(0,20),an(1,20),an(2,20),an(3,20)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den(20),cf(20)/1.550, 0.665500e+02/
      data em(20),atwt(20)/ 0.440000e-01,40.080/
      data coh(0,20),coh(1,20),coh(2,20),coh(3,20)/
     $   0.532375e+01,  0.206685e+00, -0.361664e+00,  0.193328e-01/
      data cih(0,20),cih(1,20),cih(2,20),cih(3,20)/
     $  -0.982420e-01,  0.132829e+01, -0.213747e+00,  0.773065e-02/
      data ka(20),kb(20),la(20),lb(20)/
     $   0.369100e+01,  0.401200e+01,  0.000000e+00,  0.000000e+00/
      data name(21),  ek(21),el(21)/
     $  'sc',  0.449300e+01,  0.463000e+00/
      data ak(0,21),ak(1,21),ak(2,21),ak(3,21)/
     $   0.139664e+02, -0.140872e+01, -0.414365e+00,  0.334355e-01/
      data al(0,21),al(1,21),al(2,21),al(3,21)/
     $   0.128949e+02, -0.240609e+01, -0.177791e+00,  0.000000e+00/
      data am(0,21),am(1,21),am(2,21),am(3,21)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data an(0,21),an(1,21),an(2,21),an(3,21)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den(21), cf(21)/ 0.299200e+01, 0.746500e+02/
      data em(21),atwt(21)/ 0.540000e-01,44.960/
      data coh(0,21),coh(1,21),coh(2,21),coh(3,21)/
     $   0.543942e+01,  0.200174e+00, -0.359064e+00,  0.191027e-01/
      data cih(0,21),cih(1,21),cih(2,21),cih(3,21)/
     $  -0.159831e+00,  0.139055e+01, -0.225849e+00,  0.851954e-02/
      data ka(21),kb(21),la(21),lb(21)/
     $   0.409000e+01,  0.446000e+01,  0.000000e+00,  0.000000e+00/
      data name(22),  ek(22),el(22)/
     $  'ti',  0.496500e+01,  0.531000e+00/
      data ak(0,22),ak(1,22),ak(2,22),ak(3,22)/
     $   0.143506e+02, -0.166322e+01, -0.331539e+00,  0.262065e-01/
      data al(0,22),al(1,22),al(2,22),al(3,22)/
     $   0.131075e+02, -0.253576e+01, -0.957177e-01,  0.000000e+00/
      data am(0,22),am(1,22),am(2,22),am(3,22)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data an(0,22),an(1,22),an(2,22),an(3,22)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den(22), cf(22)/ 0.454000e+01, 0.795300e+02/
      data em(22),atwt(22)/ 0.590000e-01,47.9/
      data coh(0,22),coh(1,22),coh(2,22),coh(3,22)/
     $   0.555039e+01,  0.197697e+00, -0.357694e+00,  0.189866e-01/
      data cih(0,22),cih(1,22),cih(2,22),cih(3,22)/
     $  -0.230573e+00,  0.145848e+01, -0.239160e+00,  0.938528e-02/
      data ka(22),kb(22),la(22),lb(22)/
     $   0.451000e+01,  0.493100e+01,  0.000000e+00,  0.000000e+00/
      data name(23),  ek(23),el(23)/
     $  'v ',  0.546500e+01,  0.604000e+00/
      data ak(0,23),ak(1,23),ak(2,23),ak(3,23)/
     $   0.147601e+02, -0.188867e+01, -0.271861e+00,  0.215792e-01/
      data al(0,23),al(1,23),al(2,23),al(3,23)/
     $   0.132514e+02, -0.249765e+01, -0.106383e+00,  0.000000e+00/
      data am(0,23),am(1,23),am(2,23),am(3,23)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data an(0,23),an(1,23),an(2,23),an(3,23)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den(23), cf(23)/ 0.611000e+01, 0.845900e+02/
      data em(23),atwt(23)/ 0.670000e-01,50.942/
      data coh(0,23),coh(1,23),coh(2,23),coh(3,23)/
     $   0.565514e+01,  0.199533e+00, -0.357487e+00,  0.189691e-01/
      data cih(0,23),cih(1,23),cih(2,23),cih(3,23)/
     $  -0.308103e+00,  0.152879e+01, -0.252768e+00,  0.102571e-01/
      data ka(23),kb(23),la(23),lb(23)/
     $   0.495200e+01,  0.542700e+01,  0.000000e+00,  0.000000e+00/
      data name(24),  ek(24),el(24)/
     $  'cr',  0.598900e+01,  0.682000e+00/
      data ak(0,24),ak(1,24),ak(2,24),ak(3,24)/
     $   0.148019e+02, -0.182430e+01, -0.279116e+00,  0.217324e-01/
      data al(0,24),al(1,24),al(2,24),al(3,24)/
     $   0.134236e+02, -0.251532e+01, -0.101999e+00,  0.000000e+00/
      data am(0,24),am(1,24),am(2,24),am(3,24)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data an(0,24),an(1,24),an(2,24),an(3,24)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den(24), cf(24)/ 0.719000e+01, 0.863400e+02/
      data em(24),atwt(24)/ 0.740000e-01,51.996/
      data coh(0,24),coh(1,24),coh(2,24),coh(3,24)/
     $   0.577399e+01,  0.203858e+00, -0.359699e+00,  0.192225e-01/
      data cih(0,24),cih(1,24),cih(2,24),cih(3,24)/
     $  -0.387641e+00,  0.159727e+01, -0.266240e+00,  0.111523e-01/
      data ka(24),kb(24),la(24),lb(24)/
     $   0.541400e+01,  0.594600e+01,  0.000000e+00,  0.000000e+00/
      data name(25),  ek(25),el(25)/
     $  'mn',  0.654000e+01,  0.754000e+00/
      data ak(0,25),ak(1,25),ak(2,25),ak(3,25)/
     $   0.148965e+02, -0.179872e+01, -0.283664e+00,  0.222095e-01/
      data al(0,25),al(1,25),al(2,25),al(3,25)/
     $   0.135761e+02, -0.249761e+01, -0.105493e+00,  0.000000e+00/
      data am(0,25),am(1,25),am(2,25),am(3,25)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data an(0,25),an(1,25),an(2,25),an(3,25)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den(25), cf(25)/ 0.742000e+01, 0.912200e+02/
      data em(25),atwt(25)/ 0.840000e-01,54.940/
      data coh(0,25),coh(1,25),coh(2,25),coh(3,25)/
     $   0.584604e+01,  0.213814e+00, -0.359718e+00,  0.191459e-01/
      data cih(0,25),cih(1,25),cih(2,25),cih(3,25)/
     $  -0.247059e+00,  0.149722e+01, -0.238781e+00,  0.893208e-02/
      data ka(25),kb(25),la(25),lb(25)/
     $   0.589800e+01,  0.649000e+01,  0.000000e+00,  0.000000e+00/
      data  name(26),  ek(26),el(26)/
     $   'fe',  0.711200e+01,  0.842000e+00/
      data ak(0,26),ak(1,26),ak(2,26),ak(3,26)/
     $   0.143456e+02, -0.123491e+01, -0.423491e+00,  0.321661e-01/
      data al(0,26),al(1,26),al(2,26),al(3,26)/
     $   0.136696e+02, -0.239195e+01, -0.137680e+00,  0.000000e+00/
      data den(26),cf(26)/7.860,9.274e+01/
      data em(26),atwt(26)/0.094,55.850/
      data coh(0,26),coh(1,26),coh(2,26),coh(3,26)/
     $   0.593292e+01,  0.225038e+00, -0.361748e+00,  0.193024e-01/
      data cih(0,26),cih(1,26),cih(2,26),cih(3,26)/
     $  -0.342379e+00,  0.157245e+01, -0.253198e+00,  0.985822e-02/
      data ka(26),kb(26),la(26),lb(26)/
     $   0.640300e+01,  0.705700e+01,  0.000000e+00,  0.000000e+00/
      data name(27),  ek(27),el(27)/
     $  'co',  0.770900e+01,  0.929000e+00/
      data ak(0,27),ak(1,27),ak(2,27),ak(3,27)/
     $   0.147047e+02, -0.138933e+01, -0.386631e+00,  0.303286e-01/
      data al(0,27),al(1,27),al(2,27),al(3,27)/
     $   0.138699e+02, -0.250669e+01, -0.869945e-01,  0.000000e+00/
      data am(0,27),am(1,27),am(2,27),am(3,27)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data an(0,27),an(1,27),an(2,27),an(3,27)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den(27), cf(27)/ 0.890000e+01, 0.978500e+02/
      data em(27),atwt(27)/ 0.101000e+00,58.933/
      data coh(0,27),coh(1,27),coh(2,27),coh(3,27)/
     $   0.601478e+01,  0.237959e+00, -0.364056e+00,  0.194754e-01/
      data cih(0,27),cih(1,27),cih(2,27),cih(3,27)/
     $  -0.428804e+00,  0.164129e+01, -0.266013e+00,  0.106512e-01/
      data ka(27),kb(27),la(27),lb(27)/
     $   0.693000e+01,  0.764900e+01,  0.000000e+00,  0.000000e+00/
      data name(28),  ek(28),el(28)/
     $  'ni',  0.833300e+01,  0.101200e+01/
      data ak(0,28),ak(1,28),ak(2,28),ak(3,28)/
     $   0.142388e+02, -0.967736e+00, -0.478070e+00,  0.366138e-01/
      data al(0,28),al(1,28),al(2,28),al(3,28)/
     $   0.139848e+02, -0.248080e+01, -0.888115e-01,  0.000000e+00/
      data am(0,28),am(1,28),am(2,28),am(3,28)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data an(0,28),an(1,28),an(2,28),an(3,28)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den(28), cf(28)/ 0.890000e+01, 0.974500e+02/
      data em(28),atwt(28)/ 0.113000e+00,58.690/
      data coh(0,28),coh(1,28),coh(2,28),coh(3,28)/
     $   0.609024e+01,  0.222770e+00, -0.366568e+00,  0.196586e-01/
      data cih(0,28),cih(1,28),cih(2,28),cih(3,28)/
     $  -0.504360e+00,  0.170040e+01, -0.276443e+00,  0.112628e-01/
      data ka(28),kb(28),la(28),lb(28)/
     $   0.747700e+01,  0.826400e+01,  0.000000e+00,  0.000000e+00/
      data name(29),  ek(29),el(29)/
     $  'cu',  0.897900e+01,  0.110000e+01/
      data ak(0,29),ak(1,29),ak(2,29),ak(3,29)/
     $   0.145808e+02, -0.118375e+01, -0.413850e+00,  0.312088e-01/
      data al(0,29),al(1,29),al(2,29),al(3,29)/
     $   0.142439e+02, -0.258677e+01, -0.667398e-01,  0.000000e+00/
      data am(0,29),am(1,29),am(2,29),am(3,29)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data an(0,29),an(1,29),an(2,29),an(3,29)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den(29), cf(29)/ 0.894000e+01, 0.105500e+03/
      data em(29),atwt(29)/ 0.120000e+00,63.540/
      data coh(0,29),coh(1,29),coh(2,29),coh(3,29)/
     $   0.617739e+01,  0.273123e+00, -0.372360e+00,  0.201638e-01/
      data cih(0,29),cih(1,29),cih(2,29),cih(3,29)/
     $  -0.570210e+00,  0.175042e+01, -0.2845550e+00,  0.116930e-01/
      data ka(29),kb(29),la(29),lb(29)/
     $   0.804700e+01,  0.890400e+01,  0.000000e+00,  0.000000e+00/
      data name(30),  ek(30),el(30)/
     $  'zn',  0.965900e+01,  0.119600e+01/
      data ak(0,30),ak(1,30),ak(2,30),ak(3,30)/
     $   0.144118e+02, -0.933083e+00, -0.477357e+00,  0.362829e-01/
      data al(0,30),al(1,30),al(2,30),al(3,30)/
     $   0.143221e+02, -0.262384e+01, -0.264926e-01,  0.000000e+00/
      data am(0,30),am(1,30),am(2,30),am(3,30)/
     $   0.120597e+02, -0.110258e+01,  0.000000e+00,  0.000000e+00/
      data an(0,30),an(1,30),an(2,30),an(3,30)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den(30), cf(30)/ 0.714000e+01, 0.108600e+03/
      data em(30),atwt(30)/ 0.139000e+00,65.380/
      data coh(0,30),coh(1,30),coh(2,30),coh(3,30)/
     $   0.623402e+01,  0.284312e+00, -0.372143e+00,  0.200525e-01/
      data cih(0,30),cih(1,30),cih(2,30),cih(3,30)/
     $  -0.420535e+00,  0.163400e+01, -0.253646e+00,  0.927233e-02/
      data ka(30),kb(30),la(30),lb(30)/
     $   0.863800e+01,  0.957100e+01,  0.100900e+01,  0.103200e+01/
      data name(31),  ek(31),el(31)/
     $  'ga',  0.103670e+02,  0.130200e+01/
      data ak(0,31),ak(1,31),ak(2,31),ak(3,31)/
     $   0.136182e+02, -0.318459e+00, -0.611348e+00,  0.458138e-01/
      data al(0,31),al(1,31),al(2,31),al(3,31)/
     $   0.144792e+02, -0.254469e+01, -0.757204e-01,  0.000000e+00/
      data am(0,31),am(1,31),am(2,31),am(3,31)/
     $   0.122646e+02, -0.268965e+01,  0.000000e+00,  0.000000e+00/
c  this is pretty interesting so far, isn't it!
      data an(0,31),an(1,31),an(2,31),an(3,31)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data coh(0,31),coh(1,31),coh(2,31),coh(3,31)/
     $   0.628298e+01,  0.291334e+00, -0.369391e+00,  0.197029e-01/
      data cih(0,31),cih(1,31),cih(2,31),cih(3,31)/
     $  -0.358218e+00,  0.160050e+01, -0.244908e+00,  0.861898e-02/
      data den(31), cf(31)/ 0.590300e+01, 0.115800e+03/
      data em(31),atwt(31)/ 0.158000e+00,69.720/
      data ka(31),kb(31),la(31),lb(31)/
     $   0.925100e+01,  0.102630e+02,  0.109600e+01,  0.112200e+01/
      data name(32),  ek(32),el(32)/
     $  'ge',  0.111040e+02,  0.141400e+01/
      data ak(0,32),ak(1,32),ak(2,32),ak(3,32)/
     $   0.139288e+02, -0.479613e+00, -0.572897e+00,  0.431277e-01/
      data al(0,32),al(1,32),al(2,32),al(3,32)/
     $   0.146813e+02, -0.269285e+01, -0.208355e-01,  0.000000e+00/
      data am(0,32),am(1,32),am(2,32),am(3,32)/
     $   0.124133e+02, -0.253085e+01,  0.000000e+00,  0.000000e+00/
      data an(0,32),an(1,32),an(2,32),an(3,32)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data coh(0,32),coh(1,32),coh(2,32),coh(3,32)/
     $   0.633896e+01,  0.291512e+00, -0.365643e+00,  0.192896e-01/
      data cih(0,32),cih(1,32),cih(2,32),cih(3,32)/
     $  -0.334383e+00,  0.160327e+01, -0.245555e+00,  0.871239e-02/
      data den(32), cf(32)/ 0.532300e+01, 0.120500e+03/
      data em(32),atwt(32)/ 0.181000e+00,72.590/
      data ka(32),kb(32),la(32),lb(32)/
     $   0.988500e+01,  0.109810e+02,  0.118600e+01,  0.121600e+01/
      data name(33),  ek(33),el(33)/
     $  'as',  0.118680e+02,  0.153000e+01/
      data ak(0,33),ak(1,33),ak(2,33),ak(3,33)/
     $   0.134722e+02, -0.773513e-01, -0.660456e+00,  0.492177e-01/
      data al(0,33),al(1,33),al(2,33),al(3,33)/
     $   0.146431e+02, -0.248397e+01, -0.796180e-01,  0.000000e+00/
      data am(0,33),am(1,33),am(2,33),am(3,33)/
     $   0.125392e+02, -0.241380e+01,  0.000000e+00,  0.000000e+00/
      data an(0,33),an(1,33),an(2,33),an(3,33)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data den(33), cf(33)/ 0.573000e+01, 0.124400e+03/
      data em(33),atwt(33)/ 0.206000e+00,74.92/
      data coh(0,33),coh(1,33),coh(2,33),coh(3,33)/
     $   0.639750e+01,  0.288866e+00, -0.361747e+00,  0.188788e-01/
      data cih(0,33),cih(1,33),cih(2,33),cih(3,33)/
     $  -0.339189e+00,  0.162535e+01, -0.250783e+00,  0.909103e-02/
      data ka(33),kb(33),la(33),lb(33)/
     $   0.105430e+02,  0.117250e+02,  0.128200e+01,  0.131700e+01/
      data name(34),  ek(34),el(34)/
     $  'se',  0.126580e+02,  0.165300e+01/
      data ak(0,34),ak(1,34),ak(2,34),ak(3,34)/
     $   0.130756e+02,  0.183235e+00, -0.694264e+00,  0.502280e-01/
      data al(0,34),al(1,34),al(2,34),al(3,34)/
     $   0.147048e+02, -0.238853e+01, -0.105877e+00,  0.000000e+00/
      data am(0,34),am(1,34),am(2,34),am(3,34)/
     $   0.126773e+02, -0.239750e+01,  0.000000e+00,  0.000000e+00/
      data an(0,34),an(1,34),an(2,34),an(3,34)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data atwt(34),den(34),cf(34)/ 0.789600e+02, 0.479000e+01,
     $  0.131100e+03/
      data em(34)/ 0.232000e+00/
      data coh(0,34),coh(1,34),coh(2,34),coh(3,34)/
     $   0.645637e+01,  0.286737e+00, -0.358794e+00,  0.185618e-01/
      data cih(0,34),cih(1,34),cih(2,34),cih(3,34)/
     $  -0.432927e+00,  0.172833e+01, -0.277138e+00,  0.111735e-01/
      data ka(34),kb(34),la(34),lb(34)/
     $   0.112210e+02,  0.124950e+02,  0.141900e+01,  0.137900e+01/
      data name(35),  ek(35),el(35)/
     $  'br',  0.134740e+02,  0.178200e+01/
      data ak(0,35),ak(1,35),ak(2,35),ak(3,35)/
     $   0.132273e+02,  0.137130e+00, -0.683203e+00,  0.495424e-01/
      data al(0,35),al(1,35),al(2,35),al(3,35)/
     $   0.148136e+02, -0.242347e+01, -0.914590e-01,  0.000000e+00/
      data am(0,35),am(1,35),am(2,35),am(3,35)/
     $   0.127612e+02, -0.237730e+01,  0.000000e+00,  0.000000e+00/
      data an(0,35),an(1,35),an(2,35),an(3,35)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data atwt(35),den(35),cf(35)/ 0.799200e+02, 0.312000e+01,
     $  0.132700e+03/
      data em(35)/ 0.257000e+00/
      data coh(0,35),coh(1,35),coh(2,35),coh(3,35)/
     $   0.651444e+01,  0.286324e+00, -0.357027e+00,  0.183557e-01/
      data cih(0,35),cih(1,35),cih(2,35),cih(3,35)/
     $  -0.448001e+00,  0.176082e+01, -0.285099e+00,  0.117865e-01/
      data ka(35),kb(35),la(35),lb(35)/
     $   0.119230e+02,  0.132900e+02,  0.148000e+01,  0.152600e+01/
      data name(36),  ek(36),el(36),em(36)/
     $  'kr',  0.143220e+02,  0.192000e+01,  0.288000e+00/
      data ak(0,36),ak(1,36),ak(2,36),ak(3,36)/
     $   0.135927e+02, -0.305214e-01, -0.651340e+00,  0.477616e-01/
      data al(0,36),al(1,36),al(2,36),al(3,36)/
     $   0.149190e+02, -0.242418e+01, -0.876447e-01,  0.000000e+00/
      data am(0,36),am(1,36),am(2,36),am(3,36)/
     $   0.128898e+02, -0.226021e+01,  0.000000e+00,  0.000000e+00/
      data an(0,36),an(1,36),an(2,36),an(3,36)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data atwt(36),den(36),cf(36)/
     $  0.838000e+02, 0.374000e-02, 0.139100e+03/
      data coh(0,36),coh(1,36),coh(2,36),coh(3,36)/
     $   0.657113e+01,  0.287711e+00, -0.356311e+00,  0.182470e-01/
      data cih(0,36),cih(1,36),cih(2,36),cih(3,36)/
     $  -0.391810e+00,  0.173010e+01, -0.276824e+00,  0.111280e-01/
      data ka(36),kb(36),la(36),lb(36)/
     $   0.126480e+02,  0.141120e+02,  0.158700e+01,  0.163800e+01/
      data name(37),  ek(37),el(37),em(37)/
     $  'rb',  0.152000e+02,  0.206500e+01,  0.322000e+00/
      data ak(0,37),ak(1,37),ak(2,37),ak(3,37)/
     $   0.130204e+02,  0.382736e+00, -0.732427e+00,  0.529874e-01/
      data al(0,37),al(1,37),al(2,37),al(3,37)/
     $   0.149985e+02, -0.239108e+01, -0.959473e-01,  0.000000e+00/
      data am(0,37),am(1,37),am(2,37),am(3,37)/
     $   0.130286e+02, -0.238693e+01,  0.000000e+00,  0.000000e+00/
      data an(0,37),an(1,37),an(2,37),an(3,37)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data atwt(37),den(37),cf(37)/
     $  0.854800e+02, 0.153200e+01, 0.141900e+03/
      data coh(0,37),coh(1,37),coh(2,37),coh(3,37)/
     $   0.659750e+01,  0.302389e+00, -0.356755e+00,  0.181706e-01/
      data cih(0,37),cih(1,37),cih(2,37),cih(3,37)/
     $  -0.128039e+00,  0.153044e+01, -0.227403e+00,  0.739033e-02/
      data ka(37),kb(37),la(37),lb(37)/
     $   0.133940e+02,  0.149600e+02,  0.169400e+01,  0.175200e+01/
      data name(38),  ek(38),el(38),em(38)/
     $  'sr',  0.161050e+02,  0.221600e+01,  0.358000e+00/
      data ak(0,38),ak(1,38),ak(2,38),ak(3,38)/
     $   0.135888e+02,  0.220194e-02, -0.638940e+00,  0.460070e-01/
      data al(0,38),al(1,38),al(2,38),al(3,38)/
     $   0.150114e+02, -0.228169e+01, -0.126485e+00,  0.000000e+00/
      data am(0,38),am(1,38),am(2,38),am(3,38)/
     $   0.131565e+02, -0.236655e+01,  0.000000e+00,  0.000000e+00/
      data an(0,38),an(1,38),an(2,38),an(3,38)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data atwt(38),den(38),cf(38)/
     $  0.876200e+02, 0.254000e+01, 0.145500e+03/
      data coh(0,38),coh(1,38),coh(2,38),coh(3,38)/
     $   0.662203e+01,  0.324559e+00, -0.361651e+00,  0.184800e-01/
      data cih(0,38),cih(1,38),cih(2,38),cih(3,38)/
     $   0.799161e-01,  0.138397e+01, -0.192225e+00,  0.478611e-02/
      data ka(38),kb(38),la(38),lb(38)/
     $   0.141640e+02,  0.158340e+02,  0.180600e+01,  0.187200e+01/
      data name(39),  ek(39),el(39),em(39)/
     $  'y ',  0.170800e+02,  0.237300e+01,  0.395000e+00/
      data ak(0,39),ak(1,39),ak(2,39),ak(3,39)/
     $   0.134674e+02,  0.191023e+00, -0.686616e+00,  0.497356e-01/
      data al(0,39),al(1,39),al(2,39),al(3,39)/
     $   0.151822e+02, -0.238946e+01, -0.881174e-01,  0.000000e+00/
      data am(0,39),am(1,39),am(2,39),am(3,39)/
     $   0.132775e+02, -0.243174e+01,  0.000000e+00,  0.000000e+00/
      data an(0,39),an(1,39),an(2,39),an(3,39)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data atwt(39),den(39),cf(39)/
     $  0.889050e+02, 0.440500e+01, 0.147600e+03/
      data coh(0,39),coh(1,39),coh(2,39),coh(3,39)/
     $   0.667096e+01,  0.325070e+00, -0.360613e+00,  0.183326e-01/
      data cih(0,39),cih(1,39),cih(2,39),cih(3,39)/
     $   0.629057e-01,  0.141577e+01, -0.199713e+00,  0.533312e-02/
      data ka(39),kb(39),la(39),lb(39)/
     $   0.149570e+02,  0.167360e+02,  0.192200e+01,  0.199600e+01/
      data name(40),  ek(40),el(40),em(40)/
     $  'zr',  0.179980e+02,  0.253200e+01,  0.431000e+00/
      data ak(0,40),ak(1,40),ak(2,40),ak(3,40)/
     $   0.127538e+02,  0.697409e+00, -0.789307e+00,  0.564531e-01/
      data al(0,40),al(1,40),al(2,40),al(3,40)/
     $   0.152906e+02, -0.238703e+01, -0.912292e-01,  0.000000e+00/
      data am(0,40),am(1,40),am(2,40),am(3,40)/
     $   0.134508e+02, -0.250201e+01,  0.000000e+00,  0.000000e+00/
      data an(0,40),an(1,40),an(2,40),an(3,40)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data atwt(40),den(40),cf(40)/
     $  0.912200e+02, 0.653000e+01, 0.151500e+03/
      data coh(0,40),coh(1,40),coh(2,40),coh(3,40)/
     $   0.672275e+01,  0.323964e+00, -0.359463e+00,  0.181890e-01/
      data cih(0,40),cih(1,40),cih(2,40),cih(3,40)/
     $   0.366697e-01,  0.145207e+01, -0.208122e+00,  0.595139e-02/
      data ka(40),kb(40),la(40),lb(40)/
     $   0.157740e+02,  0.176660e+02,  0.204200e+01,  0.212400e+01/
      data name(41),  ek(41),el(41),em(41)/
     $  'nb',  0.189860e+02,  0.269800e+01,  0.468000e+00/
      data ak(0,41),ak(1,41),ak(2,41),ak(3,41)/
     $   0.133843e+02,  0.281028e+00, -0.686607e+00,  0.486607e-01/
      data al(0,41),al(1,41),al(2,41),al(3,41)/
     $   0.152088e+02, -0.220278e+01, -0.136759e+00,  0.000000e+00/
      data am(0,41),am(1,41),am(2,41),am(3,41)/
     $   0.135434e+02, -0.250135e+01,  0.000000e+00,  0.000000e+00/
      data an(0,41),an(1,41),an(2,41),an(3,41)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data atwt(41),den(41),cf(41)/
     $  0.929060e+02, 0.857000e+01, 0.154300e+03/
      data coh(0,41),coh(1,41),coh(2,41),coh(3,41)/
     $   0.679013e+01,  0.311282e+00, -0.355233e+00,  0.178231e-01/
      data cih(0,41),cih(1,41),cih(2,41),cih(3,41)/
     $   0.202289e-03,  0.149347e+01, -0.217419e+00,  0.662245e-02/
      data ka(41),kb(41),la(41),lb(41)/
     $   0.166140e+02,  0.186210e+02,  0.216600e+01,  0.225700e+01/
      data name(42),  ek(42),el(42),em(42)/
     $  'mo',  0.199990e+02,  0.286600e+01,  0.505000e+00/
      data ak(0,42),ak(1,42),ak(2,42),ak(3,42)/
     $   0.139853e+02, -0.117426e+00, -0.591094e+00,  0.417843e-01/
      data al(0,42),al(1,42),al(2,42),al(3,42)/
     $   0.153494e+02, -0.226640e+01, -0.116881e+00,  0.000000e+00/
      data am(0,42),am(1,42),am(2,42),am(3,42)/
     $   0.136568e+02, -0.248480e+01,  0.000000e+00,  0.000000e+00/
      data an(0,42),an(1,42),an(2,42),an(3,42)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data atwt(42),den(42),cf(42)/
     $  0.959500e+02, 0.102200e+02, 0.159300e+03/
      data coh(0,42),coh(1,42),coh(2,42),coh(3,42)/
     $   0.684600e+01,  0.302790e+00, -0.351131e+00,  0.174403e-01/
      data cih(0,42),cih(1,42),cih(2,42),cih(3,42)/
     $  -0.562860e-01,  0.155778e+01, -0.233341e+00,  0.785506e-02/
      data ka(42),kb(42),la(42),lb(42)/
     $   0.174780e+02,  0.196070e+02,  0.229300e+01,  0.239500e+01/
      data name(43),  ek(43),el(43),em(43)/
     $  'tc',  0.210450e+02,  0.304300e+01,  0.544000e+00/
      data ak(0,43),ak(1,43),ak(2,43),ak(3,43)/
     $   0.128214e+02,  0.751993e+00, -0.787006e+00,  0.558668e-01/
      data al(0,43),al(1,43),al(2,43),al(3,43)/
     $   0.155086e+02, -0.233733e+01, -0.987857e-01,  0.000000e+00/
      data am(0,43),am(1,43),am(2,43),am(3,43)/
     $   0.137498e+02, -0.244730e+01,  0.000000e+00,  0.000000e+00/
      data an(0,43),an(1,43),an(2,43),an(3,43)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data atwt(43),den(43),cf(43)/
     $  0.990000e+02, 0.115000e+02, 0.164400e+03/
      data coh(0,43),coh(1,43),coh(2,43),coh(3,43)/
     $   0.687599e+01,  0.326165e+00, -0.358969e+00,  0.180482e-01/
      data cih(0,43),cih(1,43),cih(2,43),cih(3,43)/
     $   0.757616e-01,  0.144950e+01, -0.204890e+00,  0.564745e-02/
      data ka(43),kb(43),la(43),lb(43)/
     $   0.184100e+02,  0.205850e+02,  0.242400e+01,  0.253800e+01/
      data name(44),  ek(44),el(44),em(44)/
     $  'ru',  0.221170e+02,  0.322400e+01,  0.585000e+00/
      data ak(0,44),ak(1,44),ak(2,44),ak(3,44)/
     $   0.126658e+02,  0.885020e+00, -0.811144e+00,  0.573759e-01/
      data al(0,44),al(1,44),al(2,44),al(3,44)/
     $   0.154734e+02, -0.223080e+01, -0.119454e+00,  0.000000e+00/
      data am(0,44),am(1,44),am(2,44),am(3,44)/
     $   0.138782e+02, -0.248066e+01,  0.000000e+00,  0.000000e+00/
      data an(0,44),an(1,44),an(2,44),an(3,44)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data atwt(44),den(44),cf(44)/
     $  0.101070e+03, 0.124100e+02, 0.167800e+03/
      data coh(0,44),coh(1,44),coh(2,44),coh(3,44)/
     $   0.693136e+01,  0.334794e+00, -0.363497e+00,  0.184429e-01/
      data cih(0,44),cih(1,44),cih(2,44),cih(3,44)/
     $  -0.424981e-01,  0.154639e+01, -0.226470e+00,  0.718375e-02/
      data ka(44),kb(44),la(44),lb(44)/
     $   0.192780e+02,  0.216550e+02,  0.255800e+01,  0.268300e+01/
      data name(45),  ek(45),el(45),em(45)/
     $  'rh',  0.232200e+02,  0.341200e+01,  0.627000e+00/
      data ak(0,45),ak(1,45),ak(2,45),ak(3,45)/
     $   0.121760e+02,  0.119682e+01, -0.866697e+00,  0.606931e-01/
      data al(0,45),al(1,45),al(2,45),al(3,45)/
     $   0.155757e+02, -0.224976e+01, -0.113377e+00,  0.000000e+00/
      data am(0,45),am(1,45),am(2,45),am(3,45)/
     $   0.140312e+02, -0.261303e+01,  0.000000e+00,  0.000000e+00/
      data an(0,45),an(1,45),an(2,45),an(3,45)/
     $   0.000000e+00,  0.000000e+00,  0.000000e+00,  0.000000e+00/
      data atwt(45),den(45),cf(45)/
     $  0.102910e+03, 0.124400e+02, 0.170900e+03/
      data coh(0,45),coh(1,45),coh(2,45),coh(3,45)/
     $   0.697547e+01,  0.346394e+00, -0.367794e+00,  0.187885e-01/
      data cih(0,45),cih(1,45),cih(2,45),cih(3,45)/
     $  -0.160399e+00,  0.164861e+01, -0.250238e+00,  0.893818e-02/
      data ka(45),kb(45),la(45),lb(45)/
     $   0.202140e+02,  0.227210e+02,  0.269600e+01,  0.283400e+01/
      data name(46),  ek(46),el(46),em(46)/
     $  'pd',  0.243500e+02,  0.360500e+01,  0.670000e+00/
      data ak(0,46),ak(1,46),ak(2,46),ak(3,46