c***********************************************************************
c
      program fragyld
c
      include 'fragyld.par'
c
c   fragyld.par:
c   large model (7 MB):
c      parameter(kdz=20,kdn=40,kst=40)
c      parameter(nex=500)
c   small model (600 kB):
c      parameter(kdz=15,kdn=15,kst=15)
c      parameter(nex=100)
c   the rest:
c      parameter(kem=4)
c      parameter(ns=30,ntt=300,nst=2*ntt*(ns+1))
c      parameter(npr=3,prec=0.001,epsy=10.*prec)
c
      character symaz*5,direct*30
      dimension ylda(8,kdn+kdz),yldg(3,kdn,kdz),yldpz(8)
      common/nucs/cn(4,nex,kdn,kdz),nmax(kdn,kdz),za(2),acx(3),de,ip,p0
      common/hilo/yldp(5,kdn,kdz),hyld(kdn,kdz),nzn(2),fcor,ncor
      data amu/939.0/
c
c      direct='[carlson.frag]'
c      lnd=14
c
c      direct='D:\Myfor\Frag\Data\'
c      lnd=19
c
      direct='E:\Brett\Myfor\Frag\Data\'
      lnd=25
c
      open(unit=5,file=direct(1:lnd)//'fragyld.in',status='old')
      open(unit=8,file=direct(1:lnd)//'fragyld.out',status='unknown')
c      open(unit=9,file=direct(1:lnd)//'fragyld.pzn',status='unknown')
c
 10   time=second()
      call init(at,inc,iprpr)
      call abrade
      call ablate(inc)
c
      arat=at/(za(2)+at)
      neu=za(2)-za(1)
c
      do 30 ja=1,nzn(1)+nzn(2)
        do 20 i=1,8
          ylda(i,ja)=0.
 20     continue
 30   continue
      iz=za(1)+1
c
      if(iprpr.ne.0) then
        call geoyld(za(1),za(2),at,yldg)
        decor2=fcor*(1.-fcor)*ncor**2
        decor=(1.-fcor)*ncor
        write(8,810)
        do 60 jz=1,nzn(1)
          iz=iz-1
          ia0=iz+neu+1
          do 35 k=1,8
            yldpz(k)=0.
 35       continue
          jzn=1
          if(jz.eq.1) jzn=2
          do 50 jn=jzn,nzn(2)
            ja=jn+jz-1
            do 40 k=1,5
              ylda(k,ja)=ylda(k,ja)+yldp(k,jn,jz)
 40         continue
            ylda(6,ja)=ylda(6,ja)+yldg(1,jn,jz)
            ylda(7,ja)=ylda(7,ja)+yldg(1,jn,jz)*yldg(2,jn,jz)
            ylda(8,ja)=ylda(8,ja)+yldg(1,jn,jz)*yldg(3,jn,jz)
            do 45 k=1,5
              yldpz(k)=yldpz(k)+yldp(k,jn,jz)
 45         continue
            yldpz(6)=yldpz(6)+yldg(1,jn,jz)
            yldpz(7)=yldpz(7)+yldg(1,jn,jz)*yldg(2,jn,jz)
            yldpz(8)=yldpz(8)+yldg(1,jn,jz)*yldg(3,jn,jz)
            ia=ia0-jn
            if(yldp(1,jn,jz).gt.epsy) then
              call chsym(iz,ia,symaz)
              yldp(2,jn,jz)=yldp(2,jn,jz)/yldp(1,jn,jz)
              yldp(3,jn,jz)=yldp(3,jn,jz)/yldp(1,jn,jz)-yldp(2,jn,jz)**2
              yldp(3,jn,jz)=sqrt(yldp(3,jn,jz))
              yldp(4,jn,jz)=yldp(4,jn,jz)/yldp(1,jn,jz)
              yldp(5,jn,jz)=yldp(5,jn,jz)/yldp(1,jn,jz)-yldp(4,jn,jz)**2
              yldp(4,jn,jz)=de*(yldp(4,jn,jz)+decor)
              yldp(5,jn,jz)=de*sqrt(yldp(5,jn,jz)+decor2)
              write(8,820) symaz,(yldp(k,jn,jz),k=1,5),
     1                           (yldg(k,jn,jz),k=1,3)
             endif
 50       continue
          if(yldpz(1).gt.epsy) then
            call chsym(iz,ia,symaz)
            symaz(1:3)='   '
              yldpz(2)=yldpz(2)/yldpz(1)
              yldpz(3)=yldpz(3)/yldpz(1)-yldpz(2)**2
              yldpz(3)=sqrt(yldpz(3))
              yldpz(4)=yldpz(4)/yldpz(1)
              yldpz(5)=yldpz(5)/yldpz(1)-yldpz(4)**2
              yldpz(4)=de*(yldpz(4)+decor)
              yldpz(5)=de*sqrt(yldpz(5)+decor2)
              if(yldpz(6).gt.1.0e-6) then
                yldpz(7)=yldpz(7)/yldpz(6)
                yldpz(8)=yldpz(8)/yldpz(6)
               else
                yldpz(7)=0.
                yldpz(8)=0.
               endif
            write(8,830) symaz,(yldpz(k),k=1,8)
           endif
          if(jz.eq.1) sgpz=yldp(1,1,1)-yldpz(1)
 60     continue
        symaz='chrge'
        write(8,830) symaz,sgpz
c
        ia=za(2)
        write(8,850)
        do 70 ja=2,nzn(1)+nzn(2)
          ia=ia-1
          if(ylda(1,ja).gt.epsy) then
            ylda(2,ja)=ylda(2,ja)/ylda(1,ja)
            ylda(3,ja)=ylda(3,ja)/ylda(1,ja)-ylda(2,ja)**2
            ylda(3,ja)=sqrt(ylda(3,ja))
            ylda(4,ja)=ylda(4,ja)/ylda(1,ja)
            ylda(5,ja)=ylda(5,ja)/ylda(1,ja)-ylda(4,ja)**2
            ylda(4,ja)=de*(ylda(4,ja)+decor)
            ylda(5,ja)=de*sqrt(ylda(5,ja)+decor2)
            if(ylda(6,ja).gt.1.0e-6) then
              ylda(7,ja)=ylda(7,ja)/ylda(6,ja)
              ylda(8,ja)=ylda(8,ja)/ylda(6,ja)
             else
              ylda(7,ja)=0.
              ylda(8,ja)=0.
             endif
          write(8,860) ia,(ylda(k,ja),k=1,8)
         endif
 70     continue
        symaz='mass '
        write(8,830) symaz,yldp(1,1,1)
       endif
c
      do 90 ja=1,nzn(1)+nzn(2)
        do 80 i=1,4
          ylda(i,ja)=0.
 80     continue
 90   continue
      iz=za(1)+1
c
      write(8,910)
      do 120 jz=1,nzn(1)
        iz=iz-1
        ia0=iz+neu+1
        yldz=0.
        pavz=0.
        pparz=0.
        pperz=0.
        jzn=1
        if(jz.eq.1) jzn=2
        do 110 jn=jzn,nzn(2)
          ja=jn+jz-1
          do 100 i=1,4
            ylda(i,ja)=ylda(i,ja)+cn(i,1,jn,jz)
 100      continue
          ia=ia0-jn
          if(cn(1,1,jn,jz).gt.epsy) then
            call chsym(iz,ia,symaz)
            pav=1.-arat*(1.-cn(2,1,jn,jz)/cn(1,1,jn,jz))
            ppar=2.*amu*cn(3,1,jn,jz)/cn(1,1,jn,jz)
            pper=sqrt(cn(4,1,jn,jz)/cn(1,1,jn,jz)+ppar)
            ppar=sqrt(ppar)
            write(8,920) symaz,cn(1,1,jn,jz),pav,ppar,pper
           endif
          yldz=yldz+cn(1,1,jn,jz)
          pavz=pavz+cn(2,1,jn,jz)
          pparz=pparz+cn(3,1,jn,jz)
          pperz=pperz+cn(4,1,jn,jz)
 110    continue
        if(yldz.gt.epsy) then
          call chsym(iz,ia,symaz)
          symaz(1:3)='   '
          pavz=1.-arat*(1.-pavz/yldz)
          pparz=2.*amu*pparz/yldz
          pperz=sqrt(pperz/yldz+pparz)
          pparz=sqrt(pparz)
          write(8,930) symaz,yldz,pavz,pparz,pperz
         endif
        if(jz.eq.1) sgz=yldp(1,1,1)-yldz
 120  continue
      symaz='chrge'
      write(8,930) symaz,sgz
c
      ia=za(2)
      write(8,950)
      do 130 ja=2,nzn(1)+nzn(2)
        ia=ia-1
        if(ylda(1,ja).gt.epsy) then
          pav=1.-arat*(1.-ylda(2,ja)/ylda(1,ja))
          ppar=2.*amu*ylda(3,ja)/ylda(1,ja)
          pper=sqrt(ylda(4,ja)/ylda(1,ja)+ppar)
          ppar=sqrt(ppar)
          write(8,960) ia,ylda(1,ja),pav,ppar,pper
         endif
 130   continue
      symaz='mass '
      write(8,930) symaz,yldp(1,1,1)
      time=second()-time
      write(8,800) time
      go to 10
 780  format(1h1)
 800  format(///,' Elapsed time = ',f8.1,' seconds',///)
 810  format(///,21x,'Combinatorial/Geometrical Primary Yields',/,
     1          ' Isotope',3x,'Prim.Yld (mb)',4x,
     2          'Ave.b (Fm)',3x,'Std.d.b (Fm)',2x,
     3          'Ave.E (MeV)',2x,'Std.d.E (MeV)',/)
 820  format(2x,a5,2x,5(f9.3,5x),/,15x,2(f9.3,5x),14x,f9.3)
 830  format(/,2x,a5,2x,5(f9.3,5x),/,15x,2(f9.3,5x),14x,f9.3,//)
 850  format(///,21x,'Combinatorial/Geometrical Primary Yields',/,
     1          '  Mass  ',3x,'Prim.Yld (mb)',4x,
     2          'Ave.b (Fm)',3x,'Std.d.b (Fm)',2x,
     3          'Ave.E (MeV)',2x,'Std.d.E (MeV)',/)
 860  format(2x,i3,4x,5(f9.3,5x),/,15x,2(f9.3,5x),14x,f9.3)
 910  format(///,29x,'Secondary Yields',/,' Isotope',4x,
     1          'Sec.Yld (mb)',4x,'V_f/V_i',4x,
     2          'Par.Wd. (MeV/c)',2x,'Perp.Wd. (MeV/c)',/)
 920  format(2x,a5,5x,f9.3,5x,f9.5,6x,f9.3,8x,f9.3)
 930  format(/,2x,a5,5x,f9.3,5x,f9.5,6x,f9.3,8x,f9.3,//)
 950  format(///,29x,'Secondary Yields',/,'  Mass  ',4x,
     1          'Sec.Yld (mb)',4x,'V_f/V_i',4x,
     2          'Par.Wd. (MeV/c)',2x,'Perp.Wd. (MeV/c)',/)
 960  format(2x,i3,2x,5x,f9.3,5x,f9.5,6x,f9.3,8x,f9.3)
      end
c
c-----------------------------------------------------------------------====
c
      subroutine init(at,inc,iprpr)
c
      include 'fragyld.par'
c
      character symp*5,symt*5      
      common/nucs/cn(4,nex,kdn,kdz),nmax(kdn,kdz),za(2),acx(3),de,ip,p0
      common/hilo/yldp(5,kdn,kdz),hyld(kdn,kdz),nzn(2),fcor,ncor
      common/ind/eps(ns,2),jg(ns,2),nx(2),izn(2),nzf(2)
      common/ees/el(kst,2),ne(kst,2),nxf
      common/tind/nn(ns,2),ll(ns,2),jj(ns,2),b0,iwt
      common/dense/dc(4),w(4),r(4),a(4),itp(4),ird,iwrt
      common/thick/ipzn,jtmax,jpmax,jbmax,dr,sigpp,sigpn,tt(ntt,2)
     1                                               ,tp(ntt,ns,2)
      data amu/939.0/
      write(8,'('' Echo Print of Input ****************************'')')
c-----------------------------------------------------------------------====
c   basic parameters of calculation -- cannot be set to 0 (all others can)
c  zp,ap,zt,at - charge and mass of projectile and target
c  ep          - projectile energy / nucleon
c
      read(5,*,end=80) zp,ap,zt,at,ep
      write(*,*) zp,ap,zt,at,ep
      write(8,*) zp,ap,zt,at,ep
      if(zp+ap.lt.0.5) go to 80
      za(1)=zp
      izn(1)=zp
      za(2)=ap
      anp=ap-zp
      izn(2)=anp
      f0=1.44*(ap+at)/(ep*at*(1.07*(ap**(1./3.)+at**(1./3.))+2.7))
c-----------------------------------------------------------------------====
c 
c          flags controlling thickness calculation
c  iwt  - <0:geom. thickness; =0:s.p. thick, no ang.mom.; >0:s.p. w/ ang.mom.
c  ipzn - =0: thickness as in PRC51(1995)252 w/ iwt>=0
c                       as in PRC46(1992)R30 w/ iwt<0
c         >0: ??? 
c  ird  - !=0: read orbital psi's
c  iwrt - !=0: write detailed info on thickness
c  iprpr- !=0: prints primary distributions
c
      read(5,*) iwt,ipzn,ird,iwrt,iprpr
      if(ipzn.eq.0) ipzn=2
      if(ipzn.gt.2) ipzn=2
      write(8,*) iwt,ipzn,ird,iwrt
c-----------------------------------------------------------------------====
c        parameters describing Z, N dist. of target (always) and
c        proj. (geom. abrasion)
c  itp - =1: Wood-Saxon dist.; =2 Gaussian dist. for Zp, Np, Zt, Nt resp. 
c  r   - >0: red. radius (WS) or rms width (G); <0: -abs.radius/rms width
c  a   - diff. (WS) or coeff. of r^2 term (G)
c  w   - coeff. of r^2 (WS)
c
      read(5,*) itp(3),itp(4),itp(1),itp(2)
      read(5,*) r(3),r(4),r(1),r(2)
      read(5,*) a(3),a(4),a(1),a(2)
      read(5,*) w(3),w(4),w(1),w(2)
      zi=zt
      ai=at
      do 10 i=1,4
        if(i.gt.2) then
          zi=zp
          ai=ap
         endif
        if(itp(i).lt.1) then
          itp(i)=1
          if(zi.lt.8.5) itp(i)=2
         endif
        if(r(i).eq.0.) then
          if(itp(i).gt.1) then
            r(i)=0.72
            a(i)=amax1(0.,0.25*(zi-2.))
           else
            r(i)=1.11-1.5/ai
           endif
         endif
        if(a(i).lt.0.05.and.itp(i).eq.1) a(i)=0.57
 10   continue
      write(8,*) itp(3),itp(4),itp(1),itp(2)
      write(8,*) r(3),r(4),r(1),r(2)
      write(8,*) a(3),a(4),a(1),a(2)
      write(8,*) w(3),w(4),w(1),w(2)
c-----------------------------------------------------------------------====
c     
c         indices controlling abrasion calculation
c   nzn - max. number of abraded prot./neut., resp.ly taken into account.
c         The ablation calculation is also limited to these residuals.
c   nzf - max. number of prot./neut. for combinatorial calc. of E dist.
c         for z,n >nzf(1,2), E dist. is calculated using moments. These
c         can usually be set to about 6 with almost no loss of accuracy.
c   est - energy below which E dist. calculated combinatorially even for
c         z,n > nzf(1,2)
c
      read(5,*) nzn,nzf,est
      if(nzn(1).lt.1) nzn(1)=izn(1)
      if(nzn(2).lt.1) nzn(2)=izn(2)
      nzn(1)=max0(min0(nzn(1),izn(1)-2,kst,kdz),2)
      nzn(2)=max0(min0(nzn(2),izn(2)-2,kst,kdn),2)
      if(nzf(1).lt.1) nzf(1)=izn(1)
      if(nzf(2).lt.1) nzf(2)=izn(2)
      nzf(1)=min0(nzf(1),nzn(1))
      nzf(2)=min0(nzf(2),nzn(2))
      if(est.lt.0.1) est=20.
      write(8,*) nzn,nzf,est
c-----------------------------------------------------------------------====
c   nx    - number of p,n orbitals to be input. If either nx(1) or nx(2)
c           is !=0, then both p and n orbitals must be input. If both nx(1)
c           and nx(2) =0, then the orbitals are calculated internally.
c   spfac - factor multiplying s.p. energies, used to increase excitation E 
c
      read(5,*) nx,spfac
      if(spfac.lt.1.0e-8) spfac=1.
      write(8,*) nx,spfac
c-----------------------------------------------------------------------====
c     Input of orbitals
c   eps -- s.p. energy of orbital
c   nn  -- principal q. number of orbital, (nn=1 for 1s1/2)
c   ll  -- 2*orbital ang.mom. of orbital
c   jj  -- 2* total ang.mom. of orbital, (jj=1 for 1s1/2)
c
      if(nx(1)+nx(2).gt.0) then
        do 30 k=1,2
          if(nx(k).gt.ns) go to 60
          read(5,*) (eps(n,k),nn(n,k),ll(n,k),jj(n,k),n=1,nx(k))
          write(8,*) (eps(n,k),nn(n,k),ll(n,k),jj(n,k),n=1,nx(k))
          ioc=izn(k)
          do 20 n=nx(k),1,-1
            eps(n,k)=spfac*eps(n,k)
            jg(n,k)=jj(n,k)+1
            if(ioc.lt.jg(n,k)) jg(n,k)=ioc
            ioc=ioc-jg(n,k)
 20       continue
          if(ioc.gt.0) go to 70
 30     continue
       erpa=0.
       else
        call spni(spfac,erpa)
       endif
c-----------------------------------------------------------------------====
c   Parameters of momentum distribution of proj. fragments
c   sig0 - ave.perp.mom. transferred in one-nucleon removal(Goldberger-Watson) 
c   fcor - fraction of E dist. NOT shifted by ecor (??). Best kept =1.
c   ecor - energy shift of (1-fcor) of abraded e dist. (??). Best kept =0.
c
      read(5,*) sig0,fcor,ecor
      if(sig0.lt.0.001) sig0=80.
      if(fcor.eq.0.0) fcor=0.6
      if(fcor.lt.1.0.and.ecor.eq.0.0) ecor=erpa
      write(8,*) sig0,fcor,ecor
      sig0=0.5*sig0*sig0/amu
c-----------------------------------------------------------------------====
c   sigpp,sigpn - proton-proton and proton-neutron cross sections
c                 at ep MeV/nucleon
c
      read(5,*) sigpp,sigpn
      if(sigpp.le.0.) call signn(ep,sigpp,sigpn)
      write(8,*) sigpp,sigpn
      sigpp=sigpp/10.
      sigpn=sigpn/10.
      if(ipzn.eq.2) then
        sigpp=0.5*sigpp
        sigpn=0.5*sigpn
       endif
c-----------------------------------------------------------------------====
c     Ablation parameters
c   inc - no. of particle types emitted in stat. decay; =1:gamma;
c         =2: gamma, n; =3:gamma, n, p; =4:gamma, n, p, alpha, etc.
c   ip  - >1: prints details of ablation calculation; >2:prints more
c         details of abrasion calc. than you probably want to see.
c
      read(5,*) inc,ip
      if(inc.lt.1) inc=4
      if(inc.gt.kem) inc=kem
      write(8,*) inc,ip
c-----------------------------------------------------------------------====
c     Level density parameters
c   de     - energy bin size for abrasion/ablation calc.
c   acx(1) - level density parameter in a=A/acx(1). acx(1)=8 reasonable.
c   acx(2) - pairing gap parameter in delta = acx(2)/sqrt(A).
c   acx(3) - =0: Fermi gas level density; =1: Gilbert-Cameron level density
c
      read(5,*) de,acx
      xz=amin1(nzn(1)/zp,0.5)
      xz=1.-xz**1.5-(1.-xz)**1.5
      xn=amin1(nzn(2)/anp,0.5)
      xn=1.-xn**1.5-(1.-xn)**1.5
      xn=24.*(zp*xz+anp*xn)
      de=amax1(de,xn/(nex-1))
      if(acx(1).lt.0.) acx(2)=12.
      if(acx(1).eq.0.) then
        acx(1)=8.0
        if(ap.lt.40.) acx(1)=7.1
        if(acx(2).eq.0.) acx(2)=12.
       endif
      if(acx(3).gt.0.5) then
        acx(3)=1.
       else
        acx(3)=0.
       endif
      write(8,*) de,acx
c-----------------------------------------------------------------------====
      write(8,'('' End of Input ***********************************'')')
      nxf=est/de+1.5
      ncor=ecor/de+0.5
      if(ncor.lt.1) fcor=1.
      call chsym(int(zp),int(ap),symp)
      call chsym(int(zt),int(at),symt)
      write(8,'(//,'' Projectile fragmentation of '',a5,
     1            '' incident on '',a5,'' at'',f7.1,
     2            '' MeV/nucleon'')') symp,symt,ep
      call wrsp(de)
      call precomb(zp,ap,zt,at,f0,sig0)
      call thprep(ap,zt,at)
      call emdiss(zp,ap,zt,at,ep)
      return
 60   write(8,800)
      stop
 70   write(8,810)
      stop
 80   write(8,820)
      stop
 800  format(' Single particle array dimension was exceeded in READIT.')
 810  format(' Too few single particle levels were input.')
 820  format(' All done.')
      end
c
c-----------------------------------------------------------------------====
c
      subroutine spni(spfac,erpa)
c
      include 'fragyld.par'
      parameter(nss=40)
      dimension esp(nss),jsp(nss),lsp(nss),nsp(nss)
      common/ind/eps(ns,2),jg(ns,2),nx(2),izn(2),nzf(2)
      common/tind/nn(ns,2),ll(ns,2),jj(ns,2),b0,iwt
      data esp/1.500,2.450,2.600,3.400,3.500,3.650,4.140,4.415,4.490,
     1   4.565,4.850,5.265,5.300,5.500,5.515,5.575,6.080,6.125,6.259,
     2   6.405,6.430,6.555,6.852,6.909,6.935,7.266,7.302,7.500,7.516,
     3   7.538,7.599,7.685,8.090,8.097,8.149,8.289,8.388,8.407,8.440,
     4   8.557/
      data jsp/ 1, 3, 1, 5, 1, 3, 7, 3, 5, 1, 9, 5, 7, 1, 3,11, 7, 9,
     1   13, 3, 5, 1, 9,11,15, 5, 7, 1, 3,17,11,13, 7,19, 9,13,15, 3,
     2    5, 1/
      data lsp/ 0, 2, 2, 4, 0, 4, 6, 2, 6, 2, 8, 4, 8, 0, 4,10, 6,10,
     1   12, 2, 6, 2, 8,12,14, 4, 8, 0, 4,16,10,14, 6,18,10,12,16, 2,
     2    6, 2/
      data nsp/ 1, 1, 1, 1, 2, 1, 1, 2, 1, 2, 1, 2, 1, 3, 2, 1, 2, 1,
     1    1, 3, 2, 3, 2, 1, 1, 3, 2, 4, 3, 1, 2, 1, 3, 1, 2, 2, 1, 4,
     2    3, 4/
      hw=izn(1)+izn(2)
      b0=1.01*hw**(1./6.)
c      hw=60./hw**(1./3.)
      hw=40.*spfac/hw**(1./3.)
      erpa=10.
      do 60 k=1,2
        jsum=0
        ioc=izn(k)
        do 30 n=1,nss
          nf=n
          jsum=jsum+jsp(n)+1
          if(jsum.ge.ioc) go to 40
 30     continue
 40     if(nf.gt.ns) go to 80
        nx(k)=nf
        ef=esp(nf)
        nf=nf+1
        erpa=amin1(max0(0,ioc-jsum+2)*(esp(nf)-ef),erpa)
        do 50 n=nx(k),1,-1
          eps(n,k)=hw*(ef-esp(nf-n))
          nn(n,k)=nsp(nf-n)
          ll(n,k)=lsp(nf-n)
          jj(n,k)=jsp(nf-n)
          jg(n,k)=jsp(nf-n)+1
          if(ioc.lt.jg(n,k)) jg(n,k)=ioc
          ioc=ioc-jg(n,k)
 50     continue
 60   continue
      erpa=hw*erpa
      return
 80   write(8,800)
      stop
 800  format(' Single particle array dimension was exceeded in SPNI.')
      end
c
c-----------------------------------------------------------------------====
c
      subroutine wrsp(de)
c
      include 'fragyld.par'
      dimension aor(10)
      common/ind/eps(ns,2),jg(ns,2),nx(2),izn(2),nzf(2)
      common/tind/nn(ns,2),ll(ns,2),jj(ns,2),b0,iwt
      data aor/'s','p','d','f','g','h','i','j','k','l'/
      n2=min0(nx(1),nx(2))
      write(8,810)
      write(8,820)
     1((eps(n,k),nn(n,k),aor(1+ll(n,k)/2),jj(n,k),jg(n,k),k=1,2),n=1,n2)
      n2=n2+1
      if(nx(1).gt.nx(2)) write(8,830)
     1    (eps(n,1),nn(n,1),aor(1+ll(n,1)/2),jj(n,1),jg(n,1),n=n2,nx(1))
      if(nx(2).gt.nx(1)) write(8,840)
     1    (eps(n,2),nn(n,2),aor(1+ll(n,2)/2),jj(n,2),jg(n,2),n=n2,nx(2))
      do 20 k=1,2
        do 10 n=1,nx(k)
          eps(n,k)=eps(n,k)/de
 10     continue
 20   continue
      return
 810  format(///,30x,'Projectile Levels',/,21x,'Proton',24x,'Neutron',/,
     1  9x,2(5x,'Energy',4x,'Orbit',2x,'Occ',5x))
 820  format((8x,2(5x,f7.3,3x,i1,a1,i2,'/2',1x,i3,5x)))
 830  format((13x,f7.3,3x,i1,a1,i2,'/2',1x,i3))
 840  format((43x,f7.3,3x,i1,a1,i2,'/2',1x,i3))
      end
c
c-----------------------------------------------------------------------====
c
      subroutine signn(e0,sigpp,sigpn)
c
      e=e0
      if(e.gt.650.) e=650.
      ei=1./e
      if(e.gt.111.511) then
        sigpn=28.1451+ei*(1431.74+ei*297627.)
       else
        if(e.gt.38.22) then
          sigpn=4.923+ei*(5673.25+ei*113412.)
         else
          sigpn=-1.86+e*(0.09415+e*1.306e-4)
          sigpn=3./(1.206*e+sigpn*sigpn)
          sigpn=sigpn+1./(1.206*e+(.4233+e*0.13)**2)
          sigpn=3141.59*sigpn
         endif
       endif
      if(e.gt.158.555) then
        sigpp=11.7386+0.0189*e+1362.11*ei
       else
        if(e.gt.42.738) then
          sigpp=17.8465+ei*(454.414+ei*65760.5)
         else
          if(e.gt.10.54) then
            sigpp=-2.0331+ei*(2690.66+ei*6498.86)
           else
            ei=1./(e+0.62)
            sigpp=-229.51+ei*(6596.46-ei*3920.61)
            sigpp=e*ei*sigpp
           endif
         endif
       endif
      return
      end
c
c-----------------------------------------------------------------------====
c
      subroutine precomb(zp,ap,zt,at,f0,sig0)
c
      include 'fragyld.par'
      common/nucs/cn(4,nex,kdn,kdz),nmax(kdn,kdz),za(2),acx(3),de,ip,p0
      common/hilo/yldp(5,kdn,kdz),hyld(kdn,kdz),nzn(2),fcor,ncor
      common/ees/el(kst,2),ne(kst,2),nxf
      common/ind/eps(ns,2),jg(ns,2),nx(2),izn(2),nzf(2)
      common/tind/nn(ns,2),ll(ns,2),jj(ns,2),b0,iwt
      do 20 l=1,2
        kl=1
        jl=jg(kl,l)+1
        ell=-eps(kl,l)
        kh=nx(l)
        jh=jg(kh,l)+1
        ehh=-eps(kh,l)
        do 10 j=1,nzn(l)
          ell=ell+eps(kl,l)
          jl=jl-1
          if(jl.le.0) then
            kl=kl+1
            jl=jg(kl,l)
           endif
          el(j,l)=1.001-ell
          ehh=ehh+eps(kh,l)
          jh=jh-1
          if(jh.le.0) then
            kh=kh-1
            jh=jg(kh,l)
           endif
          ne(j,l)=ehh-ell+1.999
 10     continue
 20   continue
      if(iwt.gt.0) then
        do 50 l=1,2
          jtot=0
          do 30 k=1,nx(l)
            jtot=jtot+jg(k,l)
 30       continue
          nsmax=(jtot+1)/2
          if(nsmax.gt.ns) then
            write(8,800) nsmax,ns
            stop
           endif
          jd=2
          if(2*nsmax.gt.jtot) jd=1
          koc=nx(l)
          joc=jg(koc,l)
          do 40 k=nsmax,1,-1
            nn(k,l)=nn(koc,l)
            ll(k,l)=ll(koc,l)
            jj(k,l)=jj(koc,l)
            eps(k,l)=eps(koc,l)
            jg(k,l)=jd
            joc=joc-jd
            if(joc.le.0.and.k.gt.1) then
              koc=koc-1
              joc=jg(koc,l)
             endif
            jd=2
 40       continue
          nx(l)=nsmax
 50     continue
       endif
      iz=izn(1)+1
      neu=izn(2)
      do 80 jz=1,nzn(1)
        iz=iz-1
        ia=iz+neu+1
        do 70 jn=1,nzn(2)
          ia=ia-1
          yldp(1,jn,jz)=0.
          yldp(2,jn,jz)=0.
          yldp(3,jn,jz)=0.
          yldp(4,jn,jz)=0.
          yldp(5,jn,jz)=0.
          hyld(jn,jz)=0.
          nmax(jn,jz)=0
          do 60 n=1,nex
            cn(1,n,jn,jz)=0.
 60       continue
          cn(2,1,jn,jz)=
     1        sqrt(1.+f0*((ap+at-ia)*iz*(zp+zt-iz)/(ia*at)-zp*zt/ap))
          cn(3,1,jn,jz)=ia*(ap-ia)*sig0/(ap-1.)
 70     continue
 80   continue
      return
 800  format(' In PRECOMB, number of states, NMAX=',i2,' exceeds array',
     1 ' dimension, NS=',i2)
      end
c
c-----------------------------------------------------------------------====
c
      subroutine thprep(ap,zt,at)
c
      include 'fragyld.par'
      dimension rm(4)
      common/dense/dc(4),w(4),r(4),a(4),itp(4),ird,iwrt
      common/thick/ipzn,jtmax,jpmax,jbmax,dr,sigpp,sigpn,stuff(nst)
      xr=at**(1./3.)
      xp=ap**(1./3.)
      aamin=10.
      do 10 i=1,4
        if(i.gt.2) xr=xp
        if(r(i).gt.0.) then
          r(i)=xr*r(i)
         else
          r(i)=-r(i)
         endif
        if(itp(i).gt.1) then
          rm(i)=r(i)*sqrt(2.31*npr+alog(1.+3.*npr))
          aamin=amin1(aamin,r(i)/3.)
         else
          rm(i)=r(i)+2.31*npr*a(i)
          w(i)=w(i)/(r(i)*r(i))
          aamin=amin1(aamin,a(i))
         endif
 10   continue
      rtmax=amax1(rm(1),rm(2))
      rpmax=amax1(rm(3),rm(4))
      rmax=amax1(rtmax,rpmax)
      dr=aamin*(10.**(2-npr)/rmax)**(0.25)
      jtmax=2*(int(rtmax/dr+2.)/2)
      jpmax=2*(int(rpmax/dr+2.)/2)
      jbmax=jtmax+jpmax-2.31*npr
      jmax=max0(jtmax,jpmax)+1
      if(jmax.gt.ntt) then
        write(8,800) dr,jmax
        stop
       endif
      if(iwrt.ne.0) iwrt=(jbmax-1)/40+1
      call tarth(zt,at,rm)
      call proth(rm)
      return
 800  format(' for dr= ',f4.2,', ntt must be at least ',i4,'.')
      end
c
c-----------------------------------------------------------------------====
c
      subroutine tarth(zt,at,rm)
c
      include 'fragyld.par'
      character*7 atip(2)
      dimension rm(4),zx(2)
      common/dense/dc(4),w(4),r(4),a(4),itp(4),ird,iwrt
      common/thick/ipzn,jtmax,jpmax,jbmax,dr,sigpp,sigpn,tt(ntt,2)
     1                                               ,tp(ntt,ns,2)
      data atip/' proton','neutron'/
      data pi4/12.56637/
      zx(1)=zt
      zx(2)=at-zt
      do 50 l=1,2
        if(itp(l).gt.1) then
          dc(l)=0.25*pi4*(1.+1.5*a(l))*r(l)**3
          dc(l)=zx(l)/dc(l)
          e=1./(r(l)*r(l))
          s=-dr
          do 10 j=1,jtmax
            s=s+dr
            s2=e*s*s
            tt(j,l)=dc(l)*r(l)*exp(-s2)*(1.+a(l)*(0.5+s2))
 10       continue
         else
          s=0.
          dc(l)=0.
          e=exp(-r(l)/a(l))
          de=exp(dr/a(l))
          dws=-2.*pi4*dr/3.
          ws=-2.*dws
          kr=rm(l)/dr+2.
          kr=2*(kr/2)
          rm(l)=rm(l)**2
          do 20 k=1,kr
            s=s+dr
            s2=s*s
            e=e*de
            dc(l)=dc(l)+ws*s2*(1.+w(l)*s2)/(1.+e)
            ws=ws+dws
            dws=-dws
 20       continue
          dc(l)=zx(l)/dc(l)
          s=-dr
          do 40 j=1,jtmax
            s=s+dr
            s2=s*s
            z=0.
            ws=dr*dc(l)/3.
            tt(j,l)=ws*(1.+w(l)*s2)/(1.+exp((s-r(l))/a(l)))
            dws=-2.*ws
            ws=-2.*dws
            mr=sqrt(amax1(rm(l)-s2,0.))/dr+2.
            mr=2*(mr/2)
            do 30 m=1,mr
              z=z+dr
              x2=s2+z*z
              x=sqrt(x2)
              tt(j,l)=tt(j,l)+ws*(1.+w(l)*x2)/(1.+exp((x-r(l))/a(l)))
              ws=ws+dws
              dws=-dws
 30         continue
            tt(j,l)=2.*tt(j,l)
 40       continue
         endif
        tt(jtmax+1,l)=0.
        if(iwrt.ne.0) then
          write(8,800) atip(l),iwrt*dr
          write(8,810) (tt(j,l),j=1,jtmax,iwrt)
         endif
 50   continue
      return
 800  format(/,'   Target   ',a7,' thickness (1/fm**2)  db=',f4.2,' fm')
 810  format((10f8.4))
      end
c
c-----------------------------------------------------------------------====
c
      subroutine proth(rm)
c
      include 'fragyld.par'
      character*7 atip(2)
      dimension wf(ntt),rm(4)
      common/dense/dc(4),w(4),r(4),a(4),itp(4),ird,iwrt
      common/ind/eps(ns,2),jg(ns,2),nx(2),izn(2),nzf(2)
      common/tind/nn(ns,2),ll(ns,2),jj(ns,2),b0,iwt
      common/thick/ipzn,jtmax,jpmax,jbmax,dr,sigpp,sigpn,tt(ntt,2)
     1                                               ,tp(ntt,ns,2)
      data pi4/12.56637/,sqpi/1.772454/
      data atip/' proton','neutron'/
      if(iwt.lt.0) then
        do 80 l=3,4
          lp=l-2
          if(itp(l).gt.1) then
            dc(l)=0.25*pi4*(1.+1.5*a(l))*r(l)**3
            dc(l)=1./dc(l)
            e=1./(r(l)*r(l))
            s=-dr
            do 20 j=1,jpmax
              s=s+dr
              s2=e*s*s
              tpp=dc(l)*r(l)*exp(-s2)*(1.+a(l)*(0.5+s2))
              do 10 k=1,nx(lp)
                tp(j,k,lp)=tpp
 10           continue
 20         continue
           else
            s=0.
            dc(l)=0.
            e=exp(-r(l)/a(l))
            de=exp(dr/a(l))
            dws=-2.*pi4*dr/3.
            ws=-2.*dws
            kr=rm(l)/dr+2.
            kr=2*(kr/2)
            rm(l)=rm(l)**2
            do 30 k=1,kr
              s=s+dr
              s2=s*s
              e=e*de
              dc(l)=dc(l)+ws*s2*(1.+w(l)*s2)/(1.+e)
              ws=ws+dws
              dws=-dws
 30         continue
            dc(l)=1./dc(l)
            s=-dr
            do 60 j=1,jpmax
              s=s+dr
              s2=s*s
              z=0.
              ws=dr*dc(l)/3.
              tpp=ws*(1.+w(l)*s2)/(1.+exp((s-r(l))/a(l)))
              dws=-2.*ws
              ws=-2.*dws
              mr=sqrt(amax1(rm(l)-s2,0.))/dr+2.
              mr=2*(mr/2)
              do 40 m=1,mr
                z=z+dr
                x2=s2+z*z
                x=sqrt(x2)
                tpp=tpp+ws*(1.+w(l)*x2)/(1.+exp((x-r(l))/a(l)))
                ws=ws+dws
                dws=-dws
 40           continue
              tpp=2.*tpp
              do 50 k=1,nx(lp)
                tp(j,k,lp)=tpp
 50           continue
 60         continue
           endif
          do 70 k=1,nx(lp)
            tp(jpmax+1,k,lp)=0.
 70       continue
          if(iwrt.ne.0) then
            write(8,810) atip(lp),iwrt*dr
            write(8,820) (tp(j,1,lp),j=1,jpmax,iwrt)
           endif
 80     continue
       else
        do 120 l=1,2
          rm(l+2)=rm(l+2)**2
          if(iwrt.ne.0) write(8,810) atip(l),iwrt*dr
          jold=-2
          do 110 k=1,nx(l)
            if(jj(k,l).ne.jold) then
              jold=jj(k,l)
              mm=0
             else
              mm=mm+1
             endif
            if(ird.gt.0) then
              read(5,*) nwf,dx
              if(nwf.gt.ntt-2) then
                write(8,800) nwf,ntt-2
                stop
               endif
              read(5,*) (wf(n),n=1,nwf)
              wf(nwf+1)=0.
              wf(nwf+2)=0.
              wnorm=1./pi4
             else
              wnorm=(hon(nn(k,l),ll(k,l)))**2/(sqpi*pi4*b0**3)
             endif
            s=-dr
            do 100 j=1,jpmax
              s=s+dr
              s2=s*s
              ws=wnorm*dr/3.
              if(ird.gt.0) then
                rn=s/dx
                n=rn
                rn=rn-n
                if(n.gt.nwf) n=nwf
                wfr=(1.-rn)*wf(n+1)+rn*wf(n+2)
               else
                wfr=howf(nn(k,l),ll(k,l),s/b0)
               endif
              if(iwt.gt.0) then
                wfr=wfr**2*wlm(0.,mm,ll(k,l),jj(k,l))
               else
                wfr=wfr**2
               endif
              tpp=ws*wfr
              z=0.
              dws=-2.*ws
              ws=-2.*dws
              mr=sqrt(amax1(rm(l+2)-s2,0.))/dr+2.
              mr=2*(mr/2)
              do 90 m=1,mr
                z=z+dr
                rr=sqrt(s2+z*z)
                if(ird.gt.0) then
                  rn=rr/dx
                  n=rn
                  rn=rn-n
                  if(n.gt.nwf) n=nwf
                  wfr=(1.-rn)*wf(n+1)+rn*wf(n+2)
                 else
                  wfr=howf(nn(k,l),ll(k,l),rr/b0)
                 endif
                if(iwt.gt.0) then
                  wfr=wfr**2*wlm(z/rr,mm,ll(k,l),jj(k,l))
                 else
                  wfr=wfr**2
                 endif
                tpp=tpp+ws*wfr
                ws=ws+dws
                dws=-dws
 90           continue
              tp(j,k,l)=2.*tpp
 100        continue
            tp(jpmax+1,k,l)=0.
            if(iwrt.ne.0) then
              write(8,830) k
              write(8,820) (tp(j,k,l),j=1,jpmax,iwrt)
             endif
 110      continue
 120    continue
       endif
      return
 800  format(' In PROTH, the variable NWF=',i3,' exceeded NTT-2=',i3)
 810  format(/,' Projectile ',a7,' thickness (1/fm**2)  db=',f4.2,' fm')
 820  format((10f8.4))
 830  format('  Level ',i2)
      end
c
c-----------------------------------------------------------------------====
c
      function howf(n,l,x)
c
      howf=0.
      if(x.le.0.) return
      howf=1.
      xx=-x*x
      nn=n-1
      if(nn.gt.0) then
        nx=nn
        al=(l+3.)/2.
        dh=1.
        do 10 k=1,nn
          dh=dh*nx*xx/(k*al)
          howf=howf+dh
          nx=nx-1
          al=al+1.
 10     continue
       endif
      howf=howf*x**(l/2)*exp(0.5*xx)
      return
      end
c
c-----------------------------------------------------------------------====
c
      function hon(n,l)
c
      dimension fac(4),facac(10)
      data fac/0.00,0.00,6.9314718055995e-01,1.7917594692281e+00/
      data facac/0.00
     1    ,1.0986122886681e+00,2.7080502011022e+00,4.6539603501574e+00
     2    ,6.8511849274933e+00,9.2490802002923e+00,1.1814029557753e+01
     3    ,1.4522079758856e+01,1.7355293102912e+01,2.0299732082078e+01/
      ll=l/2
      hon=(ll-n+3.)*fac(3)+facac(ll+n)-fac(n)
      hon=exp(0.5*hon-facac(ll+1))
      return
      end
c
c-----------------------------------------------------------------------====
c
      function wlm(cost,mm,ll,jj)
c
      dimension wlnrm(11)
      data wlnrm/1.00000000,0.70710678,0.61237244,0.55901699,0.52291252,
     1           0.49607837,0.47495888,0.45768183,0.44314852,0.43066296,
     2           0.41975833/
      l=ll/2
      mp=mm+1
      sint=amax1(sqrt(1.-cost*cost),1.0e-20)
      sn=wlnrm(mm+1)*sint**mm
      spn=0.
      if(l.gt.mm) then
        so=sn
        an=sqrt(2.*mm+1.)
        sn=an*cost*so
        spn=wlnrm(mp+1)*sint**mp
        if(l.gt.mp) then
          spo=0.
          apn=0.
          xl=2.*mp-1.
          do 10 lx=mp+1,l
            xl=xl+2.
            ao=an
            an=sqrt(float(lx*lx-mm*mm))
            t=(xl*cost*sn-ao*so)/an
            so=sn
            sn=t
            apo=apn
            apn=sqrt(float(lx*lx-mp*mp))
            t=(xl*cost*spn-apo*spo)/apn
            spo=spn
            spn=t
 10       continue
         endif
       endif
      if(jj.gt.ll) then
        wlm=(l+mp)*sn*sn+(l-mm)*spn*spn
       else
        wlm=(l-mm)*sn*sn+(l+mp)*spn*spn
       endif
      return
      end
c
c-----------------------------------------------------------------------====
c
      subroutine emdiss(zp,ap,zt,at,ep)
c
      data am/939.0/,hc/197.3/,pi/3.14159/
      cap=ap**(1./3.)
      cat=at**(1./3.)
c   Guesstimate
c      b0=1.25*(cap+cat)
c   Benesh-Cook-Vary
      b0=1.34*(cap+cat-0.75*(1./cap+1./cat))
c   Kox et al.
c      b0=1.1*(cap+cat+1.85*cap*cat/(cap+cat)-1.9)
c
      egr=80./cap
      x0=egr/(hc*sqrt((1.+ep/am)**2-1.))
      bmax=amax1(1./x0,b0)
      deg=2.*zp*(ap-zp)*(1.44*zt)**2/(ap*am*(1.-(am/(am+ep))**2))
      sig=20.*pi*deg*alog(bmax/b0)/egr
      eav=deg/b0**2
      avn=eav/egr
      write(8,800) egr,sig,b0,eav,avn
      return
 800  format(//,29x,'Coulomb Dissociation',/,4x,'E1 Res.E =',f6.2,
     1 ' MeV',28x,'Diss.Yld = ',f8.2,' mb',/,8x,'At b = ',f5.2,
     2 ' Fm',6x,'Ave.Ex.E =',f6.2,' MeV',6x,'Ave.N = ',
     3 f10.4)
      end
c
c-----------------------------------------------------------------------====
c
      subroutine abrade
c
      include 'fragyld.par'
      dimension ave(5,40),pav(4)
      common/thick/ipzn,jtmax,jpmax,jbmax,dr,sigxx(2),stuff(nst)
      data pi2/6.283186/
      jav=(jbmax-1)/40+1
      nav=0
      b=0.
      dwb=-20.*pi2*dr/3.
      wb=-2.*dwb
      do 10 jb=1,jbmax
      b=b+dr
      call pzn(jb,pav)
c      write(8,'(5f10.5)') b,pav
      if(jb.eq.jav*(jb/jav)) then
            nav=nav+1
            ave(1,nav)=b
            ave(2,nav)=pav(1)
            ave(3,nav)=pav(2)
            ave(4,nav)=pav(1)+pav(2)
            ave(5,nav)=sqrt(pav(3)+pav(4))
         endif
      call combos(b,wb*b,pav)
      wb=wb+dwb
      dwb=-dwb
 10   continue
      nav=(nav-1)/2+1
      write(8,800)
      write(8,810) ((ave(j,i),j=1,5),(ave(j,i+nav),j=1,5),i=1,nav)
      return
 800  format(///,22x,'Average Number of Transmitted Nucleons',/,
     1  2(2x,'b (Fm)',3x,'Prot   Neut',3x,'Total',2x,'Width',2x))
 810  format(2(1x,f6.2,1x,2f7.2,1x,2f7.2,2x))
      end
c
c-----------------------------------------------------------------------====
c
      subroutine pzn(jb,pav)
c
      include 'fragyld.par'
      dimension dp(ns,2),pav(4),tzn(2)
      common/pees/psp(ns,2),plsp(ns,2),pmom(7,kst,2)
      common/ind/eps(ns,2),jg(ns,2),nx(2),izn(2),nzf(2)
      common/thick/ipzn,jtmax,jpmax,jbmax,dr,sigpp,sigpn,tt(ntt,2)
     1                                               ,tp(ntt,ns,2)
      data pi/3.141593/
      do 30 l=1,2
        pav(l)=0.
        pav(l+2)=0.
        do 10 i=1,7
          pmom(i,1,l)=0.
          pmom(i,2,l)=0.
 10     continue
        do 20 k=1,nx(l)
          plsp(k,l)=0.
 20     continue
 30   continue
      b=dr*jb
      b2=b*b
      dws=-2.*dr/3.
      ws=-2.*dws
      jsmin=max0(jb-jpmax,0)+1
      jsmax=min0(jtmax,jb+jpmax-1)
      jsmax=2*((jsmax-jsmin+1)/2)+jsmin-1
      s=dr*(jsmin-1)
      do 140 js=jsmin,jsmax
        s=s+dr
        bs0=b2+s*s
        bs1=2.*b*s
        th=0.
        dth=dr/amin1(b,s)
        nth=2*(int(pi/dth+2.)/2)
        dth=pi/nth
        wt=dth/3.
        dwt=-2.*wt
        jbs=iabs(jb-js)+1
        do 50 l=1,2
          do 40 k=1,nx(l)
            dp(k,l)=wt*tp(jbs,k,l)
 40       continue
 50     continue
        wt=4.*wt
        do 80 nt=1,nth
          th=th+dth
          bs=sqrt(bs0-bs1*cos(th))
          dbs=bs/dr
          jbs=dbs
          dbs=dbs-jbs
          jbs=jbs+1
          if(jbs.gt.jpmax) go to 110
          do 70 l=1,2
            do 60 k=1,nx(l)
              ddp=(1.-dbs)*tp(jbs,k,l)+dbs*tp(jbs+1,k,l)
              dp(k,l)=dp(k,l)+wt*ddp
 60         continue
 70       continue
          wt=wt+dwt
          dwt=-dwt
 80     continue
        wt=0.25*wt
        do 100 l=1,2
          do 90 k=1,nx(l)
            ddp=(1.-dbs)*tp(jbs,k,l)+dbs*tp(jbs+1,k,l)
            dp(k,l)=dp(k,l)-wt*ddp
 90       continue
 100    continue
 110    jt=js+1
        tzn(1)=sigpp*tt(jt,1)+sigpn*tt(jt,2)
        tzn(2)=sigpn*tt(jt,1)+sigpp*tt(jt,2)
        do 130 l=1,2
          if(ipzn.gt.0) then
            tzn(l)=ws*s*(1.-exp(-amin1(tzn(l),65.)))
           else
            tzn(l)=ws*s*tzn(l)
          endif
          do 120 k=1,nx(l)
            plsp(k,l)=plsp(k,l)+dp(k,l)*tzn(l)
 120      continue
 130    continue
        ws=ws+dws
        dws=-dws
 140  continue
c      xint=0.
      do 200 l=1,2
        do 160 k=1,nx(l)
          if(ipzn.gt.0) then
            dp(k,l)=amax1(amin1(1.-2.*plsp(k,l),1.),0.)**ipzn
           else
c            xint=xint+2.*jg(k,l)*plsp(k,l)
            dp(k,l)=exp(-amin1(2.*plsp(k,l),20.))
           endif
          pav(l)=pav(l)+jg(k,l)*dp(k,l)
          pav(l+2)=pav(l+2)+jg(k,l)*dp(k,l)*(1.-dp(k,l))
          ddp=amax1(dp(k,l),1.0e-7)
          pmom(1,1,l)=pmom(1,1,l)+jg(k,l)*alog(ddp)
          ddp=amax1(1.-dp(k,l),1.0e-7)/ddp
          psp(k,l)=ddp
          plsp(k,l)=alog(ddp)
          ex=jg(k,l)
          do 150 i=1,7
            pmom(i,2,l)=pmom(i,2,l)+ex*ddp
            ex=ex*eps(k,l)
 150      continue
 160    continue
        do 170 i=2,7
          pmom(i,2,l)=pmom(i,2,l)/pmom(1,2,l)
 170    continue
        ddp=pmom(1,2,l)*exp(pmom(1,1,l))
        pmom(1,2,l)=1.
        ex=pmom(2,2,l)
        do 190 i=2,7
          c=1.
          do 180 k=2,i
            c=c*ex*(i-k+1.)/(k-1.)
            pmom(i,2,l)=pmom(i,2,l)-c*pmom(i-k+1,2,l)
 180      continue
 190    continue
        pmom(1,2,l)=ddp
        pmom(2,2,l)=ex
        pmom(5,2,l)=pmom(5,2,l)-3.*pmom(3,2,l)**2
        pmom(6,2,l)=pmom(6,2,l)-10.*pmom(4,2,l)*pmom(3,2,l)
        pmom(7,2,l)=pmom(7,2,l)-10.*pmom(4,l,2)**2
     1             -15.*pmom(3,2,l)*(pmom(5,2,l)+pmom(3,2,l)**2)
 200  continue
c      write(8,'(2f10.4)') b,xint
c      write(9,'(f7.3)') b
c      write(9,'(10f7.4)') ((dp(k,l),k=1,nx(l)),l=1,2)
      return
      end
c
c-----------------------------------------------------------------------====
c
      subroutine combos(b,dsig,pav)
c
      include 'fragyld.par'
      dimension e(kst),rr(kst),pl(kst),gs(kst),js(kst),pav(4)
      dimension yldx(kst,2),cnx(nex,kst,2),nzz(2),nzm(2)
      common/ind/eps(ns,2),jg(ns,2),nx(2),izn(2),nzf(2)
      common/nucs/cn(4,nex,kdn,kdz),nmax(kdn,kdz),za(2),acx(3),de,ip,p0
      common/hilo/yldp(5,kdn,kdz),hyld(kdn,kdz),nzn(2),fcor,ncor
      common/pees/psp(ns,2),plsp(ns,2),pmom(7,kst,2)
      common/ees/el(kst,2),ne(kst,2),nxf
      yldp(1,1,1)=yldp(1,1,1)-dsig
      do 100 l=1,2
        pav(l)=izn(l)-pav(l)
        p=8.*sqrt(pav(l+2))+0.6
        nzz(l)=min0(max0(int(pav(l)+p+1.),2),nzn(l))
        nzm(l)=max0(int(pav(l)-p+1.),1)
        if(nzm(l).gt.nzz(l)) return
        do 20 j=1,nzz(l)
          yldx(j,l)=0.
          nhi=min0(nex,ne(j,l))
          if(j.gt.nzf(l)) nhi=min0(nhi,nxf)
          do 10 i=1,nhi
            cnx(i,j,l)=0.
 10       continue
 20     continue
        j=1
        js(j)=1
        e(j)=0.
        gs(j)=0.
        rr(j)=1.
        pl(j)=pmom(1,1,l)
        jj=j
 30     j=j+1
        js(j)=js(j-1)
 40     er=e(jj)+el(jj,l)
        ie=er
        if(jj.le.nzf(l).or.ie.le.nxf) then
          er=er-ie
          rx=rr(jj)*exp(amax1(pl(jj),-80.))
          if(ie.le.nex) then
            cnx(ie,jj,l)=cnx(ie,jj,l)+rx*(1.-er)
           endif
          if(ie.lt.nex) cnx(ie+1,jj,l)=cnx(ie+1,jj,l)+rx*er
          yldx(jj,l)=yldx(jj,l)+rx
         else
          js(j)=nx(l)+1
         endif
 50     if(js(j).le.nx(l)) go to 70
        j=j-1
        if(j.le.1) go to 80
 60     js(j)=js(j)+1
        go to 50
 70     jj=j
        oc=0.
        if(js(j).eq.js(j-1)) oc=gs(j-1)
        gs(j)=oc+1.
        rr(j)=rr(j-1)*(jg(js(j),l)-oc)/gs(j)
        if(rr(j).lt.0.5) go to 60
        e(j)=e(j-1)+eps(js(j),l)
        pl(j)=pl(j-1)+plsp(js(j),l)
        if(j.lt.nzz(l)) go to 30
        js(j)=js(j)+1
        go to 40
 80     if(nzm(l).le.nzz(l)) then
          do 90 j=nzm(l),nzz(l)
            if(j.gt.2) call moms(j,l,pav(l))
            pmom(2,j,l)=pmom(2,j,l)+el(j,l)-1.001
            if(j.gt.nzf(l).and.ne(j,l).gt.nxf) yldx(j,l)=pmom(1,j,l)
 90       continue
         endif
 100  continue
      do 140 jz=nzm(1),nzz(1)
        do 130 jn=nzm(2),nzz(2)
          nhi=min0(ne(jz,1)+ne(jn,2)-1,nex)
          if(jz.gt.nzf(1).or.jn.gt.nzf(2)) nhi=min0(nhi,nxf)
          yld=0.
          do 120 i=1,nhi
            p=0.
            do 110 k=max0(i-ne(jz,1),0)+1,min0(ne(jn,2),i)
              p=p+cnx(k,jn,2)*cnx(i-k+1,jz,1)
 110        continue
            yld=yld+p
            cn(1,i,jn,jz)=cn(1,i,jn,jz)+p*dsig
 120      continue
          p=yldx(jn,2)*yldx(jz,1)*dsig
          yldp(1,jn,jz)=yldp(1,jn,jz)+p
          yldp(2,jn,jz)=yldp(2,jn,jz)+p*b
          yldp(3,jn,jz)=yldp(3,jn,jz)+p*b*b
          yldp(4,jn,jz)=yldp(4,jn,jz)+p*(pmom(2,jz,1)+pmom(2,jn,2))
          yldp(5,jn,jz)=yldp(5,jn,jz)
     1    +p*(pmom(3,jz,1)+pmom(3,jn,2)+(pmom(2,jz,1)+pmom(2,jn,2))**2)
          yld=yld*dsig
          if(jz.gt.nzf(1).or.jn.gt.nzf(2)) call cndis(jn,jz,p,yld)
          hyld(jn,jz)=hyld(jn,jz)+p-yld
 130    continue
 140  continue
      return
      end
c
c-----------------------------------------------------------------------====
c
      subroutine moms(j,l,df0)
c
      include 'fragyld.par'
      dimension dm(7)
      common/ind/eps(ns,2),jg(ns,2),nx(2),izn(2),nzf(2)
      common/pees/psp(ns,2),plsp(ns,2),pmom(7,kst,2)
      data rtpi/-0.5723649/
      x=(j-1.)/df0
      do 20 it=1,10
        df=j-1.
        ddf=0.
        dm(1)=0.
        do 10 n=1,nx(l)
          p=psp(n,l)*x
          p=p/(1.+p)
          df=df-jg(n,l)*p
          ddf=ddf+jg(n,l)*p*(1.-p)
          dm(1)=dm(1)+jg(n,l)*eps(n,l)*p*(1.-p)
 10     continue
        dx=-df/ddf
        if(abs(df).lt.1.0e-5.and.abs(dx).lt.1.0e-4) go to 30
        x=x*exp(amin1(-amin1(dx,3.),3.))
 20   continue
      write(8,'('' No convergence in MOMS'',2i5,3f15.5)') l,j-1,df0,x,df
 30   pmom(1,j,l)=rtpi+(1.-j)*alog(x)-0.5*alog(2.*ddf)
      dm(1)=dm(1)/ddf
      do 40 i=2,7
        pmom(i,j,l)=0.
        dm(i)=0.
 40   continue
      do 50 n=1,nx(l)
        p=psp(n,l)*x
        pmom(1,j,l)=pmom(1,j,l)+jg(n,l)*alog((1.+p)/(1.+psp(n,l)))
        p=p/(1.+p)
        pmom(2,j,l)=pmom(2,j,l)+jg(n,l)*eps(n,l)*p
        ex=eps(n,l)-dm(1)
        pmom(3,j,l)=pmom(3,j,l)+jg(n,l)*ex*ex*p*(1.-p)
        df=jg(n,l)*p*(1.-p*(3.-p*2.))
        dm(5)=dm(5)+df
        df=ex*df
        dm(6)=dm(6)+df
        df=ex*df
        dm(2)=dm(2)+df
        df=ex*df
        pmom(4,j,l)=pmom(4,j,l)+df
        ee=jg(n,l)*ex*ex
        df=ee*p*(1.-p*(7.-p*(12.-p*6.)))
        dm(7)=dm(7)+df
        df=ex*df
        dm(3)=dm(3)+df
        df=ex*df
        pmom(5,j,l)=pmom(5,j,l)+df
        ee=ex*ex*ee
        df=p*(1.-p*(15.-p*(50.-p*(60.-p*24.))))
        dm(4)=dm(4)+ee*df
        ee=ex*ee
        pmom(6,j,l)=pmom(6,j,l)+ee*df
        ee=ex*ee
        pmom(7,j,l)=pmom(7,j,l)
     1             +ee*p*(1.-p*(31.-p*(180.-p*(390.-p*(360.-p*120.)))))
 50   continue
      dm(2)=-dm(2)/ddf
      dm(3)=-(dm(3)+3.*dm(2)*dm(6))/ddf
      dm(4)=-(dm(4)+3.*dm(2)*(2.*dm(7)+dm(2)*dm(5))+4.*dm(3)*dm(6))/ddf
      pmom(1,j,l)=exp(pmom(1,j,l))
      pmom(5,j,l)=pmom(5,j,l)-3.*ddf*dm(2)**2
      pmom(6,j,l)=pmom(6,j,l)-5.*dm(2)*(2.*ddf*dm(3)+3.*dm(6)*dm(2))
      pmom(7,j,l)=pmom(7,j,l)-10.*dm(3)*(ddf*dm(3)+6.*dm(2)*dm(6))
      pmom(7,j,l)=pmom(7,j,l)
     1           -15.*dm(2)*(ddf*dm(4)+dm(2)*(3.*dm(7)+2.*dm(2)*dm(5)))
      return
      end
c
c-----------------------------------------------------------------------====
c
      subroutine cndis(jn,jz,p,yld)
c
      include 'fragyld.par'
      dimension cc(7)
      common/nucs/cn(4,nex,kdn,kdz),nmax(kdn,kdz),za(2),acx(3),de,ip,p0
      common/pees/psp(ns,2),plsp(ns,2),pmom(7,kst,2)
      common/ees/el(kst,2),ne(kst,2),nxf
      nhi=min0(ne(jz,1)+ne(jn,2),nex)
      if(nxf.ge.nhi) return
      cyld=p-yld
      if(cyld.lt.1.0e-6) return
      cc(1)=1.
      cc(2)=pmom(2,jz,1)+pmom(2,jn,2)
      ee=pmom(3,jz,1)+pmom(3,jn,2)
      cc(3)=sqrt(ee)
      do 10 i=4,7
        ee=ee*cc(3)
        cc(i)=(pmom(i,jz,1)+pmom(i,jn,2))/ee
 10   continue
      cc(7)=cc(7)+10.*cc(4)**2
      ee=nxf
      roo=rowi(ee,cc)
      anorm=rowi(real(ne(jz,1)+ne(jn,2)),cc)-roo
      cc(1)=cyld/anorm
      roo=cc(1)*roo
      yld=yld-roo
      do 20 i=nxf+1,nhi
        ee=ee+1.
        ron=rowi(ee,cc)
        cn(1,i,jn,jz)=cn(1,i,jn,jz)+ron-roo
        roo=ron
 20   continue
      yld=yld+roo
      return
      end
c
c-----------------------------------------------------------------------====
c
      function rowi(ee,cc)
c
      dimension cc(7),he(7)
      data rtpi/-0.5723649/,rt2/0.7071067/
      rowi=0.
      ye=(ee-cc(2))/cc(3)
      he(2)=1.
      he(3)=ye
      do 10 i=4,7
        he(i)=(ye*he(i-1)-he(i-2))/(i-2)
        rowi=rowi+cc(i)*he(i)/(i-1)
 10   continue
      arg=rtpi-0.5*ye*ye
      rowi=cc(1)*(0.5*erf(rt2*ye)-rowi*exp(arg))
      return
      end
c
c-----------------------------------------------------------------------====
c
      subroutine ablate(inc)
c
      include 'fragyld.par'
      common/nucs/cn(4,nex,kdn,kdz),nmax(kdn,kdz),za(2),acx(3),de,ip,p0
      common/hilo/yldp(5,kdn,kdz),hyld(kdn,kdz),nzn(2),fcor,ncor
      data hc/197.3/
c      data be/8.0/
      yldp(1,1,1)=-yldp(1,1,1)
      do 30 jz=1,nzn(1)
        jzn=1
        if(jz.eq.1) jzn=2
        do 20 jn=jzn,nzn(2)
          if(yldp(1,jn,jz).gt.prec) then
            sigp=yldp(2,jn,jz)/yldp(1,jn,jz)
            sigp=yldp(3,jn,jz)/yldp(1,jn,jz)-sigp**2
            sigp=0.25*hc*hc/sigp
c            nxx=min0(int(be*(za(2)-jn-jz+2.)/de+2.),nex)
            nxx=nex
            vf=cn(2,1,jn,jz)
            sigh=cn(3,1,jn,jz)
            maxn=0
            gcor=1.-fcor
            do 10 n=nxx,1,-1
              nnc=n+ncor
              if(cn(1,n,jn,jz).gt.prec.and.maxn.eq.0) maxn=nnc
              if(nnc.le.nex) then
                if(nnc.le.nxx) then
                  cn(1,nnc,jn,jz)=cn(1,nnc,jn,jz)+gcor*cn(1,n,jn,jz)
                  cn(2,nnc,jn,jz)=cn(2,nnc,jn,jz)+vf*cn(1,n,jn,jz)
                  cn(3,nnc,jn,jz)=cn(3,nnc,jn,jz)+sigh*cn(1,n,jn,jz)
                  cn(4,nnc,jn,jz)=cn(4,nnc,jn,jz)+sigp*cn(1,n,jn,jz)
                 else
                  cn(1,nnc,jn,jz)=gcor*cn(1,n,jn,jz)
                  cn(2,nnc,jn,jz)=vf*cn(1,n,jn,jz)
                  cn(3,nnc,jn,jz)=sigh*cn(1,n,jn,jz)
                  cn(4,nnc,jn,jz)=sigp*cn(1,n,jn,jz)
                 endif
               else
                hyld(jn,jz)=hyld(jn,jz)+gcor*cn(1,n,jn,jz)
               endif
              cn(1,n,jn,jz)=fcor*cn(1,n,jn,jz)
              cn(2,n,jn,jz)=vf*cn(1,n,jn,jz)
              cn(3,n,jn,jz)=sigh*cn(1,n,jn,jz)
              cn(4,n,jn,jz)=sigp*cn(1,n,jn,jz)
 10         continue
            nmax(jn,jz)=min0(maxn,nex)
            hyld(jn,jz)=100.*hyld(jn,jz)/yldp(1,jn,jz)
            if(maxn.gt.nex.and.hyld(jn,jz).gt.1.) write(8,800)
     1                        hyld(jn,jz),za(1)-jz+1.,za(2)-jn-jz+2.
           else
            cn(2,1,jn,jz)=0.
            cn(3,1,jn,jz)=0.
            cn(4,1,jn,jz)=0.
           endif

          if(jn.le.0.and.jz.le.0) then
            write(8,820) za(1)-jz+1.,za(2)-jn-jz+2.,nmax(jn,jz)
            write(8,850) (cn(1,n,jn,jz),n=1,nmax(jn,jz))
          endif
 820  format(/,' p0(n)      z=',f4.0,3x,'a=',f4.0,3x,'nmax= ',i3)
 850  format(10f9.4)

 20     continue
 30   continue
      do 50 jz=1,nzn(1)
        jzn=1
        if(jz.eq.1) jzn=2
        do 40 jn=jzn,nzn(2)
          if(nmax(jn,jz).gt.0) call decay(jn,jz,inc)
 40     continue
 50   continue
      return
 800  format(f6.1,' % of the population of the nucleus Z=',f3.0,'  A=',
     1 f4.0,  ' lies above emax.')
      end
c
c-----------------------------------------------------------------------====
c
      subroutine decay(jn,jz,inc)
c
      include 'fragyld.par'
      double precision ro(nex,kem),gam(nex),dpav
      dimension tau(2,nex,kem),nb(kem),emc(kem)
      dimension iz(19),in(19),pav(4)
      common/nucs/cn(4,nex,kdn,kdz),nmax(kdn,kdz),za(2),acx(3),de,ip,p0
c      common/pemtab/emt(kem),emd(kem,kdn,kdz)
      common/pemtab/emt(kem)
      data iz/0,0,1,2,1,1,2,2,3,3,4,2,3,4,5,3,4,5,6/
      data in/0,1,0,2,1,2,1,4,3,4,3,6,5,4,3,6,5,4,3/
      data gamin/1.0e-6/
c decay uses the e distribution and moments accumulated for each residue
c and the partial decay widths calculated in gamem and partem to update the
c e distribution of the appropriate residual nuclei.
      call gamem(jn,jz,gam,ro,tau)
      if(nmax(jn,jz).lt.1) return
      nb(1)=0
      call partem(jn,jz,inc,nb,gam,ro,tau)
      nmin0=2
      if(sngl(gam(1)/ro(1,1)).gt.gamin) nmin0=1
      do 70 i=1,inc
        jbn=jn+in(i)
        jbz=jz+iz(i)
        emc(i)=0.
        nmin=max0(nb(i)+1,nmin0)
        nrmax=nmax(jn,jz)-nb(i)+1
        do 50 n=nmax(jn,jz),nmin,-1
          nrmax=nrmax-1
          if(cn(1,n,jn,jz).lt.prec) go to 40
          do 20 l=2,4
            pav(l)=cn(l,n,jn,jz)/cn(1,n,jn,jz)
 20       continue
          dpav=cn(1,n,jn,jz)/gam(n)
          nrmr=nrmax+1
          do 30 nr=1,nrmax
            nrmr=nrmr-1
            dpb=tau(1,nrmr,i)*sngl(ro(nr,i)*dpav)
            cn(1,nr,jbn,jbz)=cn(1,nr,jbn,jbz)+dpb
            cn(2,nr,jbn,jbz)=cn(2,nr,jbn,jbz)+dpb*pav(2)
            cn(3,nr,jbn,jbz)=cn(3,nr,jbn,jbz)+dpb*(tau(2,nrmr,i)+pav(3))
            cn(4,nr,jbn,jbz)=cn(4,nr,jbn,jbz)+dpb*pav(4)
            emc(i)=emc(i)+dpb
 30       continue
 40       continue
 50     continue
        emt(i)=emt(i)+emc(i)
c       emd(i,jn,jz)=emd(i,jn,jz)+emc(i)
 70   continue
      if(sngl(gam(1)/ro(1,1)).gt.gamin) then
        do 80 l=1,4
          cn(l,1,jn,jz)=0.
 80       continue
       endif
      if(ip.gt.1) write(8,800) (i,emc(i)/p0,i=1,inc)
      return
 800  format(16x,4(2x,'br',i1,'=',f7.3))
      end
c
c-----------------------------------------------------------------------====
c
      subroutine gamem(jn,jz,gam,ro,tau)
c
      include 'fragyld.par'
      double precision ro(nex,kem),gam(nex)
      dimension tau(2,nex,kem),pav(6)
      common/nucs/cn(4,nex,kdn,kdz),nmax(kdn,kdz),za(2),acx(3),de,ip,p0
      data pi12/37.69911/,pi2i/0.159155/,hi/33.00/,amu/939.0/
      data er0/43.4/,erx0/-0.215/,gr0/0.3/,sigr0/4.0e-6/
c gamem calculates gamma emission level densities, emission factors and partial
c widths, assuming dominance of the Brink-Axel E1 mode. It also calculates
c average properties of the CN distribution (and prints them on request).
      znc=za(1)-jz+1.
      anc=za(2)-jz-jn+2.
      hnc=hi/anc**(5./3.)
      if(acx(1).le.0.) then
        aa=acn(znc,anc,pair)
       else
        aa=anc/acx(1)
        pair=2.*(1+int(znc+0.5)/2+int(anc-znc+0.5)/2)-anc
        pair=pair*acx(2)/sqrt(anc)
       endif
      ro0=aa*(hnc*aa)**1.5/pi12
      t0=100./(20.+anc)
      if(acx(3).gt.0.5) then
        egc=aa*t0
        ep=0.25*egc*egc
        egc=1./egc
        do 10 it=1,10
          sqp=sqrt(ep)
          f=1./sqp+0.5/ep-2.5/(ep+2)-egc
          df=0.5*(1./sqp+1./ep)/ep-2.5/(ep+2)**2
          dep=f/df
          if(abs(f).le.1.0e-4.and.abs(dep).lt.1.0e-4) go to 20
          ep=amax1(ep+dep,0.5*ep)
 10     continue
 20     egc=ep/aa+pair
        rogc=ro0*sqp*exp(2.*sqp-egc/t0)/(ep+2.)**2.5
       endif
      er=er0*exp(erx0*alog(anc))
      gr=gr0*er
      er=er*er
      gr=gr*gr
      sigr=sigr0*gr*anc
      do 30 l=1,6
        pav(l)=0.
 30   continue
      maxn=0
      ex=-de
      do 70 n=1,nmax(jn,jz)
        ex=ex+de
        if(acx(3).gt.0.5) then
          if(ex.gt.egc) then
            ep=aa*(ex-pair)
            sqp=sqrt(ep)
            ro(n,1)=ro0*sqp*dexp(dble(2.*sqp))/(ep+2.)**2.5
           else
            ro(n,1)=rogc*exp(ex/t0)
           endif
         else
          ep=aa*amax1(ex,1.0e-8)
          sqp=sqrt(ep)
          t=aa*(1./sqp+0.5/ep-2.5/(ep+2.))
          t=amax1(t0,1./t)
          ro(n,1)=ro0*sqp*dexp(dble(2.*sqp-pair/t))/(ep+2.)**2.5
         endif
        ex2=ex*ex
        xtau=sigr*ex2*ex2/((ex2-er)**2+gr*ex2)
        tau(1,n,1)=xtau
        tau(2,n,1)=0.
        gam(n)=0.0d0
        do 50 j=1,n
          gam(n)=gam(n)+ro(j,1)*tau(1,n+1-j,1)
 50     continue
        do 60 l=1,4
          pav(l)=pav(l)+cn(l,n,jn,jz)
 60     continue
        pav(5)=pav(5)+ex*cn(1,n,jn,jz)
        pav(6)=pav(6)+ex2*cn(1,n,jn,jz)
        if(cn(1,n,jn,jz).ge.prec) maxn=n
 70   continue
      nmax(jn,jz)=maxn
      if(maxn.lt.1) return
      if(ip.lt.2) return
      p0=pav(1)
      do 80 l=2,6
        pav(l)=pav(l)/pav(1)
 80   continue
      pav(3)=2.*amu*amax1(0.,pav(3))
      pav(4)=sqrt(amax1(0.,pav(4))+pav(3))
      pav(3)=sqrt(pav(3))
      pav(6)=sqrt(amax1(0.,pav(6)-pav(5)**2))
      write(8,800) znc,anc,pav(1),pav(5),pav(2),pav(6),pav(3),pav(4)
      if(ip.lt.3) return
      write(8,820) nmax(jn,jz)
      write(8,850) (cn(1,n,jn,jz),n=1,nmax(jn,jz))
      write(8,830) nmax(jn,jz)
      write(8,860) (ro(n,1),n=1,nmax(jn,jz))
      n=1
      write(8,840) n,n,nmax(jn,jz)
      write(8,860) (pi2i*gam(n)/ro(n,1),n=1,nmax(jn,jz))
      return
 800  format(/,2x,'z=',f4.0,3x,'a=',f4.0,3x,'p=',f8.3,3x,'e=',f7.2,
     1 3x,'dv=',f7.4,/,32x,'se=',f7.2,3x,'sp=',f7.2,2x,'spp=',f7.2)
 820  format(/,' p(n)      nmax= ',i3)
 830  format(/,' ro(n)     nmax= ',i3)
 840  format(/,' gam(n,i)',5x,'i= ',i3,'    nmin= ',i3,'    nmax= ',i3)
 850  format(10f9.4)
 860  format(1p10d9.2)
      end
c
c-----------------------------------------------------------------------====
c
      subroutine partem(jn,jz,inc,nb,gam,ro,tau)
c
      include 'fragyld.par'
      double precision ro(nex,kem),gam(nex),gamx(nex)
      dimension tau(2,nex,kem),nb(kem)
      dimension iz(19),in(19),ba(19),bs(19),bex(19)
      common/nucs/cn(4,nex,kdn,kdz),nmax(kdn,kdz),za(2),acx(3),de,ip,p0
      data iz/0,0,1,2,1,1,2,2,3,3,4,2,3,4,5,3,4,5,6/
      data in/0,1,0,2,1,2,1,4,3,4,3,6,5,4,3,6,5,4,3/
      data ba/0.,1.,1.,4.,2.,3.,3.,6.,6.,7.,7.,8.,8.,8.,8.,9.,9.,9.,9./
      data bs/3.,2.,2.,1.,3.,2.,2.,1.,3.,4.,4.,1.,5.,1.,5.,4.,4.,4.,4./
      data bex/ 0.000, 8.071, 7.289, 2.425,13.136,14.950,14.931,17.592,
     1         14.086,14.907,15.769,31.598,20.945, 4.942,22.920,24.954,
     2         11.348,12.416,28.913/
      data pi12/37.69911/,pi2i/0.159155/,hi/33.00/
c partem calculates particle emission level densities, transmission factors,
c partial widths and related quantities.
      znc=za(1)-jz+1.
      anc=za(2)-jz-jn+2.
      bb0=excess(znc,anc,acx(2))
      ecmax=de*nmax(jn,jz)
      strtch=1.5
      do 70 i=2,inc
        znb=znc-iz(i)
        anb=anc-ba(i)
        pmu=ba(i)*anb/anc
        hnb=hi/anb**(5./3.)
        if(acx(1).le.0.) then
          aa=acn(znb,anb,pair)
         else
          aa=anb/acx(1)
          pair=2.*(1+int(znb+0.5)/2+int(anb-znb+0.5)/2)-anb
          pair=pair*acx(2)/sqrt(anb)
         endif
        ro0=aa*(hnb*aa)**1.5/pi12
        t0=100./(20.+anb)
        if(acx(3).gt.0.5) then
          egc=aa*t0
          ep=0.25*egc*egc
          egc=1./egc
          do 10 it=1,10
            sqp=sqrt(ep)
            f=1./sqp+0.5/ep-2.5/(ep+2)-egc
            df=0.5*(1./sqp+1./ep)/ep-2.5/(ep+2)**2
            dep=f/df
            if(abs(f).le.1.0e-4.and.abs(dep).lt.1.0e-4) go to 20
            ep=amax1(ep+dep,0.5*ep)
 10       continue
 20       egc=ep/aa+pair
          rogc=ro0*sqp*exp(2.*sqp-egc/t0)/(ep+2.)**2.5
         endif
        bb=excess(znb,anb,acx(2))+bex(i)-bb0
        vb0=vv(znb,anb,i-1,hl,ci,fco,tau0)
        sig0=ci/(hnb+hl)
        vb=amax1(0.,-vb0)
        vc=strtch*abs(vb0)
        ebmax=ecmax-bb+amin1(0.,vc-vb-vb0)
        nbmax=ebmax/de+1
        if(nbmax.gt.nex) go to 100
        nb(i)=nmax(jn,jz)-nbmax
        jbn=jn+in(i)
        jbz=jz+iz(i)
        jbd=min0(kdn-jbn,kdz-jbz)+1
        if(jbd.gt.0) then
          nbmax0=nmax(jbn,jbz)
          nmax(jbn,jbz)=max0(nbmax,nbmax0)
         else
          nbmax0=nbmax
         endif
        ex=-de
        eb=de*(nb(i)-1)-bb-vb0
        nbn=nb(i)
        do 50 n=1,nbmax
          if(jbd.gt.0.and.n.gt.nbmax0) then
            do 30 l=1,4
            cn(l,n,jbn,jbz)=0.
 30         continue
           endif
          ex=ex+de
          if(acx(3).gt.0.5) then
            if(ex.gt.egc) then
              ep=aa*(ex-pair)
              sqp=sqrt(ep)
              t=1./(aa*(1./sqp+0.5/ep-2.5/(ep+2.)))
              ro(n,i)=ro0*sqp*dexp(dble(2.*sqp))/(ep+2.)**2.5
             else
              ro(n,i)=rogc*exp(ex/t0)
              t=t0
             endif
           else
            ep=aa*amax1(ex,1.0e-8)
            sqp=sqrt(ep)
            t=aa*(1./sqp+0.5/ep-2.5/(ep+2.))
            t=amax1(t0,1./t)
            ro(n,i)=ro0*sqp*dexp(dble(2.*sqp-pair/t))/(ep+2.)**2.5
           endif
        if(i.eq.9) ro(n,i)=ro(n,i)*(1.+exp(-3.56/t)/3.)
        if(i.eq.10) ro(n,i)=ro(n,i)*(1.+0.5*exp(-0.478/t))
        if(i.eq.11) ro(n,i)=ro(n,i)*(1.+0.5*exp(-0.431/t))
        if(i.eq.13) ro(n,i)=ro(n,i)*(1.+1.2*exp(-0.975/t))
        if(i.eq.17) ro(n,i)=ro(n,i)*(1.+1.5*exp(-2.43/t))
        eb=eb+de
        chi=-eb/ci
        if(eb.lt.vb) then
          y=(vb-eb)/vc
          if(y.lt.1.) then
            sy=sqrt(y)
            chi=chi+fco*(asin(sy)/sqrt(1.-y)-sy*(1.+2.*y/3.))
            chi=amin1(chi,80.)
           else
            chi=80.
           endif
         endif
        sig=sig0*(alog(1.+exp(chi+1.0))-chi)
        tau(1,n,i)=bs(i)*tau0*sig/(1.+exp(chi))
        tau(2,n,i)=pmu*(eb+vb0)
        if(tau(2,n,i).lt.0.) tau(1,n,i)=0.
        nbn=nbn+1
        if(nbn.gt.0) then
          gamx(nbn)=0.0d0
          do 40 j=1,n
            gamx(nbn)=gamx(nbn)+ro(j,i)*tau(1,n+1-j,i)
 40       continue
          gam(nbn)=gam(nbn)+gamx(nbn)
         endif
 50     continue
        if(ip.lt.3.or.nb(i).gt.nmax(jn,jz)) go to 60
        write(8,800) i,nb(i)+1,nmax(jn,jz)
        write(8,820)
     1       (pi2i*gamx(n)/ro(n,1),n=max0(nb(i)+1,1),nmax(jn,jz))
 60     strtch=1.0
        if(jbd.lt.1) nb(i)=nmax(jn,jz)
 70   continue
      if(ip.lt.3) return
      write(8,810) nmax(jn,jz)
      write(8,820) (pi2i*gam(n)/ro(n,1),n=1,nmax(jn,jz))
      return
 100  write(8,830) znc,anc,i,nbmax
      stop
 800  format(/,' gam(n,i)',12x,'i= ',i3,'    nmin= ',i3,'    nmax= ',i3)
 810  format(/,' gam(n,*)',9x,'nmax= ',i3)
 820  format(1p10d9.2)
 830  format(/,' Dimension of tau and ro arrays exceeded in PARTEM for'
     1,/,'  Z=',f3.0,'  A=',f4.0,'  Particle type=',i1, '  nbmax=',i3)
      end
c
c-----------------------------------------------------------------------====
c
      function vv(z1,a1,i,hl,c1,fco,tau)
c
      dimension bz(18),ba(18),rr(8,7)
      data bz/0.,1.,2.,1.,1.,2.,2.,3.,3.,4.,2.,3.,4.,5.,3.,4.,5.,6./
      data ba/1.,1.,4.,2.,3.,3.,6.,6.,7.,7.,8.,8.,8.,8.,9.,9.,9.,9./
      data rr/
     *3.99470,1.28115,1.54345,0.86526,0.84958,-0.07239,1.33996, 0.00918,
     *4.34979,1.26223,3.85508,1.34366,0.46392, 0.10148,1.33686, 0.04788,
     *3.78748,1.25858,3.98672,1.29554,0.25695, 0.10731,1.04415, 0.01944,
     *5.10514,1.23542,4.80914,1.20796, 0.44833,0.07287,1.24211, 0.01813,
     *5.52216,1.14864,6.80373,0.90997, 0.11559,0.11715,1.22437,-0.00776,
     *5.28751,1.31716,5.72916,1.01125,-0.24798,0.29497,1.33270,-0.03607,
     *1.23300,-0.9780,2.72000,1.07000, 1.00000,0.00000,1.00000, 0.00000/
      data hm/20.734/
c vv estimates the barrier height and curvature for particle emission channels
      at=a1**(1./3.)
      fco=hm*(a1+ba(i))/(a1*ba(i))
      zz=z1*bz(i)
      if(zz.lt.1.0e-3) zz=-1.
      if(i.lt.7) then
        r1=rr(1,i)+at*rr(2,i)
        r2=rr(3,i)+at*rr(4,i)
        hl=fco/(r1*r1)
        vv=1.44*zz/r2
        c1=rr(5,i)+at*rr(6,i)
        tau=1./(rr(7,i)+at*rr(8,i))
       else
        r1=rr(1,7)*at+rr(2,7)/at
        ab=ba(i)**(1./3.)
        r2=rr(1,7)*ab+rr(2,7)/ab
        vv=50.*r1*r2/(r1+r2)
        r1=r1+r2
        r2=rr(3,7)+(at+ab)*rr(4,7)
        vv=1.44*zz/r2-vv*exp((r1-r2)/0.63)
        hl=fco/(r2*r2)
        c1=rr(5,7)
        tau=rr(7,7)
       endif
      fco=sqrt(3.*abs(vv)/fco)*r2
      return
      end
c
c-----------------------------------------------------------------------====
c
      function erf(x)
c
      data sqpi/0.564190/
c error function
      x2=x*x
      if(x2.lt.12.5) go to 10
      erf=sign(1.,x)
      return
 10   ex=exp(-x2)
      x2=2.*x2
      erf=1.
      drf=1.
      xd=1.
      do 20 i=1,30
        xd=xd+2.
        drf=x2*drf/xd
        erf=erf+drf
        if(drf/erf.lt.1.0e-6) go to 30
 20   continue
 30   erf=2.*sqpi*ex*x*erf
      return
      end
c
c-----------------------------------------------------------------------====
c
      function algam(x)
      data alpi/0.918938533/
c natural log of the gamma function
      z=x
      if(x.gt.10.) go to 10
      nz=10-int(z)
      z=z+nz
 10   zz=z*z
      algam=(1.-(1./6.-(1./3.-0.25/zz)/(7.*zz))/(5.*zz))/(12.*z)
      algam=algam+(z-0.5)*alog(z)-z+alpi
      if(x.gt.10.) return
      do 20 n=1,nz
        z=z-1
      algam=algam-alog(z)
 20   continue
      return
      end
c
c-----------------------------------------------------------------------====
c
      subroutine chsym(iz,iaa,simaz)
c
      character amt(112)*2,nn(11)*1,simaz*5
      data amt/
     *'H ','He','Li','Be','B ','C ','N ','O ','F ','Ne','Na','Mg','Al',
     *'Si','P ','S ','Cl','Ar','K ','Ca','Sc','Ti','V ','Cr','Mn','Fe',
     *'Co','Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr','Rb','Sr','Y ',
     *'Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb','Te',
     *'I ','Xe','Cs','Ba','La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb',
     *'Dy','Ho','Er','Tm','Yb','Lu','Hf','Ta','Wr','Re','Os','Ir','Pt',
     *'Au','Hg','Tl','Pb','Bi','Po','At','Rn','Fr','Ra','Ac','Th','Pa',
     *'U ','Np','Pu','Am','Cm','Bk','Cf','Es','Fm','Md','No','Lr','Rf',
     *'Db','Sg','Bh','Hs','Mt','10','11','12'/
      data nn/'0','1','2','3','4','5','6','7','8','9',' '/
c chsym returns the chemical symbol as a 5 character string in the form nnnaa
      ia=iaa
      do 10 i=1,2
        na=ia
        ia=ia/10
        na=na-ia*10
        if(na.eq.0.and.ia.eq.0) na=10
        simaz(4-i:4-i)=nn(na+1)
 10   continue
      if(ia.eq.0) ia=10
      simaz(1:1)=nn(ia+1)
      simaz(4:5)=amt(iz)
      return
      end
c
c-----------------------------------------------------------------------====
c
      function acn(zc,ac,pair)
c
      dimension defcon(2),sz(100),sn(150),pz(100),pn(150)
      data  pz/11*0.,2.46,0.,2.09,0.,1.62,0.,1.62,0.,1.83,0.,1.73,
     1 0.,1.35,0.,1.54,0.,1.28,0.26,0.88,0.19,1.35,-.05,1.52,-.09,1.17,
     2 .04,1.24,0.29,1.09,.26,1.17,.23,1.15,-.08,1.35,0.34,1.05,.28,1.27
     3 ,0.,1.05,0.,1.,.09,1.2,.2,1.4,.93,1.,-.2,1.19,.09,.97,0.,.92,.11,
     4 .68,.05,.68,-.22,.79,.09,.69,.01,.72,0.,.4,.16,.73,0.,.46,.17,
     5 .89,0.,.79,0.,.89,0.,.81,-.06,.69,-.2,.71,-.12,.72,0.,.77,2*0./
      data pn/11*0.,2.67,0.,1.8,0.,1.67,0.,1.86,0.,2.04,0.,1.64,0.,1.44,
     1 0.,1.54,0.,1.3,0.,1.27,0.,1.29,.08,1.41,-.08,1.5,-.05,2.24,-.47,
     2 1.43,-.15,1.44,.06,1.56,.25,1.57,-.16,1.46,0.,.93,.01,.62,-.5,
     3 1.42,.13,1.52,-.65,.8,-.08,1.29,-.47,1.25,-.44,.97,.08,1.65,-.11,
     4 1.26,-.46,1.06,0.22,1.55,-.07,1.37,0.1,1.2,-.27,.92,-.35,1.19,0.,
     5 1.05,-.25,1.61,-.21,.9,-.21,.74,-.38,.72,-.34,.92,-.26,.94,.01,
     6 .65,-.36,.83,.11,.67,.05,1.,.51,1.04,.33,.68,-.27,.81,.09,.75,
     7 .17,.86,.14,1.1,-.22,.84,-.47,.48,.02,.88,.24,.52,.27,.41,-.05,
     8  .38,.15,.67,0.,.61,0.,.78,0.,.67,0.,.67,0.,.79,
     9 0.,.6,.04,.64,-.06,.45,.05,.26,-.22,.39,0.0,.39/
      data sz/10*0.,-2.91,-4.17,-5.72,-7.8,-8.97,-9.7,-10.1,-10.7,-11.38
     1 ,-12.07,-12.55,-13.24,-13.93,-14.71,-15.53,-16.37,-17.36,-18.6,
     2 -18.7,-18.01,-17.87,-17.08,-16.6,-16.75,-16.5,-16.35,-16.22,
     3 -16.41,-16.89,-16.43,-16.68,-16.73,-17.45,-17.29,-17.44,-17.82,
     4 -18.62,-18.27,-19.39,-19.91,-19.14,-18.26,-17.4,-16.42,-15.77,
     5 -14.37,-13.91,-13.1,-13.11,-11.43,-10.89,-10.75,-10.62,-10.41,
     6 -10.21,-9.85,-9.47,-9.03,-8.61,-8.13,-7.46,-7.48,-7.2,-7.13,-7.06
     7 ,-6.78,-6.64,-6.64,-7.68,-7.89,-8.41,-8.49,-7.88,-6.3,-5.47,-4.78
     8 ,-4.37,-4.17,-4.13,-4.32,-4.55,-5.04,-5.28,-6.06,-6.28,-6.87,
     9  -7.20,-7.74,2*0./
      data sn/10*0.,6.8,7.53,7.55,7.21,7.44,8.07,8.94,9.81,10.6,11.39,
     1 12.54,13.68,14.34,14.19,13.83,13.5,13.,12.13,12.6,13.26,14.13,
     2 14.92,15.52,16.38,17.16,17.55,18.03,17.59,19.03,18.71,18.8,18.99,
     3 18.46,18.25,17.76,17.38,16.72,15.62,14.38,12.88,13.23,13.81,14.9,
     4 14.86,15.76,16.2,17.62,17.73,18.16,18.67,19.69,19.51,20.17,19.48,
     5 19.98,19.83,20.2,19.72,19.87,19.24,18.44,17.61,17.1,16.16,15.9,
     6 15.33,14.76,13.54,12.63,10.65,10.1,8.89,10.25,9.79,11.39,11.72,
     7 12.43,12.96,13.43,13.37,12.96,12.11,11.92,11.,10.8,10.42,10.39,
     8 9.69,9.27,8.93,8.57,8.02,7.59,7.33,7.23,7.05,7.42,6.75,6.6,6.38,
     9  6.36,6.49,6.25,5.85,5.48,4.53,4.3,3.39,2.35,1.66,.81,
     a 0.46,-.96,-1.69,-2.53,-3.16,-1.87,-.41,.71,1.66,2.62,3.22,3.76,
     b 4.1,4.46,4.83,5.09,5.18,5.17,5.1,5.01,4.97,5.09,5.03,4.93,5.28,
     c 5.49,5.50,5.37,5.30/
      data defcon/0.142,0.120/
c Gilbert-Cameron tables of a and delta
      ia=ac
      iz=zc
      in=ia-iz
      idefcn=idef(iz,in)+1.01
      acn= ac*(0.00917*(sz(iz)+sn(in)) + defcon(idefcn))
      pair = pz(iz)+pn(in)
      return
      end
c
c-----------------------------------------------------------------------====
c
      function idef(iz,in)
c
      dimension m(6),i(2)
      data m/54,78,86,122,130,180/
c idef determines if a nucleus lies in a region of deformed nuclei
      i(1)=iz
      i(2)=in
      idef=0
      do 20 j=1,2
        do 10 k=1,5,2
          if(i(j).ge.m(k).and.i(j).le.m(k+1)) idef=1
          if(idef.eq.1) return
 10     continue
 20   continue
      return
      end
c
c-----------------------------------------------------------------------====
c
      subroutine geoyld(zp,ap,at,yldg)
c
      include 'fragyld.par'
      dimension yldg(3,kdn,kdz)
      data pi/3.14159/,r0/1.2/,es0/0.95/
      kya=min0(kdz+kdn,int(ap))
      nz=min0(int(zp)+1,kdz)
      np=ap-zp
      nn=min0(np+1,kdn)
      algp=algam(zp+1.)+algam(np+1.)-algam(ap+1.)
      rp=r0*ap**(1./3.)
      rt=r0*at**(1./3.)
      amin=0.
      if(rp.gt.rt) amin=ap*abvol(b,rp,rt)
      kamax=kya+1
      ab=ap+0.5
      b=rp+rt
      do 50 ka=1,kya
        ab=ab-1.
        b1=b
        bh=b
        bl=b
        if(ab.gt.amin) go to 10
        if(ka.gt.kamax) go to 30
        b=0.
        kamax=ka
        go to 30
 10     bl=amax1(bl-0.25,0.)
        da=ap*abvol(bl,rp,rt)-ab
        if(da.gt.0.) go to 10
        do 20 i=1,20
          b=0.5*(bh+bl)
          da=ap*abvol(b,rp,rt)-ab
          if(abs(da).lt.1.0e-6) go to 30
          if(da.lt.0.) then
            bl=b
           else
            bh=b
           endif
 20     continue
 30     if(ka.eq.1.or.ka.gt.kamax) then
          algb=algp+algam(ap+1.)
          yld=10.*pi*b*b
          bs=0.
          esf=0.
         else
          algb=algp+algam(ab+1.5)+algam(ap-ab+0.5)
          yld=10.*pi*(b1*b1-b*b)
          bs=0.5*(b+b1)
          rb=r0*(ab+0.5)**(1./3.)
          esf=es0*(absur(bs,rp,rt)-4.*pi*rb*rb)
         endif
        knmi=max0(ka-nz,0)+1
        knma=min0(nn,ka)
        kz=ka+2-knmi
        do 40 kn=knmi,knma
          kz=kz-1
          zayld=yld*exp(algb-algam(float(kn))-algam(float(kz))
     1              -algam(float(np-kn+2))-algam(zp-kz+2.))
          yldg(1,kn,kz)=zayld
          yldg(2,kn,kz)=bs
          yldg(3,kn,kz)=esf
 40     continue
 50   continue
      if(int(zp).lt.kdz.and.np.lt.kdn) then
        yldg(1,nn,nz)=10.*pi**b*b
        yldg(2,nn,nz)=0.
        yldg(3,nn,nz)=0
       endif
      return
      end
c
c-----------------------------------------------------------------------====
c
      function abvol(b0,rp,rt)
c
      data pi/3.14159/,dr0/0.1/,eps/1.0e-2/
      b=amax1(b0,eps)
      b2=b*b
      bi=1./b
      rp2=rp*rp
      rt2=rt*rt
      rl=abs(b-rt)
      rh=amin1(rp,b+rt)
      nl=(rh-rl)/dr0
      nl=2*((nl/2)+1)
      dr=(rh-rl)/nl
      wt=4.*dr/3.
      if(b.lt.rt) then
        sqr=sqrt(rp2-rl*rl)
        abvol=4.*pi*sqr**3/3.-wt*pi*rl*sqr
       else
        abvol=4.*pi*rp**3/3.
       endif
      if(dr.gt.1.0e-4) then
        dwt=-2.*wt
        wt=4.*wt
        do 10 n=2,nl
          rl=rl+dr
          rl2=rl*rl
          phi=0.5*bi*(rl2-rt2+b2)/rl
          phi=acos(phi)
          abvol=abvol-wt*phi*rl*sqrt(rp2-rl2)
          wt=wt+dwt
          dwt=-dwt
 10     continue
       endif
      abvol=3.*abvol/(4.*pi*rp**3)
      return
      end
c
c-----------------------------------------------------------------------====
c
      function absur(b0,rp,rt)
c
      data pi/3.14159/,dr0/0.1/,eps/1.0e-2/
      b=amax1(b0,eps)
      b2=b*b
      bi=1./b
      rp2=rp*rp
      rt2=rt*rt
      rl=abs(b-rt)
      rh=amin1(rp,b+rt)
      nl=max0(int((rh-rl)/dr0)+1,5)
      dr=(rh-rl)/nl
      nl=nl-3
      if(b.lt.rt) then
        sqr=sqrt(rp2-rl*rl)
        absur=4.*pi*rp*sqr
        if(dr.le.1.0e-4) absur=absur+4.*pi*rt*sqr
       else
        absur=4.*pi*rp2
       endif
      if(dr.gt.1.0e-4) then
        rl=rl+dr
        rl2=rl*rl
        sqr=sqrt(rp2-rl2)
        phi=0.5*bi*(rl2-rt2+b2)/rl
        phi=acos(phi)
        absur=absur+6.*dr*(rt*bi*sqr/sin(phi)-rl*rp*phi/sqr)
        wt=4.*dr
        if(nl.gt.1) then
          do 10 n=1,nl
            rl=rl+dr
            rl2=rl*rl
            sqr=sqrt(rp2-rl2)
            phi=0.5*bi*(rl2-rt2+b2)/rl
            phi=acos(phi)
            absur=absur+wt*(rt*bi*sqr/sin(phi)-rl*rp*phi/sqr)
 10       continue
         endif
        rl=rl+dr
        rl2=rl*rl
        sqr=sqrt(rp2-rl2)
        phi=0.5*bi*(rl2-rt2+b2)/rl
        phi=acos(phi)
        absur=absur+6.*dr*rt*bi*sqr/sin(phi)
        if(rh.lt.rp) then
          absur=absur-6.*dr*rl*rp*phi/sqr
         else
          absur=absur-2.*rp*phi*(dr*rl/sqr+2.*sqr)
         endif
       endif
      return
      end
c
c-----------------------------------------------------------------------====
c
      function second()
c
c  Returns seconds since midnight. The LINUX g77 compiler has a routine with
c  the same name and function. When using it, comment out the entire routine.
c
c - IBM PC Microsoft FORTRAN
      call gettim (ihr, imin, isec, i100)
      second = real ((ihr*3600 + imin*60 + isec)*100 + i100)/100
c - IBM PC Lahey FORTRAN
c      call timer(i100)
c      second=real(i100)/100
c - VAX and CDC
c      second=secnds(0.)
      return
      end
c
c-----------------------------------------------------------------------====
c
      function excess(zz,aa,apx)
c
c   Audi-Wapstra table
      parameter(numz=112,numiso=2931)
c
c   Moller-Nix table
c      parameter(numz=137,numiso=8983)
c
      dimension z(3,3),zm(3),c(9)
c
      common/xcess/iz(numz),ia(numz),p(numiso)
c
c  Audi & Wapstra -- 1998 -- all Z
      data c/-2.308,-7.263,30.496,19.344,52.151,0.311,0.720,
     1        2.034,11.236/
c
c  Audi & Wapstra -- 1998 -- Z>=8
c      data c/-2.094,-7.017,28.636,17.161,55.235,4.037,0.709,
c     1        1.929,11.505/
c
c  Moller & Nix -- 1998
c      data c/-1.652,-6.757,30.381,14.242,69.727,9.487,0.690,
c     1        2.295,12.990/
c
c  Excess returns the mass excess from Audi & Wapstra (Moller & Nix), if
c  possible. If the nucleus is not in Audi & Wapstra (Moller & Nix), it
c  tries to extrapolate from the table. If this too fails, the routine 
c  uses a liquid drop fit to the table. Note that the liquid drop fit is
c  always used when the liquid drop pairing parameter apx is less than
c  or equal to zero.
      apr=apx
      if(apx.lt.0.) apr=c(9)
      if(apx.le.0.) go to 110
      izeta=zz+0.5
      iaa=aa+0.5
      excess=0.
      ierr=0
      if(izeta.lt.8.or.izeta.ge.numz) go to 110
      ind1=iaa-ia(izeta)
      index=iz(izeta)+ind1
      if(index.gt.numiso) go to 110
      if(index.ge.iz(izeta+1).or.ind1.lt.0) go to 20
      excess=p(index)
      return
 20   ized=izeta
      izd=2
      if(ind1.lt.0) izd=-2
      z(1,1)=0.
      z(2,1)=0.
      z(3,1)=0.
      z(3,2)=0.
      z(3,3)=0.
      zm(1)=0.
      zm(2)=0.
      zm(3)=0.
 30   ized=ized+izd
      if(ized.lt.8.or.ized.ge.numz) go to 50
      ind1=iaa-ia(ized)
      index=iz(ized)+ind1
      if(index.lt.iz(ized+1).and.ind1.ge.0) go to 40
      if(z(1,1).gt.0.or.index.gt.numiso) go to 50
      go to 30
 40   z(1,1)=z(1,1)+1.
      z(2,1)=z(2,1)+ized
      z(3,1)=z(3,1)+ized**2
      z(3,2)=z(3,2)+ized**3
      z(3,3)=z(3,3)+ized**4
      ex=p(index)
      zm(1)=zm(1)+ex
      zm(2)=zm(2)+ex*ized
      zm(3)=zm(3)+ex*ized**2
      go to 30
 50   if(z(1,1).lt.3) go to 110
      z(1,2)=z(2,1)
      z(2,2)=z(3,1)
      z(1,3)=z(3,1)
      z(2,3)=z(3,2)
      do 80 i=1,2
      if(abs(z(i,i)).lt.1.0e-35) go to 110
      i1=i+1
      do 70 j=i1,3
      fac=z(j,i)/z(i,i)
      zm(j)=zm(j)-fac*zm(i)
      do 60 k=i1,3
      z(j,k)=z(j,k)-fac*z(i,k)
 60   continue
 70   continue
 80   continue
      if(abs(z(3,3)).lt.1.0e-35) go to 110
      zm(3)=zm(3)/z(3,3)
      excess=zm(3)
      do 100 i=2,1,-1
      i1=i+1
      do 90 j=i1,3
      zm(i)=zm(i)-z(i,j)*zm(j)
 90   continue
      zm(i)=zm(i)/z(i,i)
      excess=zm(i)+izeta*excess
 100  continue
      ierr=1
      return
 110  a3=aa**(-1./3.)
      d=(1.-2.*zz/aa)/(1.+c(8)*a3)
      d=d*d
      excess=zz*c(1)+aa*(c(2)+c(3)*d+a3*(c(4)+c(5)*d+c(6)*a3))
      excess=excess+c(7)*zz*(zz-1.)*a3
      d=2.*(1+int(zz+.5)/2+int(aa-zz+.5)/2)-aa
      excess=excess-apr*d/sqrt(aa)
      ierr=-1
      return
      end
c
c-----------------------------------------------------------------------====
c
      block data
      parameter(numz=112,numiso=2931)
c
c Wapstra and Audi - 1998
      common/xcess/iz(numz),ia(numz),p(numiso)
c
      data p(   1),p(   2),p(   3)/    8.07132,    7.28897,   13.13572/
      data p(   4),p(   5),p(   6)/   14.94979,   25.92778,   36.83398/
      data p(   7),p(   8),p(   9)/   41.86376,   14.93120,    2.42491/
      data p(  10),p(  11),p(  12)/   11.38624,   17.59412,   26.11027/
      data p(  13),p(  14),p(  15)/   31.59798,   40.81837,   48.81001/
      data p(  16),p(  17),p(  18)/   25.32018,   11.67888,   14.08631/
      data p(  19),p(  20),p(  21)/   14.90767,   20.94620,   24.95390/
      data p(  22),p(  23),p(  24)/   33.05023,   40.79586,   50.09600/
      data p(  25),p(  26),p(  27)/   37.99600,   18.37447,   15.76949/
      data p(  28),p(  29),p(  30)/    4.94166,   11.34758,   12.60658/
      data p(  31),p(  32),p(  33)/   20.17397,   25.07640,   33.65845/
      data p(  34),p(  35),p(  36)/   39.88240,   27.86786,   22.92100/
      data p(  37),p(  38),p(  39)/   12.41570,   12.05076,    8.66798/
      data p(  40),p(  41),p(  42)/   13.36890,   16.56221,   23.66373/
      data p(  43),p(  44),p(  45)/   28.96694,   37.08169,   43.71631/
      data p(  46),p(  47),p(  48)/   52.32200,   59.36400,   35.09406/
      data p(  49),p(  50),p(  51)/   28.91365,   15.69857,   10.65053/
      data p(  52),p(  53),p(  54)/     .00000,    3.12501,    3.01989/
      data p(  55),p(  56),p(  57)/    9.87314,   13.69412,   21.03659/
      data p(  58),p(  59),p(  60)/   24.92404,   32.83339,   37.56006/
      data p(  61),p(  62),p(  63)/   45.96000,   52.58300,   39.69900/
      data p(  64),p(  65),p(  66)/   24.96052,   17.33808,    5.34546/
      data p(  67),p(  68),p(  69)/    2.86342,     .10144,    5.68343/
      data p(  70),p(  71),p(  72)/    7.87082,   13.11714,   15.86045/
      data p(  73),p(  74),p(  75)/   21.76649,   25.23191,   32.08089/
      data p(  76),p(  77),p(  78)/   37.73500,   47.04000,   32.04784/
      data p(  79),p(  80),p(  81)/   23.11073,    8.00646,    2.85539/
      data p(  82),p(  83),p(  84)/   -4.73700,    -.80900,    -.78206/
      data p(  85),p(  86),p(  87)/    3.33356,    3.79691,    8.06174/
      data p(  88),p(  89),p(  90)/    9.28435,   14.61637,   18.97446/
      data p(  91),p(  92),p(  93)/   27.14400,   35.16400,   33.60800/
      data p(  94),p(  95),p(  96)/   16.77700,   10.68026,    1.95170/
      data p(  97),p(  98),p(  99)/     .87343,   -1.48741,    -.01740/
      data p( 100),p( 101),p( 102)/    -.04758,    2.79378,    3.32952/
      data p( 103),p( 104),p( 105)/    7.54451,   11.26638,   18.28816/
      data p( 106),p( 107),p( 108)/   25.05003,   33.22600,   40.29600/
      data p( 109),p( 110),p( 111)/   23.99240,   16.48517,    5.30678/
      data p( 112),p( 113),p( 114)/    1.75106,   -7.04193,   -5.73172/
      data p( 115),p( 116),p( 117)/   -8.02434,   -5.15364,   -5.94752/
      data p( 118),p( 119),p( 120)/   -2.05870,     .42988,    7.09351/
      data p( 121),p( 122),p( 123)/   11.27859,   18.02059,   22.23662/
      data p( 124),p( 125),p( 126)/   30.84200,   37.17600,   25.31800/
      data p( 127),p( 128),p( 129)/   12.92862,    6.84486,   -2.18426/
      data p( 130),p( 131),p( 132)/   -5.18210,   -9.52949,   -8.41760/
      data p( 133),p( 134),p( 135)/   -9.35746,   -6.90247,   -5.58086/
      data p( 136),p( 137),p( 138)/   -1.03358,    2.61871,    8.59442/
      data p( 139),p( 140),p( 141)/   12.66376,   18.30366,   25.50989/
      data p( 142),p( 143),p( 144)/   32.50900,   41.15300,   17.57053/
      data p( 145),p( 146),p( 147)/   10.91168,    -.39677,   -5.47267/
      data p( 148),p( 149),p( 150)/  -13.93338,  -13.19273,  -16.21448/
      data p( 151),p( 152),p( 153)/  -14.58650,  -15.01875,  -10.66119/
      data p( 154),p( 155),p( 156)/   -8.88223,   -3.21509,    -.79560/
      data p( 157),p( 158),p( 159)/    5.20423,    8.45092,   16.29200/
      data p( 160),p( 161),p( 162)/   20.91200,   29.10000,   26.11900/
      data p( 163),p( 164),p( 165)/   18.18300,    6.76721,    -.05504/
      data p( 166),p( 167),p( 168)/   -8.91574,  -12.21034,  -17.19683/
      data p( 169),p( 170),p( 171)/  -16.85055,  -18.21551,  -15.87237/
      data p( 172),p( 173),p( 174)/  -14.95418,  -11.06207,   -8.50492/
      data p( 175),p( 176),p( 177)/   -2.86224,    -.05808,    5.91638/
      data p( 178),p( 179),p( 180)/    9.60370,   15.74200,   20.40000/
      data p( 181),p( 182),p( 183)/   32.16400,   23.77200,   10.75476/
      data p( 184),p( 185),p( 186)/    3.82531,   -7.14462,  -12.38443/
      data p( 187),p( 188),p( 189)/  -21.49279,  -21.89503,  -24.43288/
      data p( 190),p( 191),p( 192)/  -22.94896,  -24.08086,  -20.49238/
      data p( 193),p( 194),p( 195)/  -19.95656,  -14.35976,  -12.40064/
      data p( 196),p( 197),p( 198)/   -6.52419,   -3.74461,    2.14200/
      data p( 199),p( 200),p( 201)/    5.40300,   11.83000,   14.99700/
      data p( 202),p( 203),p( 204)/   31.99700,   18.87200,   10.97300/
      data p( 205),p( 206),p( 207)/    -.75298,   -7.16102,  -16.95191/
      data p( 208),p( 209),p( 210)/  -20.20056,  -24.44099,  -24.30532/
      data p( 211),p( 212),p( 213)/  -26.33773,  -24.55755,  -24.85761/
      data p( 214),p( 215),p( 216)/  -20.25084,  -18.99471,  -14.46610/
      data p( 217),p( 218),p( 219)/  -12.64969,   -8.33687,   -4.84377/
      data p( 220),p( 221),p( 222)/     .08400,    3.08300,    9.20300/
      data p( 223),p( 224),p( 225)/   14.10300,   22.19700,   25.97000/
      data p( 226),p( 227),p( 228)/   17.50700,    4.07311,   -3.15888/
      data p( 229),p( 230),p( 231)/  -14.06281,  -19.04493,  -26.01598/
      data p( 232),p( 233),p( 234)/  -26.58624,  -29.93185,  -28.84637/
      data p( 235),p( 236),p( 237)/  -30.66396,  -26.89622,  -26.86108/
      data p( 238),p( 239),p( 240)/  -23.16134,  -22.84955,  -18.60193/
      data p( 241),p( 242),p( 243)/  -17.24195,  -12.48202,  -10.88000/
      data p( 244),p( 245),p( 246)/   -4.82500,    -.40100,    7.09800/
      data p( 247),p( 248),p( 249)/   12.10000,   20.50200,   26.55700/
      data p( 250),p( 251),p( 252)/   13.14300,    4.44300,   -7.06444/
      data p( 253),p( 254),p( 255)/  -13.33070,  -21.00351,  -24.44057/
      data p( 256),p( 257),p( 258)/  -29.01351,  -29.52189,  -31.76152/
      data p( 259),p( 260),p( 261)/  -29.79798,  -29.80065,  -27.55773/
      data p( 262),p( 263),p( 264)/  -27.33915,  -24.98733,  -24.02939/
      data p( 265),p( 266),p( 267)/  -19.99106,  -18.90932,  -14.79200/
      data p( 268),p( 269),p( 270)/  -11.22500,   -4.79700,    -.10200/
      data p( 271),p( 272),p( 273)/    7.20000,   12.60300,   20.08300/
      data p( 274),p( 275),p( 276)/   11.29600,   -2.17908,   -9.38134/
      data p( 277),p( 278),p( 279)/  -18.37827,  -23.04821,  -30.23044/
      data p( 280),p( 281),p( 282)/  -30.94803,  -34.71476,  -33.24184/
      data p( 283),p( 284),p( 285)/  -35.03989,  -33.06726,  -34.42207/
      data p( 286),p( 287),p( 288)/  -31.97753,  -32.26204,  -29.71933/
      data p( 289),p( 290),p( 291)/  -29.72074,  -25.90835,  -23.22200/
      data p( 292),p( 293),p( 294)/  -16.59900,  -13.09700,   -6.29700/
      data p( 295),p( 296),p( 297)/   -1.70500,    5.80000,   20.41800/
      data p( 298),p( 299),p( 300)/    6.76300,   -1.48100,  -11.16711/
      data p( 301),p( 302),p( 303)/  -17.42508,  -24.79924,  -28.80169/
      data p( 304),p( 305),p( 306)/  -33.80684,  -33.53502,  -35.55887/
      data p( 307),p( 308),p( 309)/  -35.02132,  -36.59304,  -35.81021/
      data p( 310),p( 311),p( 312)/  -36.60803,  -35.41893,  -35.69689/
      data p( 313),p( 314),p( 315)/  -32.12448,  -30.32005,  -25.35263/
      data p( 316),p( 317),p( 318)/  -22.00200,  -16.19900,  -11.99800/
      data p( 319),p( 320),p( 321)/   -5.59800,    -.57000,   13.15300/
      data p( 322),p( 323),p( 324)/    4.43900,   -6.43920,  -13.16061/
      data p( 325),p( 326),p( 327)/  -22.05905,  -27.27626,  -34.84611/
      data p( 328),p( 329),p( 330)/  -35.13749,  -38.54677,  -38.40844/
      data p( 331),p( 332),p( 333)/  -41.46909,  -40.81253,  -43.13491/
      data p( 334),p( 335),p( 336)/  -42.33970,  -44.21474,  -41.29005/
      data p( 337),p( 338),p( 339)/  -39.57146,  -35.88651,  -32.50914/
      data p( 340),p( 341),p( 342)/  -27.89800,  -23.58500,  -18.11800/
      data p( 343),p( 344),p( 345)/  -13.23700,   -7.12000,   13.89800/
      data p( 346),p( 347),p( 348)/    2.84100,   -4.93700,  -14.16801/
      data p( 349),p( 350),p( 351)/  -20.52639,  -28.64221,  -32.12093/
      data p( 352),p( 353),p( 354)/  -36.18762,  -37.81581,  -41.06934/
      data p( 355),p( 356),p( 357)/  -41.75864,  -44.33163,  -44.49281/
      data p( 358),p( 359),p( 360)/  -46.55228,  -44.53751,  -43.21880/
      data p( 361),p( 362),p( 363)/  -40.38026,  -37.96800,  -34.46527/
      data p( 364),p( 365),p( 366)/  -30.33900,  -25.46700,  -21.38700/
      data p( 367),p( 368),p( 369)/  -15.77000,  -11.14000,    9.10100/
      data p( 370),p( 371),p( 372)/    1.23200,   -8.85021,  -15.71300/
      data p( 373),p( 374),p( 375)/  -25.12088,  -29.32032,  -37.54830/
      data p( 376),p( 377),p( 378)/  -39.00691,  -44.12534,  -44.93173/
      data p( 379),p( 380),p( 381)/  -48.48700,  -48.55804,  -51.42585/
      data p( 382),p( 383),p( 384)/  -49.72685,  -49.46403,  -46.82461/
      data p( 385),p( 386),p( 387)/  -45.76429,  -41.80545,  -39.13206/
      data p( 388),p( 389),p( 390)/  -34.55800,  -31.56800,  -26.11900/
      data p( 391),p( 392),p( 393)/  -22.69100,  -16.75000,   10.33000/
      data p( 394),p( 395),p( 396)/    -.24200,   -8.16900,  -18.02400/
      data p( 397),p( 398),p( 399)/  -23.84600,  -31.87359,  -37.07393/
      data p( 400),p( 401),p( 402)/  -42.00393,  -44.47466,  -47.95618/
      data p( 403),p( 404),p( 405)/  -49.21754,  -52.19749,  -51.43741/
      data p( 406),p( 407),p( 408)/  -51.84461,  -49.88673,  -49.14730/
      data p( 409),p( 410),p( 411)/  -46.23936,  -44.37637,  -40.38026/
      data p( 412),p( 413),p( 414)/  -37.91180,  -33.06803,  -30.35700/
      data p( 415),p( 416),p( 417)/  -25.02000,  -21.65700,    5.99000/
      data p( 418),p( 419),p( 420)/   -2.13600,  -13.53500,  -19.41200/
      data p( 421),p( 422),p( 423)/  -29.47093,  -34.55236,  -42.81532/
      data p( 424),p( 425),p( 426)/  -45.32544,  -50.25446,  -51.44476/
      data p( 427),p( 428),p( 429)/  -55.41280,  -55.28064,  -56.92832/
      data p( 430),p( 431),p( 432)/  -55.10330,  -55.28860,  -52.39299/
      data p( 433),p( 434),p( 435)/  -51.93078,  -47.85084,  -46.82620/
      data p( 436),p( 437),p( 438)/  -42.76488,  -41.17203,  -35.52700/
      data p( 439),p( 440),p( 441)/  -33.34700,  -27.60000,    6.39900/
      data p( 442),p( 443),p( 444)/   -5.11400,  -12.37000,  -22.26300/
      data p( 445),p( 446),p( 447)/  -28.99700,  -37.61054,  -42.62147/
      data p( 448),p( 449),p( 450)/  -48.23696,  -50.70114,  -54.68363/
      data p( 451),p( 452),p( 453)/  -55.55127,  -57.70639,  -56.90555/
      data p( 454),p( 455),p( 456)/  -57.48486,  -55.90226,  -55.47310/
      data p( 457),p( 458),p( 459)/  -52.91444,  -51.73517,  -48.46563/
      data p( 460),p( 461),p( 462)/  -46.75168,  -43.10022,  -40.89258/
      data p( 463),p( 464),p( 465)/  -36.49600,  -33.70100,   13.56300/
      data p( 466),p( 467),p( 468)/     .75500,   -6.62300,  -18.10800/
      data p( 469),p( 470),p( 471)/  -24.58200,  -34.47150,  -40.21731/
      data p( 472),p( 473),p( 474)/  -48.32914,  -50.94128,  -56.24841/
      data p( 475),p( 476),p( 477)/  -57.47501,  -60.60101,  -60.17571/
      data p( 478),p( 479),p( 480)/  -62.14885,  -60.65842,  -61.40693/
      data p( 481),p( 482),p( 483)/  -58.91749,  -58.89790,  -55.77931/
      data p( 484),p( 485),p( 486)/  -55.07923,  -51.28805,  -50.31930/
      data p( 487),p( 488),p( 489)/  -46.57469,  -44.23700,  -39.40200/
      data p( 490),p( 491),p( 492)/    1.63900,   -9.57600,  -17.19500/
      data p( 493),p( 494),p( 495)/  -27.27400,  -33.91600,  -42.63914/
      data p( 496),p( 497),p( 498)/  -48.00533,  -54.02371,  -56.03501/
      data p( 499),p( 500),p( 501)/  -59.33967,  -59.84143,  -62.22361/
      data p( 502),p( 503),p( 504)/  -61.64422,  -62.89505,  -61.42810/
      data p( 505),p( 506),p( 507)/  -61.83702,  -59.78931,  -59.16423/
      data p( 508),p( 509),p( 510)/  -56.05226,  -55.32143,  -51.82832/
      data p( 511),p( 512),p( 513)/  -51.04586,  -46.75200,  -44.96300/
      data p( 514),p( 515),p( 516)/  -40.60400,   -3.79100,  -11.43900/
      data p( 517),p( 518),p( 519)/  -22.65400,  -29.37900,  -39.20610/
      data p( 520),p( 521),p( 522)/  -45.32991,  -53.89965,  -56.07548/
      data p( 523),p( 524),p( 525)/  -60.22302,  -61.15113,  -64.46810/
      data p( 526),p( 527),p( 528)/  -64.21677,  -66.74269,  -65.50922/
      data p( 529),p( 530),p( 531)/  -67.09590,  -65.12259,  -66.02873/
      data p( 532),p( 533),p( 534)/  -63.74246,  -63.48603,  -60.37772/
      data p( 535),p( 536),p( 537)/  -59.48520,  -55.88964,  -54.67870/
      data p( 538),p( 539),p( 540)/  -50.22600,  -48.52200,  -43.80800/
      data p( 541),p( 542),p( 543)/  -41.61000,  -36.48700,  -33.72000/
      data p( 544),p( 545),p( 546)/   -2.62700,  -13.46000,  -21.69400/
      data p( 547),p( 548),p( 549)/  -31.62400,  -38.60100,  -47.30527/
      data p( 550),p( 551),p( 552)/  -51.65997,  -56.35155,  -58.34121/
      data p( 553),p( 554),p( 555)/  -61.97957,  -62.79452,  -65.57616/
      data p( 556),p( 557),p( 558)/  -65.42081,  -67.25972,  -66.25433/
      data p( 559),p( 560),p( 561)/  -67.30016,  -65.54189,  -65.73991/
      data p( 562),p( 563),p( 564)/  -62.96034,  -62.76423,  -60.06300/
      data p( 565),p( 566),p( 567)/  -59.15900,  -55.70300,  -54.30600/
      data p( 568),p( 569),p( 570)/  -50.31000,  -48.48400,  -43.95700/
      data p( 571),p( 572),p( 573)/  -41.65600,  -35.49900,   -6.56700/
      data p( 574),p( 575),p( 576)/  -14.92300,  -25.72800,  -32.68600/
      data p( 577),p( 578),p( 579)/  -42.29311,  -47.25741,  -54.18311/
      data p( 580),p( 581),p( 582)/  -56.34243,  -61.16735,  -62.20930/
      data p( 583),p( 584),p( 585)/  -65.99953,  -65.90778,  -68.89631/
      data p( 586),p( 587),p( 588)/  -67.87716,  -70.00404,  -68.41493/
      data p( 589),p( 590),p( 591)/  -69.55943,  -67.32168,  -68.12842/
      data p( 592),p( 593),p( 594)/  -65.41000,  -65.70920,  -62.46842/
      data p( 595),p( 596),p( 597)/  -62.04289,  -58.60414,  -57.22206/
      data p( 598),p( 599),p( 600)/  -53.39800,  -51.77735,  -46.12800/
      data p( 601),p( 602),p( 603)/  -42.06600,   -4.74100,  -15.90100/
      data p( 604),p( 605),p( 606)/  -23.98600,  -34.12100,  -39.99800/
      data p( 607),p( 608),p( 609)/  -47.34800,  -51.99635,  -56.68930/
      data p( 610),p( 611),p( 612)/  -58.83473,  -62.65291,  -63.72130/
      data p( 613),p( 614),p( 615)/  -66.87669,  -67.08293,  -69.32092/
      data p( 616),p( 617),p( 618)/  -68.90471,  -70.13683,  -68.58650/
      data p( 619),p( 620),p( 621)/  -69.70385,  -68.05401,  -68.46420/
      data p( 622),p( 623),p( 624)/  -66.20290,  -65.87415,  -63.66207/
      data p( 625),p( 626),p( 627)/  -62.48799,  -59.06775,  -57.98275/
      data p( 628),p( 629),p( 630)/  -52.94600,  -49.49000,  -44.39500/
      data p( 631),p( 632),p( 633)/   -8.37400,  -17.00000,  -27.76800/
      data p( 634),p( 635),p( 636)/  -33.72900,  -42.24300,  -46.91000/
      data p( 637),p( 638),p( 639)/  -54.42473,  -56.41056,  -61.62130/
      data p( 640),p( 641),p( 642)/  -62.65376,  -66.97696,  -67.09364/
      data p( 643),p( 644),p( 645)/  -70.56033,  -69.90490,  -72.58556/
      data p( 646),p( 647),p( 648)/  -71.29713,  -73.42201,  -71.85591/
      data p( 649),p( 650),p( 651)/  -73.21289,  -71.21414,  -71.86207/
      data p( 652),p( 653),p( 654)/  -69.48800,  -69.44775,  -66.30274/
      data p( 655),p( 656),p( 657)/  -65.62344,  -61.00400,  -58.39500/
      data p( 658),p( 659),p( 660)/  -53.38400,  -50.04900,   -6.39900/
      data p( 661),p( 662),p( 663)/  -18.05200,  -24.96400,  -33.82300/
      data p( 664),p( 665),p( 666)/  -39.52100,  -47.05600,  -51.82100/
      data p( 667),p( 668),p( 669)/  -56.64376,  -58.87696,  -63.08062/
      data p( 670),p( 671),p( 672)/  -64.34032,  -67.89219,  -68.22946/
      data p( 673),p( 674),p( 675)/  -70.95628,  -70.85960,  -73.03246/
      data p( 676),p( 677),p( 678)/  -72.28958,  -73.91618,  -72.81621/
      data p( 679),p( 680),p( 681)/  -73.63599,  -72.11797,  -72.53275/
      data p( 682),p( 683),p( 684)/  -70.32344,  -69.88009,  -66.08000/
      data p( 685),p( 686),p( 687)/  -63.51900,  -59.40100,  -56.28100/
      data p( 688),p( 689),p( 690)/  -51.64200,  -47.29200,  -32.91900/
      data p( 691),p( 692),p( 693)/  -41.72200,  -46.49100,  -54.14800/
      data p( 694),p( 695),p( 696)/  -56.29748,  -61.94000,  -63.09200/
      data p( 697),p( 698),p( 699)/  -67.89443,  -68.21628,  -72.21262/
      data p( 700),p( 701),p( 702)/  -72.16882,  -75.25156,  -74.59905/
      data p( 703),p( 704),p( 705)/  -77.02567,  -75.91694,  -77.75941/
      data p( 706),p( 707),p( 708)/  -76.38908,  -77.59344,  -75.34009/
      data p( 709),p( 710),p( 711)/  -75.94980,  -72.42860,  -70.54095/
      data p( 712),p( 713),p( 714)/  -66.58249,  -63.87814,  -59.59700/
      data p( 715),p( 716),p( 717)/  -56.43000,  -50.88800,  -47.19900/
      data p( 718),p( 719),p( 720)/  -32.79800,  -38.89200,  -46.40900/
      data p( 721),p( 722),p( 723)/  -51.59000,  -56.59200,  -59.15277/
      data p( 724),p( 725),p( 726)/  -63.53264,  -65.30596,  -69.13882/
      data p( 727),p( 728),p( 729)/  -70.28869,  -73.23393,  -73.45190/
      data p( 730),p( 731),p( 732)/  -76.06798,  -75.88885,  -77.97437/
      data p( 733),p( 734),p( 735)/  -77.49595,  -79.00911,  -77.77631/
      data p( 736),p( 737),p( 738)/  -78.61060,  -75.63995,  -73.85749/
      data p( 739),p( 740),p( 741)/  -70.73215,  -68.56982,  -64.61308/
      data p( 742),p( 743),p( 744)/  -61.51092,  -56.58336,  -53.00200/
      data p( 745),p( 746),p( 747)/  -47.80400,  -32.30400,  -40.97600/
      data p( 748),p( 749),p( 750)/  -46.10000,  -54.11277,  -56.88529/
      data p( 751),p( 752),p( 753)/  -62.16956,  -64.24160,  -68.97871/
      data p( 754),p( 755),p( 756)/  -70.17141,  -74.15971,  -74.44221/
      data p( 757),p( 758),p( 759)/  -77.89335,  -77.69365,  -80.58857/
      data p( 760),p( 761),p( 762)/  -79.98184,  -82.43104,  -81.48061/
      data p( 763),p( 764),p( 765)/  -83.26595,  -80.70999,  -79.69215/
      data p( 766),p( 767),p( 768)/  -76.72482,  -74.96309,  -71.31293/
      data p( 769),p( 770),p( 771)/  -68.78826,  -64.02600,  -61.14100/
      data p( 772),p( 773),p( 774)/  -56.03900,  -53.03000,  -47.91600/
      data p( 775),p( 776),p( 777)/  -32.30400,  -38.11700,  -46.23400/
      data p( 778),p( 779),p( 780)/  -51.72551,  -57.22242,  -60.48055/
      data p( 781),p( 782),p( 783)/  -64.82583,  -66.93577,  -70.79659/
      data p( 784),p( 785),p( 786)/  -72.17278,  -75.45644,  -76.18903/
      data p( 787),p( 788),p( 789)/  -79.07270,  -79.75015,  -82.16769/
      data p( 790),p( 791),p( 792)/  -82.74732,  -84.59505,  -82.60622/
      data p( 793),p( 794),p( 795)/  -81.71070,  -79.35495,  -77.74792/
      data p( 796),p( 797),p( 798)/  -74.77526,  -72.62601,  -68.55112/
      data p( 799),p( 800),p( 801)/  -65.83868,  -61.21409,  -58.36474/
      data p( 802),p( 803),p( 804)/  -54.30278,  -50.84036,  -46.69600/
      data p( 805),p( 806),p( 807)/  -43.59764,  -37.99600,  -31.69900/
      data p( 808),p( 809),p( 810)/  -40.69700,  -46.64900,  -54.39000/
      data p( 811),p( 812),p( 813)/  -57.97477,  -63.17451,  -65.47743/
      data p( 814),p( 815),p( 816)/  -70.30489,  -71.52654,  -76.00873/
      data p( 817),p( 818),p( 819)/  -76.79699,  -80.64429,  -81.10267/
      data p( 820),p( 821),p( 822)/  -84.52157,  -84.87836,  -87.91967/
      data p( 823),p( 824),p( 825)/  -86.20705,  -85.94186,  -83.63898/
      data p( 826),p( 827),p( 828)/  -82.87511,  -80.08761,  -78.84177/
      data p( 829),p( 830),p( 831)/  -75.11733,  -72.95416,  -68.79198/
      data p( 832),p( 833),p( 834)/  -66.62866,  -62.11663,  -60.21947/
      data p( 835),p( 836),p( 837)/  -55.40765,  -53.07765,  -47.55300/
      data p( 838),p( 839),p( 840)/  -44.40400,  -46.92900,  -52.62900/
      data p( 841),p( 842),p( 843)/  -58.35743,  -61.16500,  -66.01617/
      data p( 844),p( 845),p( 846)/  -68.19274,  -72.32810,  -74.15831/
      data p( 847),p( 848),p( 849)/  -77.84767,  -79.28156,  -83.01675/
      data p( 850),p( 851),p( 852)/  -84.29707,  -87.70210,  -86.48786/
      data p( 853),p( 854),p( 855)/  -86.34630,  -84.81548,  -84.22421/
      data p( 856),p( 857),p( 858)/  -82.34965,  -81.20419,  -78.34070/
      data p( 859),p( 860),p( 861)/  -76.26046,  -72.45203,  -70.20229/
      data p( 862),p( 863),p( 864)/  -67.29447,  -64.91265,  -61.89265/
      data p( 865),p( 866),p( 867)/  -58.74000,  -54.53900,  -51.14800/
      data p( 868),p( 869),p( 870)/  -46.37000,  -47.35700,  -55.37700/
      data p( 871),p( 872),p( 873)/  -58.85617,  -64.19273,  -66.46011/
      data p( 874),p( 875),p( 876)/  -71.49200,  -73.15467,  -77.80503/
      data p( 877),p( 878),p( 879)/  -79.34784,  -83.62377,  -84.86942/
      data p( 880),p( 881),p( 882)/  -88.76794,  -87.89114,  -88.45457/
      data p( 883),p( 884),p( 885)/  -87.11739,  -87.26630,  -85.65763/
      data p( 886),p( 887),p( 888)/  -85.44065,  -82.94886,  -81.27623/
      data p( 889),p( 890),p( 891)/  -77.76942,  -76.60447,  -73.45765/
      data p( 892),p( 893),p( 894)/  -71.74265,  -68.37439,  -66.34100/
      data p( 895),p( 896),p( 897)/  -62.36400,  -59.69900,  -55.08900/
      data p( 898),p( 899),p( 900)/  -51.90300,  -47.46000,  -52.97400/
      data p( 901),p( 902),p( 903)/  -58.96011,  -61.87900,  -67.15467/
      data p( 904),p( 905),p( 906)/  -69.82703,  -74.18284,  -76.42400/
      data p( 907),p( 908),p( 909)/  -80.57841,  -82.65694,  -86.63775/
      data p( 910),p( 911),p( 912)/  -86.44896,  -87.20875,  -86.36489/
      data p( 913),p( 914),p( 915)/  -86.78246,  -85.60423,  -85.60695/
      data p( 916),p( 917),p( 918)/  -83.52642,  -82.32742,  -79.93948/
      data p( 919),p( 920),p( 921)/  -78.94265,  -76.34764,  -75.31939/
      data p( 922),p( 923),p( 924)/  -72.22853,  -70.85500,  -66.89101/
      data p( 925),p( 926),p( 927)/  -64.91600,  -60.53800,  -58.09700/
      data p( 928),p( 929),p( 930)/  -53.39300,  -47.74800,  -55.80600/
      data p( 931),p( 932),p( 933)/  -59.06600,  -64.55703,  -67.69461/
      data p( 934),p( 935),p( 936)/  -72.70056,  -75.00336,  -80.16794/
      data p( 937),p( 938),p( 939)/  -82.20364,  -86.80547,  -86.80386/
      data p( 940),p( 941),p( 942)/  -88.41034,  -87.70808,  -88.79102/
      data p( 943),p( 944),p( 945)/  -87.54083,  -88.11201,  -85.96608/
      data p( 946),p( 947),p( 948)/  -86.18447,  -83.51165,  -83.55765/
      data p( 949),p( 950),p( 951)/  -80.84939,  -80.33353,  -77.34000/
      data p( 952),p( 953),p( 954)/  -76.25742,  -72.94089,  -71.18501/
      data p( 955),p( 956),p( 957)/  -67.24500,  -65.45600,  -61.00400/
      data p( 958),p( 959),p( 960)/  -58.83300,  -53.99900,  -47.56200/
      data p( 961),p( 962),p( 963)/  -53.20700,  -59.12200,  -62.56800/
      data p( 964),p( 965),p( 966)/  -67.49336,  -71.20728,  -75.98363/
      data p( 967),p( 968),p( 969)/  -78.93511,  -83.60300,  -84.15459/
      data p( 970),p( 971),p( 972)/  -86.01745,  -85.81779,  -87.22057/
      data p( 973),p( 974),p( 975)/  -86.42802,  -87.32331,  -86.01639/
      data p( 976),p( 977),p( 978)/  -86.33609,  -84.56760,  -84.59939/
      data p( 979),p( 980),p( 981)/  -82.48853,  -82.28999,  -79.77742/
      data p( 982),p( 983),p( 984)/  -79.10088,  -75.93541,  -74.86700/
      data p( 985),p( 986),p( 987)/  -71.36201,  -69.81500,  -65.91300/
      data p( 988),p( 989),p( 990)/  -63.96600,  -59.72700,  -57.49200/
      data p( 991),p( 992),p( 993)/  -47.33900,  -55.49800,  -59.51300/
      data p( 994),p( 995),p( 996)/  -65.40900,  -68.57897,  -74.40800/
      data p( 997),p( 998),p( 999)/  -77.26601,  -82.56802,  -83.45000/
      data p(1000),p(1001),p(1002)/  -86.07220,  -86.11238,  -88.22448/
      data p(1003),p(1004),p(1005)/  -87.61697,  -89.21880,  -87.94959/
      data p(1006),p(1007),p(1008)/  -89.09785,  -87.25893,  -88.09125/
      data p(1009),p(1010),p(1011)/  -85.92999,  -86.32442,  -83.92088/
      data p(1012),p(1013),p(1014)/  -83.65541,  -80.85221,  -80.13997/
      data p(1015),p(1016),p(1017)/  -76.79201,  -75.86700,  -72.15401/
      data p(1018),p(1019),p(1020)/  -70.79401,  -66.77901,  -65.05600/
      data p(1021),p(1022),p(1023)/  -60.74300,  -58.65600,  -47.15200/
      data p(1024),p(1025),p(1026)/  -53.21600,  -59.10300,  -63.36000/
      data p(1027),p(1028),p(1029)/  -69.17300,  -72.93800,  -78.34000/
      data p(1030),p(1031),p(1032)/  -79.62576,  -82.58938,  -83.16710/
      data p(1033),p(1034),p(1035)/  -85.57439,  -85.58880,  -87.40810/
      data p(1036),p(1037),p(1038)/  -86.77532,  -88.02228,  -86.95000/
      data p(1039),p(1040),p(1041)/  -87.84691,  -86.36382,  -86.86130/
      data p(1042),p(1043),p(1044)/  -85.01673,  -85.01221,  -82.94997/
      data p(1045),p(1046),p(1047)/  -82.28800,  -79.53700,  -78.78600/
      data p(1048),p(1049),p(1050)/  -75.59400,  -74.40338,  -71.06100/
      data p(1051),p(1052),p(1053)/  -69.53600,  -65.73600,  -63.93800/
      data p(1054),p(1055),p(1056)/  -59.82100,  -57.67800,  -47.05900/
      data p(1057),p(1058),p(1059)/  -55.49800,  -59.69900,  -66.35001/
      data p(1060),p(1061),p(1062)/  -70.15100,  -76.17576,  -77.79938/
      data p(1063),p(1064),p(1065)/  -81.30009,  -82.18780,  -85.22741/
      data p(1066),p(1067),p(1068)/  -85.42811,  -87.92584,  -87.47920/
      data p(1069),p(1070),p(1071)/  -89.39089,  -88.41363,  -89.90491/
      data p(1072),p(1073),p(1074)/  -88.37227,  -89.52173,  -87.60371/
      data p(1075),p(1076),p(1077)/  -88.34998,  -86.02910,  -86.33711/
      data p(1078),p(1079),p(1080)/  -83.69347,  -83.49416,  -80.40338/
      data p(1081),p(1082),p(1083)/  -79.96103,  -76.53201,  -75.46599/
      data p(1084),p(1085),p(1086)/  -72.02300,  -70.76601,  -66.90000/
      data p(1087),p(1088),p(1089)/  -65.39101,  -61.23600,  -53.30000/
      data p(1090),p(1091),p(1092)/  -60.10000,  -64.57101,  -70.79401/
      data p(1093),p(1094),p(1095)/  -72.88009,  -76.75780,  -78.18085/
      data p(1096),p(1097),p(1098)/  -81.22427,  -81.97147,  -84.79160/
      data p(1099),p(1100),p(1101)/  -85.11224,  -87.06838,  -86.93965/
      data p(1102),p(1103),p(1104)/  -88.40527,  -87.60355,  -88.71966/
      data p(1105),p(1106),p(1107)/  -87.45753,  -88.21742,  -86.62508/
      data p(1108),p(1109),p(1110)/  -87.03347,  -84.94488,  -84.98738/
      data p(1111),p(1112),p(1113)/  -82.56802,  -82.26565,  -79.56599/
      data p(1114),p(1115),p(1116)/  -78.55656,  -75.64792,  -74.65824/
      data p(1117),p(1118),p(1119)/  -71.42700,  -69.95500,  -66.57401/
      data p(1120),p(1121),p(1122)/  -64.70200,  -61.01300,  -58.79600/
      data p(1123),p(1124),p(1125)/  -56.10400,  -60.60300,  -67.46001/
      data p(1126),p(1127),p(1128)/  -69.85300,  -74.30505,  -75.74774/
      data p(1129),p(1130),p(1131)/  -79.38447,  -80.64971,  -83.97596/
      data p(1132),p(1133),p(1134)/  -84.33018,  -87.13380,  -86.98827/
      data p(1135),p(1136),p(1137)/  -89.25257,  -88.50536,  -90.34972/
      data p(1138),p(1139),p(1140)/  -89.25423,  -90.58105,  -89.04993/
      data p(1141),p(1142),p(1143)/  -90.02132,  -88.09087,  -88.71973/
      data p(1144),p(1145),p(1146)/  -86.42564,  -86.70891,  -83.90656/
      data p(1147),p(1148),p(1149)/  -83.97292,  -81.05824,  -80.57401/
      data p(1150),p(1151),p(1152)/  -77.31058,  -76.71011,  -73.35778/
      data p(1153),p(1154),p(1155)/  -72.32677,  -68.52552,  -67.29054/
      data p(1156),p(1157),p(1158)/  -63.09900,  -61.49700,  -53.80300/
      data p(1159),p(1160),p(1161)/  -60.91000,  -64.13425,  -68.40900/
      data p(1162),p(1163),p(1164)/  -70.13447,  -74.59972,  -76.06727/
      data p(1165),p(1166),p(1167)/  -79.48118,  -80.61043,  -83.56227/
      data p(1168),p(1169),p(1170)/  -84.09557,  -86.48541,  -86.47172/
      data p(1171),p(1172),p(1173)/  -88.38882,  -87.99512,  -89.36639/
      data p(1174),p(1175),p(1176)/  -88.56946,  -89.53675,  -88.24974/
      data p(1177),p(1178),p(1179)/  -88.94301,  -87.23009,  -87.70357/
      data p(1180),p(1181),p(1182)/  -85.73330,  -85.83824,  -83.57633/
      data p(1183),p(1184),p(1185)/  -83.42558,  -80.87611,  -80.47977/
      data p(1186),p(1187),p(1188)/  -77.81277,  -76.99352,  -74.36054/
      data p(1189),p(1190),p(1191)/  -72.97514,  -69.99717,  -68.21571/
      data p(1192),p(1193),p(1194)/  -62.48554,  -57.43600,  -51.54900/
      data p(1195),p(1196),p(1197)/  -56.86400,  -59.56000,  -64.74800/
      data p(1198),p(1199),p(1200)/  -66.94601,  -71.55227,  -73.22436/
      data p(1201),p(1202),p(1203)/  -77.42533,  -78.55595,  -82.00375/
      data p(1204),p(1205),p(1206)/  -82.63573,  -85.83466,  -85.94391/
      data p(1207),p(1208),p(1209)/  -88.65883,  -88.33043,  -90.55814/
      data p(1210),p(1211),p(1212)/  -90.03264,  -91.52473,  -90.39797/
      data p(1213),p(1214),p(1215)/  -91.65311,  -90.06719,  -91.10330/
      data p(1216),p(1217),p(1218)/  -89.20277,  -89.94492,  -87.81947/
      data p(1219),p(1220),p(1221)/  -88.23611,  -85.89777,  -86.01978/
      data p(1222),p(1223),p(1224)/  -83.50752,  -83.33614,  -80.63013/
      data p(1225),p(1226),p(1227)/  -80.24617,  -77.38931,  -76.62054/
      data p(1228),p(1229),p(1230)/  -70.96671,  -66.63574,  -60.79900/
      data p(1231),p(1232),p(1233)/  -56.50400,  -50.49600,  -55.77800/
      data p(1234),p(1235),p(1236)/  -59.34800,  -63.78070,  -66.35700/
      data p(1237),p(1238),p(1239)/  -70.65401,  -72.50700,  -76.25573/
      data p(1240),p(1241),p(1242)/  -77.53500,  -80.84400,  -81.60386/
      data p(1243),p(1244),p(1245)/  -84.41389,  -84.67664,  -87.00264/
      data p(1246),p(1247),p(1248)/  -86.81781,  -88.64134,  -87.99648/
      data p(1249),p(1250),p(1251)/  -89.47328,  -88.42270,  -89.59290/
      data p(1252),p(1253),p(1254)/  -88.32852,  -89.22250,  -87.61862/
      data p(1255),p(1256),p(1257)/  -88.26109,  -86.39777,  -86.70852/
      data p(1258),p(1259),p(1260)/  -84.61008,  -84.62614,  -82.39394/
      data p(1261),p(1262),p(1263)/  -82.02131,  -79.72353,  -78.95672/
      data p(1264),p(1265),p(1266)/  -74.00574,  -69.70559,  -64.59000/
      data p(1267),p(1268),p(1269)/  -60.25800,  -54.99500,  -50.57100/
      data p(1270),p(1271),p(1272)/  -58.03000,  -60.51300,  -65.68258/
      data p(1273),p(1274),p(1275)/  -67.57384,  -72.27725,  -73.47569/
      data p(1276),p(1277),p(1278)/  -77.25660,  -78.31400,  -81.91901/
      data p(1279),p(1280),p(1281)/  -82.36398,  -85.30597,  -85.10670/
      data p(1282),p(1283),p(1284)/  -87.72326,  -87.18027,  -89.40488/
      data p(1285),p(1286),p(1287)/  -88.55730,  -90.31107,  -89.16916/
      data p(1288),p(1289),p(1290)/  -90.52307,  -89.02779,  -90.07030/
      data p(1291),p(1292),p(1293)/  -88.28952,  -88.99364,  -87.00564/
      data p(1294),p(1295),p(1296)/  -87.35294,  -85.21131,  -85.20953/
      data p(1297),p(1298),p(1299)/  -82.95972,  -82.39944,  -77.82559/
      data p(1300),p(1301),p(1302)/  -74.42342,  -69.55952,  -65.93100/
      data p(1303),p(1304),p(1305)/  -60.79900,  -57.10100,  -51.80000/
      data p(1306),p(1307),p(1308)/  -47.97200,  -52.82400,  -57.57409/
      data p(1309),p(1310),p(1311)/  -60.34800,  -64.94701,  -67.09600/
      data p(1312),p(1313),p(1314)/  -71.12492,  -72.79601,  -76.45901/
      data p(1315),p(1316),p(1317)/  -77.56082,  -80.43665,  -80.69045/
      data p(1318),p(1319),p(1320)/  -83.66600,  -83.78989,  -86.28795/
      data p(1321),p(1322),p(1323)/  -86.07706,  -87.93494,  -87.36349/
      data p(1324),p(1325),p(1326)/  -88.84202,  -87.91496,  -88.98708/
      data p(1327),p(1328),p(1329)/  -87.74183,  -88.50358,  -86.93258/
      data p(1330),p(1331),p(1332)/  -87.44476,  -85.70254,  -85.87772/
      data p(1333),p(1334),p(1335)/  -83.94944,  -83.78759,  -79.49825/
      data p(1336),p(1337),p(1338)/  -76.50112,  -72.29914,  -68.84354/
      data p(1339),p(1340),p(1341)/  -64.07700,  -60.70500,  -55.72200/
      data p(1342),p(1343),p(1344)/  -52.09800,  -46.93800,  -51.72100/
      data p(1345),p(1346),p(1347)/  -54.36900,  -59.92733,  -62.05348/
      data p(1348),p(1349),p(1350)/  -66.93201,  -68.43000,  -72.90100/
      data p(1351),p(1352),p(1353)/  -73.99381,  -77.71368,  -78.66066/
      data p(1354),p(1355),p(1356)/  -81.82989,  -82.54294,  -85.18661/
      data p(1357),p(1358),p(1359)/  -85.25894,  -87.65751,  -87.18948/
      data p(1360),p(1361),p(1362)/  -89.17297,  -88.32465,  -89.86081/
      data p(1363),p(1364),p(1365)/  -88.69736,  -89.88180,  -88.41561/
      data p(1366),p(1367),p(1368)/  -89.27953,  -87.64830,  -88.12444/
      data p(1369),p(1370),p(1371)/  -86.43565,  -86.42445,  -82.37859/
      data p(1372),p(1373),p(1374)/  -80.11914,  -75.64954,  -72.99590/
      data p(1375),p(1376),p(1377)/  -68.32854,  -65.48125,  -60.65000/
      data p(1378),p(1379),p(1380)/  -57.53800,  -52.47100,  -49.09000/
      data p(1381),p(1382),p(1383)/  -43.77100,  -46.26600,  -51.66483/
      data p(1384),p(1385),p(1386)/  -54.56600,  -59.67200,  -62.49006/
      data p(1387),p(1388),p(1389)/  -66.47189,  -68.41368,  -72.31105/
      data p(1390),p(1391),p(1392)/  -73.88776,  -77.14294,  -78.13190/
      data p(1393),p(1394),p(1395)/  -81.04913,  -81.74257,  -84.09073/
      data p(1396),p(1397),p(1398)/  -84.34869,  -86.23994,  -85.93225/
      data p(1399),p(1400),p(1401)/  -87.50140,  -86.90263,  -88.06322/
      data p(1402),p(1403),p(1404)/  -87.16007,  -88.07566,  -86.89588/
      data p(1405),p(1406),p(1407)/  -87.58660,  -86.34414,  -86.55116/
      data p(1408),p(1409),p(1410)/  -82.89314,  -80.70657,  -77.05590/
      data p(1411),p(1412),p(1413)/  -74.47855,  -70.52125,  -67.69139/
      data p(1414),p(1415),p(1416)/  -63.31609,  -60.18547,  -55.73870/
      data p(1417),p(1418),p(1419)/  -52.28994,  -47.59977,  -44.04100/
      data p(1420),p(1421),p(1422)/  -39.15100,  -35.39700,  -45.69800/
      data p(1423),p(1424),p(1425)/  -48.70800,  -54.32500,  -56.95200/
      data p(1426),p(1427),p(1428)/  -62.00000,  -64.22471,  -68.88776/
      data p(1429),p(1430),p(1431)/  -70.34092,  -74.27700,  -75.59100/
      data p(1432),p(1433),p(1434)/  -79.09460,  -79.53073,  -82.67554/
      data p(1435),p(1436),p(1437)/  -82.78994,  -85.40973,  -85.06985/
      data p(1438),p(1439),p(1440)/  -87.27122,  -86.69340,  -88.43961/
      data p(1441),p(1442),p(1443)/  -87.55822,  -88.95455,  -87.85594/
      data p(1444),p(1445),p(1446)/  -88.89236,  -87.72678,  -88.26717/
      data p(1447),p(1448),p(1449)/  -84.91929,  -83.27603,  -79.72988/
      data p(1450),p(1451),p(1452)/  -77.82802,  -73.94461,  -71.78049/
      data p(1453),p(1454),p(1455)/  -68.07003,  -65.10523,  -61.48557/
      data p(1456),p(1457),p(1458)/  -58.04849,  -53.59800,  -50.65500/
      data p(1459),p(1460),p(1461)/  -45.92300,  -42.70000,  -37.62300/
      data p(1462),p(1463),p(1464)/  -46.56500,  -49.77000,  -54.96700/
      data p(1465),p(1466),p(1467)/  -57.68700,  -62.40100,  -64.54300/
      data p(1468),p(1469),p(1470)/  -68.70700,  -70.30000,  -73.89500/
      data p(1471),p(1472),p(1473)/  -75.10600,  -78.09600,  -78.75973/
      data p(1474),p(1475),p(1476)/  -81.34985,  -81.67300,  -83.73339/
      data p(1477),p(1478),p(1479)/  -83.73161,  -85.32822,  -85.24137/
      data p(1480),p(1481),p(1482)/  -86.65594,  -86.02236,  -87.12667/
      data p(1483),p(1484),p(1485)/  -86.52943,  -87.23612,  -84.32577/
      data p(1486),p(1487),p(1488)/  -82.94299,  -80.03909,  -78.19086/
      data p(1489),p(1490),p(1491)/  -74.89987,  -72.99338,  -69.20986/
      data p(1492),p(1493),p(1494)/  -67.23557,  -63.16349,  -61.13400/
      data p(1495),p(1496),p(1497)/  -57.22200,  -54.43700,  -50.19800/
      data p(1498),p(1499),p(1500)/  -47.08700,  -42.47600,  -39.00200/
      data p(1501),p(1502),p(1503)/  -44.00400,  -49.70500,  -52.47100/
      data p(1504),p(1505),p(1506)/  -57.74300,  -60.07200,  -64.72000/
      data p(1507),p(1508),p(1509)/  -66.56500,  -70.70000,  -71.95800/
      data p(1510),p(1511),p(1512)/  -75.57201,  -76.30000,  -79.46201/
      data p(1513),p(1514),p(1515)/  -79.71339,  -82.44701,  -82.39101/
      data p(1516),p(1517),p(1518)/  -84.74137,  -84.63036,  -86.49519/
      data p(1519),p(1520),p(1521)/  -85.90457,  -87.57386,  -86.95812/
      data p(1522),p(1523),p(1524)/  -88.08762,  -85.44491,  -84.54263/
      data p(1525),p(1526),p(1527)/  -81.61642,  -80.44131,  -77.10173/
      data p(1528),p(1529),p(1530)/  -75.74003,  -72.18056,  -70.42584/
      data p(1531),p(1532),p(1533)/  -66.79816,  -64.99368,  -61.44100/
      data p(1534),p(1535),p(1536)/  -59.26200,  -55.34900,  -52.79700/
      data p(1537),p(1538),p(1539)/  -48.40000,  -45.40100,  -40.66900/
      data p(1540),p(1541),p(1542)/  -41.57900,  -45.03800,  -50.33800/
      data p(1543),p(1544),p(1545)/  -53.13200,  -57.91100,  -60.25800/
      data p(1546),p(1547),p(1548)/  -64.43100,  -66.32201,  -69.99200/
      data p(1549),p(1550),p(1551)/  -71.37100,  -74.46339,  -75.33900/
      data p(1552),p(1553),p(1554)/  -78.05901,  -78.55100,  -80.91036/
      data p(1555),p(1556),p(1557)/  -81.36885,  -83.20258,  -83.13686/
      data p(1558),p(1559),p(1560)/  -84.82912,  -84.69962,  -86.02558/
      data p(1561),p(1562),p(1563)/  -83.79732,  -83.07787,  -80.75996/
      data p(1564),p(1565),p(1566)/  -79.63631,  -76.76626,  -75.47057/
      data p(1567),p(1568),p(1569)/  -72.48584,  -70.98817,  -68.00368/
      data p(1570),p(1571),p(1572)/  -66.85530,  -63.71400,  -61.80500/
      data p(1573),p(1574),p(1575)/  -58.32100,  -55.89900,  -52.05200/
      data p(1576),p(1577),p(1578)/  -49.21100,  -44.91700,  -41.70300/
      data p(1579),p(1580),p(1581)/  -53.03000,  -55.42400,  -60.18400/
      data p(1582),p(1583),p(1584)/  -62.17300,  -66.34100,  -67.90340/
      data p(1585),p(1586),p(1587)/  -71.61301,  -72.46101,  -75.78101/
      data p(1588),p(1589),p(1590)/  -76.15900,  -79.15784,  -79.51257/
      data p(1591),p(1592),p(1593)/  -82.03700,  -82.04212,  -84.47735/
      data p(1594),p(1595),p(1596)/  -84.20258,  -85.95952,  -84.01179/
      data p(1597),p(1598),p(1599)/  -83.75748,  -81.44158,  -80.93551/
      data p(1600),p(1601),p(1602)/  -78.15625,  -77.41784,  -74.38520/
      data p(1603),p(1604),p(1605)/  -73.69368,  -70.95679,  -70.15786/
      data p(1606),p(1607),p(1608)/  -67.35210,  -65.68587,  -62.75516/
      data p(1609),p(1610),p(1611)/  -60.36100,  -56.57000,  -54.14800/
      data p(1612),p(1613),p(1614)/  -49.93700,  -47.14300,  -42.54100/
      data p(1615),p(1616),p(1617)/  -48.19500,  -52.94600,  -55.47000/
      data p(1618),p(1619),p(1620)/  -59.80200,  -61.71100,  -65.46500/
      data p(1621),p(1622),p(1623)/  -66.61100,  -70.21900,  -71.30785/
      data p(1624),p(1625),p(1626)/  -73.85600,  -75.03700,  -77.53772/
      data p(1627),p(1628),p(1629)/  -78.43025,  -80.47489,  -81.08585/
      data p(1630),p(1631),p(1632)/  -82.97043,  -81.42583,  -81.27854/
      data p(1633),p(1634),p(1635)/  -79.46373,  -79.05225,  -76.87825/
      data p(1636),p(1637),p(1638)/  -76.07585,  -73.60714,  -73.39921/
      data p(1639),p(1640),p(1641)/  -71.26808,  -70.68810,  -68.42101/
      data p(1642),p(1643),p(1644)/  -66.97716,  -64.21686,  -62.22400/
      data p(1645),p(1646),p(1647)/  -58.97300,  -56.70000,  -53.10400/
      data p(1648),p(1649),p(1650)/  -50.43100,  -46.30500,  -43.29600/
      data p(1651),p(1652),p(1653)/  -47.85100,  -50.40300,  -55.12600/
      data p(1654),p(1655),p(1656)/  -57.07300,  -61.46000,  -63.01600/
      data p(1657),p(1658),p(1659)/  -66.78800,  -67.95555,  -71.22200/
      data p(1660),p(1661),p(1662)/  -72.37521,  -75.45939,  -75.94608/
      data p(1663),p(1664),p(1665)/  -78.99695,  -79.52763,  -81.97637/
      data p(1666),p(1667),p(1668)/  -80.66215,  -81.00573,  -79.27640/
      data p(1669),p(1670),p(1671)/  -79.34660,  -77.14677,  -77.06113/
      data p(1672),p(1673),p(1674)/  -74.58625,  -74.77265,  -72.56905/
      data p(1675),p(1676),p(1677)/  -72.46529,  -70.20116,  -69.37186/
      data p(1678),p(1679),p(1680)/  -66.73734,  -65.21581,  -62.22400/
      data p(1681),p(1682),p(1683)/  -60.41700,  -56.97900,  -54.75300/
      data p(1684),p(1685),p(1686)/  -50.89700,  -48.17700,  -43.79900/
      data p(1687),p(1688),p(1689)/  -42.70000,  -47.59900,  -50.00300/
      data p(1690),p(1691),p(1692)/  -54.28700,  -56.35500,  -60.35100/
      data p(1693),p(1694),p(1695)/  -61.99100,  -65.35500,  -66.98939/
      data p(1696),p(1697),p(1698)/  -69.96835,  -71.35240,  -74.25251/
      data p(1699),p(1700),p(1701)/  -75.66142,  -78.00211,  -77.12797/
      data p(1702),p(1703),p(1704)/  -77.55506,  -76.23926,  -76.45150/
      data p(1705),p(1706),p(1707)/  -74.80055,  -74.66294,  -72.89834/
      data p(1708),p(1709),p(1710)/  -73.37730,  -71.74796,  -71.82803/
      data p(1711),p(1712),p(1713)/  -70.09412,  -69.47134,  -67.21481/
      data p(1714),p(1715),p(1716)/  -66.05736,  -63.37200,  -61.77700/
      data p(1717),p(1718),p(1719)/  -58.64700,  -56.62600,  -53.10400/
      data p(1720),p(1721),p(1722)/  -50.56100,  -46.60300,  -43.73400/
      data p(1723),p(1724),p(1725)/  -49.30400,  -51.55800,  -55.91800/
      data p(1726),p(1727),p(1728)/  -57.67800,  -61.52900,  -63.14600/
      data p(1729),p(1730),p(1731)/  -66.85201,  -68.24251,  -71.92101/
      data p(1732),p(1733),p(1734)/  -72.94762,  -76.09808,  -75.36769/
      data p(1735),p(1736),p(1737)/  -76.28024,  -75.13763,  -75.77195/
      data p(1738),p(1739),p(1740)/  -74.19882,  -74.71709,  -72.89286/
      data p(1741),p(1742),p(1743)/  -73.71631,  -72.08012,  -72.54516/
      data p(1744),p(1745),p(1746)/  -70.83389,  -70.69990,  -68.57185/
      data p(1747),p(1748),p(1749)/  -67.95191,  -65.51598,  -64.29058/
      data p(1750),p(1751),p(1752)/  -61.48800,  -59.74600,  -56.46700/
      data p(1753),p(1754),p(1755)/  -54.39900,  -50.70100,  -48.10200/
      data p(1756),p(1757),p(1758)/  -43.90100,  -43.90100,  -48.41000/
      data p(1759),p(1760),p(1761)/  -50.72900,  -54.80900,  -56.95200/
      data p(1762),p(1763),p(1764)/  -60.78000,  -62.84800,  -66.24800/
      data p(1765),p(1766),p(1767)/  -67.83080,  -70.75891,  -70.51537/
      data p(1768),p(1769),p(1770)/  -71.49995,  -71.11569,  -71.63359/
      data p(1771),p(1772),p(1773)/  -70.72710,  -71.32369,  -70.15430/
      data p(1774),p(1775),p(1776)/  -71.25890,  -70.10074,  -70.77383/
      data p(1777),p(1778),p(1779)/  -69.47989,  -69.54242,  -67.84627/
      data p(1780),p(1781),p(1782)/  -67.47156,  -65.68448,  -64.60474/
      data p(1783),p(1784),p(1785)/  -62.08663,  -60.65900,  -57.70600/
      data p(1786),p(1787),p(1788)/  -55.84300,  -52.49900,  -50.09600/
      data p(1789),p(1790),p(1791)/  -46.34200,  -43.50100,  -43.04400/
      data p(1792),p(1793),p(1794)/  -45.46600,  -50.05200,  -52.32200/
      data p(1795),p(1796),p(1797)/  -56.75600,  -58.72800,  -62.67080/
      data p(1798),p(1799),p(1800)/  -64.38626,  -67.83337,  -67.68796/
      data p(1801),p(1802),p(1803)/  -69.32203,  -68.76322,  -70.12856/
      data p(1804),p(1805),p(1806)/  -69.15330,  -70.40040,  -69.16440/
      data p(1807),p(1808),p(1809)/  -70.53432,  -69.43239,  -70.41662/
      data p(1810),p(1811),p(1812)/  -69.17677,  -69.68159,  -68.06464/
      data p(1813),p(1814),p(1815)/  -68.19026,  -66.38987,  -65.97663/
      data p(1816),p(1817),p(1818)/  -63.62119,  -62.59337,  -59.94254/
      data p(1819),p(1820),p(1821)/  -58.47000,  -55.60679,  -53.40300/
      data p(1822),p(1823),p(1824)/  -49.85400,  -47.40400,  -43.37000/
      data p(1825),p(1826),p(1827)/  -37.39000,  -42.20600,  -45.04700/
      data p(1828),p(1829),p(1830)/  -49.48100,  -52.07100,  -56.03900/
      data p(1831),p(1832),p(1833)/  -58.43300,  -61.67426,  -62.08200/
      data p(1834),p(1835),p(1836)/  -63.63892,  -63.58322,  -65.02339/
      data p(1837),p(1838),p(1839)/  -64.64915,  -66.06240,  -65.47401/
      data p(1840),p(1841),p(1842)/  -66.89239,  -66.18662,  -67.33906/
      data p(1843),p(1844),p(1845)/  -66.39159,  -67.20573,  -66.04998/
      data p(1846),p(1847),p(1848)/  -66.38731,  -64.98979,  -64.90727/
      data p(1849),p(1850),p(1851)/  -63.07959,  -62.29254,  -60.08469/
      data p(1852),p(1853),p(1854)/  -58.80679,  -56.24830,  -54.52851/
      data p(1855),p(1856),p(1857)/  -51.40000,  -49.09900,  -45.50300/
      data p(1858),p(1859),p(1860)/  -42.80200,  -36.71000,  -39.62600/
      data p(1861),p(1862),p(1863)/  -44.60000,  -47.21700,  -51.75400/
      data p(1864),p(1865),p(1866)/  -53.86400,  -57.97400,  -58.25600/
      data p(1867),p(1868),p(1869)/  -60.47403,  -60.46036,  -62.61754/
      data p(1870),p(1871),p(1872)/  -62.21981,  -64.25910,  -63.39233/
      data p(1873),p(1874),p(1875)/  -65.28700,  -64.57049,  -66.06255/
      data p(1876),p(1877),p(1878)/  -65.20332,  -66.34572,  -65.17731/
      data p(1879),p(1880),p(1881)/  -65.95257,  -64.53129,  -64.93447/
      data p(1882),p(1883),p(1884)/  -63.29925,  -62.99900,  -60.93081/
      data p(1885),p(1886),p(1887)/  -60.11831,  -57.72851,  -56.49310/
      data p(1888),p(1889),p(1890)/  -53.65400,  -51.84700,  -48.50300/
      data p(1891),p(1892),p(1893)/  -46.30500,  -42.50400,  -31.21000/
      data p(1894),p(1895),p(1896)/  -36.25300,  -39.54200,  -44.10600/
      data p(1897),p(1898),p(1899)/  -46.88200,  -50.82800,  -51.88400/
      data p(1900),p(1901),p(1902)/  -54.00091,  -54.56400,  -56.64269/
      data p(1903),p(1904),p(1905)/  -56.81470,  -58.91133,  -58.68700/
      data p(1906),p(1907),p(1908)/  -60.72505,  -60.46255,  -62.03932/
      data p(1909),p(1910),p(1911)/  -61.50640,  -62.73831,  -61.99001/
      data p(1912),p(1913),p(1914)/  -62.93875,  -61.89485,  -62.55089/
      data p(1915),p(1916),p(1917)/  -61.31989,  -61.28194,  -59.80389/
      data p(1918),p(1919),p(1920)/  -59.21896,  -57.38364,  -56.26192/
      data p(1921),p(1922),p(1923)/  -53.87331,  -52.31932,  -49.37718/
      data p(1924),p(1925),p(1926)/  -47.46900,  -44.11600,  -41.60100/
      data p(1927),p(1928),p(1929)/  -30.96300,  -34.01800,  -39.13200/
      data p(1930),p(1931),p(1932)/  -41.68500,  -46.41900,  -47.31100/
      data p(1933),p(1934),p(1935)/  -50.07500,  -50.49400,  -53.23757/
      data p(1936),p(1937),p(1938)/  -53.41312,  -56.02200,  -55.74643/
      data p(1939),p(1940),p(1941)/  -58.16300,  -57.88900,  -59.84800/
      data p(1942),p(1943),p(1944)/  -59.36831,  -60.99400,  -60.17675/
      data p(1945),p(1946),p(1947)/  -61.59073,  -60.59660,  -61.57690/
      data p(1948),p(1949),p(1950)/  -60.37280,  -60.77192,  -59.31539/
      data p(1951),p(1952),p(1953)/  -59.26379,  -57.56003,  -56.95331/
      data p(1954),p(1955),p(1956)/  -54.70432,  -53.49718,  -50.99265/
      data p(1957),p(1958),p(1959)/  -49.70135,  -46.41600,  -44.40400/
      data p(1960),p(1961),p(1962)/  -40.84600,  -25.46000,  -30.60200/
      data p(1963),p(1964),p(1965)/  -33.89700,  -38.48000,  -39.96100/
      data p(1966),p(1967),p(1968)/  -42.63200,  -43.86600,  -46.48011/
      data p(1969),p(1970),p(1971)/  -47.34900,  -49.72770,  -50.28300/
      data p(1972),p(1973),p(1974)/  -52.58900,  -52.88800,  -54.76831/
      data p(1975),p(1976),p(1977)/  -54.75800,  -56.25675,  -56.11073/
      data p(1978),p(1979),p(1980)/  -57.46660,  -57.10190,  -58.07980/
      data p(1981),p(1982),p(1983)/  -57.31278,  -57.83654,  -56.74452/
      data p(1984),p(1985),p(1986)/  -56.88922,  -55.57896,  -55.17434/
      data p(1987),p(1988),p(1989)/  -53.39100,  -52.39188,  -50.34597/
      data p(1990),p(1991),p(1992)/  -49.06717,  -46.68650,  -44.74000/
      data p(1993),p(1994),p(1995)/  -41.72200,  -39.52300,  -36.17000/
      data p(1996),p(1997),p(1998)/  -33.30100,  -34.68900,  -37.96100/
      data p(1999),p(2000),p(2001)/  -39.00400,  -42.24600,  -42.84600/
      data p(2002),p(2003),p(2004)/  -45.90999,  -46.26651,  -49.18010/
      data p(2005),p(2006),p(2007)/  -49.31600,  -51.77000,  -51.66100/
      data p(2008),p(2009),p(2010)/  -53.79400,  -53.46800,  -55.30300/
      data p(2011),p(2012),p(2013)/  -54.81044,  -56.21600,  -55.43300/
      data p(2014),p(2015),p(2016)/  -56.39452,  -55.28400,  -55.85223/
      data p(2017),p(2018),p(2019)/  -54.48961,  -54.58384,  -52.89021/
      data p(2020),p(2021),p(2022)/  -52.44522,  -50.47293,  -49.78951/
      data p(2023),p(2024),p(2025)/  -47.41385,  -46.05968,  -43.28558/
      data p(2026),p(2027),p(2028)/  -41.50003,  -38.39600,  -36.40300/
      data p(2029),p(2030),p(2031)/  -26.37100,  -29.67300,  -31.32800/
      data p(2032),p(2033),p(2034)/  -34.54500,  -35.99500,  -38.77530/
      data p(2035),p(2036),p(2037)/  -39.91700,  -42.55376,  -43.24900/
      data p(2038),p(2039),p(2040)/  -45.81300,  -46.13700,  -48.46300/
      data p(2041),p(2042),p(2043)/  -48.63600,  -50.37500,  -50.21700/
      data p(2044),p(2045),p(2046)/  -51.73500,  -51.47452,  -52.59400/
      data p(2047),p(2048),p(2049)/  -52.00723,  -52.49000,  -51.47384/
      data p(2050),p(2051),p(2052)/  -51.72421,  -50.53322,  -50.36204/
      data p(2053),p(2054),p(2055)/  -48.93542,  -48.44109,  -46.43272/
      data p(2056),p(2057),p(2058)/  -45.29558,  -42.84003,  -41.39644/
      data p(2059),p(2060),p(2061)/  -38.61033,  -36.87800,  -33.80400/
      data p(2062),p(2063),p(2064)/  -24.27600,  -25.82100,  -29.46400/
      data p(2065),p(2066),p(2067)/  -30.65600,  -34.14700,  -34.90100/
      data p(2068),p(2069),p(2070)/  -38.20633,  -38.80980,  -41.89869/
      data p(2071),p(2072),p(2073)/  -42.22300,  -44.83900,  -44.93600/
      data p(2074),p(2075),p(2076)/  -47.23600,  -47.07800,  -48.97500/
      data p(2077),p(2078),p(2079)/  -48.59400,  -50.15200,  -49.58300/
      data p(2080),p(2081),p(2082)/  -50.68300,  -49.72300,  -50.44192/
      data p(2083),p(2084),p(2085)/  -49.30236,  -49.64328,  -48.25320/
      data p(2086),p(2087),p(2088)/  -48.24625,  -46.36561,  -45.70603/
      data p(2089),p(2090),p(2091)/  -43.38844,  -42.51133,  -39.90672/
      data p(2092),p(2093),p(2094)/  -38.66915,  -35.47854,  -34.29803/
      data p(2095),p(2096),p(2097)/  -17.24700,  -20.80900,  -22.62900/
      data p(2098),p(2099),p(2100)/  -26.11200,  -27.64700,  -30.69247/
      data p(2101),p(2102),p(2103)/  -31.85500,  -34.87200,  -35.76100/
      data p(2104),p(2105),p(2106)/  -38.35000,  -38.97100,  -41.40800/
      data p(2107),p(2108),p(2109)/  -41.65100,  -43.72200,  -43.67600/
      data p(2110),p(2111),p(2112)/  -45.27700,  -45.11200,  -46.32300/
      data p(2113),p(2114),p(2115)/  -45.78192,  -46.59236,  -45.84098/
      data p(2116),p(2117),p(2118)/  -46.51453,  -45.44624,  -45.80961/
      data p(2119),p(2120),p(2121)/  -44.22334,  -43.82143,  -41.92977/
      data p(2122),p(2123),p(2124)/  -41.21787,  -39.01815,  -37.97854/
      data p(2125),p(2126),p(2127)/  -35.56803,  -34.35015,  -31.70800/
      data p(2128),p(2129),p(2130)/  -15.07200,  -16.72200,  -20.56100/
      data p(2131),p(2132),p(2133)/  -21.91400,  -25.59200,  -26.49700/
      data p(2134),p(2135),p(2136)/  -29.96345,  -30.66831,  -33.93459/
      data p(2137),p(2138),p(2139)/  -34.42800,  -37.18700,  -37.45400/
      data p(2140),p(2141),p(2142)/  -39.93900,  -39.98000,  -41.96400/
      data p(2143),p(2144),p(2145)/  -41.87500,  -43.45584,  -42.89400/
      data p(2146),p(2147),p(2148)/  -44.38500,  -43.52452,  -44.53825/
      data p(2149),p(2150),p(2151)/  -43.67800,  -44.25452,  -42.80864/
      data p(2152),p(2153),p(2154)/  -42.99929,  -41.22054,  -41.13850/
      data p(2155),p(2156),p(2157)/  -38.98780,  -38.70803,  -36.39538/
      data p(2158),p(2159),p(2160)/  -35.88203,  -33.39584,  -32.43526/
      data p(2161),p(2162),p(2163)/  -29.69240,  -28.29641,  -11.56900/
      data p(2164),p(2165),p(2166)/  -13.50100,  -17.19300,  -18.66200/
      data p(2167),p(2168),p(2169)/  -21.99176,  -23.25700,  -26.28800/
      data p(2170),p(2171),p(2172)/  -27.34600,  -30.08000,  -30.92200/
      data p(2173),p(2174),p(2175)/  -33.27400,  -33.98900,  -36.17000/
      data p(2176),p(2177),p(2178)/  -36.25100,  -38.05200,  -37.95500/
      data p(2179),p(2180),p(2181)/  -39.45607,  -39.00380,  -40.22800/
      data p(2182),p(2183),p(2184)/  -39.69252,  -40.43600,  -39.16829/
      data p(2185),p(2186),p(2187)/  -39.71813,  -38.32915,  -38.45535/
      data p(2188),p(2189),p(2190)/  -36.70803,  -36.70906,  -34.83583/
      data p(2191),p(2192),p(2193)/  -34.53635,  -32.53186,  -31.69240/
      data p(2194),p(2195),p(2196)/  -29.45393,  -28.28342,  -25.82100/
      data p(2197),p(2198),p(2199)/  -24.41710,  -11.14600,  -12.64900/
      data p(2200),p(2201),p(2202)/  -16.46300,  -17.46500,  -21.07399/
      data p(2203),p(2204),p(2205)/  -21.89044,  -25.32613,  -25.82500/
      data p(2206),p(2207),p(2208)/  -28.87600,  -29.38600,  -31.94100/
      data p(2209),p(2210),p(2211)/  -32