c/* @(#)bsmain.f        1.7      latest revision 4/13/87 11:38:02 */
* %w%      latest revision %g% %u% */
 


      program bsmain
c   1=output,tape4=64,tape8=64,tape7,debug=output)
c      bound states via det(1-gv)=0
c      ridge32 version july 85
c      nes =no 2 body energies for t s, or or swtch for k-p potents
c      nwaves lt use non reltv reduced mass in sig pi channel
c               zero,  then reset
c      nifty(19)= 0     real*8 k2 in v2......dont redefine k2 in v2 as e change
c                 1      ..                do ...
c                 2     comple ....        do ...
       implicit real*8 (a-h,o-z)
       external detcal
      real*8 xeri(2),a(2,2),w(31),drr(10),dri(10),detri(2)
      common/det/ar11,alfanz
c     the open stmts are c'ed out so files spefied on run card
      open(5,file='runbs',status='unknown')
      open(6,file='outbs',status='unknown')
      open(7,status='scratch',access='sequential',
     1                      form='unformatted')
      i=signal(119,dumsub,1)
      write(6,879)
      write(6,879)
      write(6,879)
 879  format(' BOUND STATE calculation via det(1-GV)=0, ridge32, ',
     1'8/86, program bsdetr, detcal')
  10  read(5,*,end=233)drr(1),eim,(drr(i),i=2,7)
      read(5,*)tol,aitmax
  880 format(8f10.6)
      write(6,882)drr(1),eim,(drr(i),i=2,7)
      write(6,882)tol,aitmax
  882 format(8f15.6)
      maxit = aitmax
      nrei=2
       if(eim.ge.0.d+00)nrei=1
      l=0
c      calculate point coulomb b.e. for sch eqtn or kge
c      no dirac yet, i.e. no spin
c      fj=l +- 1/2....i.e.fj=l for these cases
      fj=l
      ekge=0.
 226  format(24h nbohr,e(se),e(kge1,2)= ,i3,3(1pe20.10),6h eta= ,
     1 e20.8)
c--------do loop over energies
      do 213 ix=1,7
      if(drr(ix).eq.0.d+00) go to 213
      xeri(1)=drr(ix)
      xeri(2)=eim/ix
      write(6,880)xeri(1),xeri(2)
      call snsqe(detcal,jac,2,2,xeri,detri,tol,-1,info,w,31,maxit)
	if (info .eq. 0)then
	   write(6,100)
 100       format(1h0,'improper input parameters in search')
	else if (info .eq. 1) then
		write(6,101) tol
 101       format(1h0,'normal termination,  relative error is at most',
     1 f13.10)
	     else if (info .eq. 2)then
		     write(6,102) maxit
 102                 format(1h0,'iterations exceed ',i4)
		  else if (info .eq. 3) then
			  write(6,103) tol
 103                      format(1h0,f13.10,'tol too small.no further 
     1improvement possible in xeri')
		       else if (info .eq. 4) then
				write(6,104)
 104                         format(1h0,'iteration not making good pr
     1ogress')
	endif
      nmin=l+1
      if(ix.ne.1) go to 211
      do 223 nbohr=nmin,4
      esch=(-ar11/2.)*(alfanz/nbohr)**2
      dri(nbohr)=esch
      eta=nbohr-(fj+.5d+00)+sqrt((fj+.5d+00)**2-(alfanz)**2)
      wron1=(alfanz/eta)**2
      eta=sqrt(1.+wron1)
      wron2=ar11*(-.5*wron1+wron1**2/8.-5.*wron1**3/16.)
      ekge=ar11*(1./eta-1.)
      write(6,226) nbohr,esch,ekge,wron2,eta
 223   continue
  211 epsev=(dri(ix)-xeri(1))
      gamev=-xeri(2)*2.
      esch=-epsev
      write(6,881) xeri(1),xeri(2),epsev,gamev,esch,xeri(2)
  881 format('0E(MeV)=',2f16.14,' eps,gamma(eV)=',2(6pf10.1)/
     1 ' Del E(ev)=',2(6pf10.1))
  213 continue
      go to 10
  233 write(6,841)
  841 format(' **************eof in main***********')
      stop
      end
      subroutine caxpy(n,ca,cx,incx,cy,incy)
c***begin prologue  caxpy
c***revision date  811015   (yymmdd)
c***category no.  f1a
c***keywords  blas,complex
c***date written  october 1979
c
c***changed to double precision june 1985 by m.s.
c
c***author lawson c. (jpl),hanson r. (sla),
c                            kincaid d. (u texas), krogh f. (jpl)
c***purpose
c  complex computation y = a*x + y
c***description
c                b l a s  subprogram
c    description of parameters
c
c     --input--
c        n  number of elements in input vector(s)
c       ca  complex scalar multiplier
c       cx  complex vector with n elements
c     incx  storage spacing between elements of cx
c       cy  complex vector with n elements
c     incy  storage spacing between elements of cy
c
c     --output--
c       cy  complex result (unchanged if n.le.0)
c
c     overwrite complex cy with complex  ca*cx + cy.
c     for i = 0 to n-1, replace  cy(ly+i*incy) with ca*cx(lx+i*incx) +
c       cy(ly+i*incy), where lx = 1 if incx .ge. 0, else lx = (-incx)*n
c       and ly is defined in a similar way using incy.
c
c
c***references
c  lawson c.l., hanson r.j., kincaid d.r., krogh f.t.,
c   *basic linear algebra subprograms for fortran usage*,
c  algorithm no. 539, transactions on mathematical software,
c  volume 5, number 3, september 1979, 308-323
c***routines called    (none)
c***end prologue  caxpy
c
      implicit real*8 (a-h,o-z)
      complex*16 cx(1),cy(1),ca
c***first executable statement  caxpy
      canorm = abs(dble(ca)) + abs(dimag(ca))
      if(n.le.0.or.canorm.eq.0.e0) return
      if(incx.eq.incy.and.incx.gt.0) go to 20
      kx = 1
      ky = 1
      if(incx.lt.0) kx = 1+(1-n)*incx
      if(incy.lt.0) ky = 1+(1-n)*incy
	  do 10 i = 1,n
	  cy(ky) = cy(ky) + ca*cx(kx)
	  kx = kx + incx
	  ky = ky + incy
   10 continue
      return
   20 continue
      ns = n*incx
	  do 30 i=1,ns,incx
	  cy(i) = ca*cx(i) + cy(i)
   30     continue
      return
      end
      subroutine cgedi(a,lda,n,ipvt,det,work,job)
c***begin prologue  cgedi
c***revision date  811015   (yymmdd)
c***category no.  f4f, f3
c***keywords  linpack,matrix,linear equations,complex,inverse,
c             determinant
c***date written  august 14, 1978
c
c***changed to double precision june 1985 by m.s.
c
c***author  moler c.b. (unm)
c***purpose
c  computes the determinant and inverse of a complex matrix
c   using the factors computed by cgeco or cgefa.
c***description
c     cgedi computes the determinant and inverse of a matrix
c     using the factors computed by cgeco or cgefa.
c
c     on entry
c
c        a       complex(lda, n)
c                the output from cgeco or cgefa.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c        ipvt    integer(n)
c                the pivot vector from cgeco or cgefa.
c
c        work    complex(n)
c                work vector.  contents destroyed.
c
c        job     integer
c                = 11   both determinant and inverse.
c                = 01   inverse only.
c                = 10   determinant only.
c
c     on return
c
c        a       inverse of original matrix if requested.
c                otherwise unchanged.
c
c        det     complex(2)
c                determinant of original matrix if requested.
c                otherwise not referenced.
c                determinant = det(1) * 10.0**det(2)
c                with  1.0 .le. cabs1(det(1)) .lt. 10.0
c                or  det(1) .eq. 0.0 .
c
c     error condition
c
c        a division by zero will occur if the input factor contains
c        a zero on the diagonal and the inverse is requested.
c        it will not occur if the subroutines are called correctly
c        and if cgeco has set rcond .gt. 0.0 or cgefa has set
c        info .eq. 0 .
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas caxpy,cscal,cswap
c     fortran abs,aimag,cmplx,mod,real
c
c***references
c  dongarra j.j., bunch j.r., moler c.b., stewart g.w.,
c   *linpack users  guide*, siam, 1979.
c***routines called  cswap,caxpy,cscal
c***end prologue  cgedi
      implicit real*8 (a-h,o-z)
      integer lda,n,ipvt(1),job
      complex*16 a(lda,1),det(2),work(1)
c
      complex*16 t
      real*8 ten
      integer i,j,k,kb,kp1,l,nm1
      complex*16 zdum
      real*8 cabs1
      cabs1(zdum) = abs(dble(zdum)) + abs(dimag(zdum))
c
c     compute determinant
c
c***first executable statement  cgedi
      if (job/10 .eq. 0) go to 70
	 det(1) = (1.0e0,0.0e0)
	 det(2) = (0.0e0,0.0e0)
	 ten = 10.0e0
	 do 50 i = 1, n
	    if (ipvt(i) .ne. i) det(1) = -det(1)
	    det(1) = a(i,i)*det(1)
c        ...exit
	    if (cabs1(det(1)) .eq. 0.0e0) go to 60
   10       if (cabs1(det(1)) .ge. 1.0e0) go to 20
	       det(1) = cmplx(ten,0.0e0)*det(1)
	       det(2) = det(2) - (1.0e0,0.0e0)
	    go to 10
   20       continue
   30       if (cabs1(det(1)) .lt. ten) go to 40
	       det(1) = det(1)/cmplx(ten,0.0e0)
	       det(2) = det(2) + (1.0e0,0.0e0)
	    go to 30
   40       continue
   50    continue
   60    continue
   70 continue
c
c     compute inverse(u)
c
      if (mod(job,10) .eq. 0) go to 150
	 do 100 k = 1, n
	    a(k,k) = (1.0e0,0.0e0)/a(k,k)
	    t = -a(k,k)
	    call cscal(k-1,t,a(1,k),1)
	    kp1 = k + 1
	    if (n .lt. kp1) go to 90
	    do 80 j = kp1, n
	       t = a(k,j)
	       a(k,j) = (0.0e0,0.0e0)
	       call caxpy(k,t,a(1,k),1,a(1,j),1)
   80       continue
   90       continue
  100    continue
c
c        form inverse(u)*inverse(l)
c
	 nm1 = n - 1
	 if (nm1 .lt. 1) go to 140
	 do 130 kb = 1, nm1
	    k = n - kb
	    kp1 = k + 1
	    do 110 i = kp1, n
	       work(i) = a(i,k)
	       a(i,k) = (0.0e0,0.0e0)
  110       continue
	    do 120 j = kp1, n
	       t = work(j)
	       call caxpy(n,t,a(1,j),1,a(1,k),1)
  120       continue
	    l = ipvt(k)
	    if (l .ne. k) call cswap(n,a(1,k),1,a(1,l),1)
  130    continue
  140    continue
  150 continue
      return
      end
      subroutine cgefa(a,lda,n,ipvt,info)
c***begin prologue  cgefa
c***revision date  811015   (yymmdd)
c***category no.  f4f
c***keywords  linpack,matrix,linear equations,complex,factor,
c             gauss elimination
c***date written  august 14, 1978
c
c***changed to double precision june 1985 by m.s.
c
c***author  moler c.b. (unm)
c***purpose
c   factors a complex matrix by gaussian elimination.
c***description
c     cgefa factors a complex matrix by gaussian elimination.
c
c     cgefa is usually called by cgeco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c     (time for cgeco) = (1 + 9/n)*(time for cgefa) .
c
c     on entry
c
c        a       complex(lda, n)
c                the matrix to be factored.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c     on return
c
c        a       an upper triangular matrix and the multipliers
c                which were used to obtain it.
c                the factorization can be written  a = l*u  where
c                l  is a product of permutation and unit lower
c                triangular matrices and  u  is upper triangular.
c
c        ipvt    integer(n)
c                an integer vector of pivot indices.
c
c        info    integer
c                = 0  normal value.
c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
c                     condition for this subroutine, but it does
c                     indicate that cgesl or cgedi will divide by zero
c                     if called.  use  rcond  in cgeco for a reliable
c                     indication of singularity.
c
c     linpack. this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas caxpy,cscal,icamax
c     fortran abs,aimag,real
c
c***references
c  dongarra j.j., bunch j.r., moler c.b., stewart g.w.,
c   *linpack users  guide*, siam, 1979.
c***routines called  caxpy,cscal,icamax
c***end prologue  cgefa
      implicit real*8 (a-h,o-z)
      integer lda,n,ipvt(1),info
      complex*16 a(lda,1)
c
      complex*16 t
      integer icamax,j,k,kp1,l,nm1
      complex*16 zdum
      real*8 cabs1
      cabs1(zdum) = abs(dble(zdum)) + abs(dimag(zdum))
c
c     gaussian elimination with partial pivoting
c
c***first executable statement  cgefa
      info = 0
      nm1 = n - 1
      if (nm1 .lt. 1) go to 70
      do 60 k = 1, nm1
	 kp1 = k + 1
c
c        find l = pivot index
c
	 l = icamax(n-k+1,a(k,k),1) + k - 1
	 ipvt(k) = l
c
c        zero pivot implies this column already triangularized
c
	 if (cabs1(a(l,k)) .eq. 0.0e0) go to 40
c
c           interchange if necessary
c
	    if (l .eq. k) go to 10
	       t = a(l,k)
	       a(l,k) = a(k,k)
	       a(k,k) = t
   10       continue
c
c           compute multipliers
c
	    t = -(1.0e0,0.0e0)/a(k,k)
	    call cscal(n-k,t,a(k+1,k),1)
c
c           row elimination with column indexing
c
	    do 30 j = kp1, n
	       t = a(l,j)
	       if (l .eq. k) go to 20
		  a(l,j) = a(k,j)
		  a(k,j) = t
   20          continue
	       call caxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
   30       continue
	 go to 50
   40    continue
	    info = k
   50    continue
   60 continue
   70 continue
      ipvt(n) = n
      if (cabs1(a(n,n)) .eq. 0.0e0) info = n
      return
      end
c/* @(#)cmatin.f        1.1      latest revision 7/30/86 14:16:36 */
c/* %w%      latest revision %g% %u% */
      subroutine cmatin(a,b,indexa,indexb,ipivot,n,ndim,detr,deti)
c
c
c     finds the inverse of complex*16 matrix h=(rea,aima)=(a,b)
c
c     code provided by norman bardsley to kwon+tabakin s prog bopit
c      n=no of elements in a assigned(and in the sub=dimension)
c      matrix a,b(real+imag parts) destroyed by call
c      the other arguments are dummy just to save space
      implicit real*8 (a-h,o-z)
      dimension  indexa(n),indexb(n),ipivot(n)
      dimension a(ndim,ndim),b(ndim,ndim)
      equivalence (irow,jrow),(icolum,jcolum)
c
      detr=1.0
c      rhl add to define tempi for real*8 potential
      tempi=0.
      deti=0.0
      do 10 j=1,n
  10  ipivot(j)=0
      do 70 i=1,n
      xmax=0.0
      do 35 j=1,n
      if(ipivot(j)-1) 15, 35, 15
  15  do 30 k=1,n
      if (ipivot(k)-1) 20, 30, 90
  20  if(abs(a(j,k)).lt.1.0e-20) go to 30
      xjk=abs(a(j,k))+abs(b(j,k))
      if(xmax-xjk)  25,30,30
  25  irow=j
      icolum=k
      xmax=xjk
  30  continue
  35  continue
      ipivot(icolum)=ipivot(icolum)+1
      if(irow.eq.icolum) go to 50
      detr=-deti
      deti=-deti
      do 45 l=1,n
      swap=a(irow,l)
      swapi=b(irow,l)
      a(irow,l)=a(icolum,l)
      b(irow,l)=b(icolum,l)
      b(icolum,l)=swapi
  45  a(icolum,l)=swap
  50  indexa(i)=irow
      indexb(i)=icolum
      pivotr=a(icolum,icolum)
      pivoti=b(icolum,icolum)
      temp=detr*pivotr-deti*pivoti
      deti=detr*pivoti+deti*pivotr
      detr=temp
      a(icolum,icolum)=1.0
      b(icolum,icolum)=0.0
      if(pivoti.eq.0.0) tempr=1.0/pivotr
      if(pivoti.ne.0.0) tempr=pivotr/(pivotr*pivotr+pivoti*pivoti)
      if(pivoti.ne.0.0) tempi=-pivoti/(pivotr*pivotr+pivoti*pivoti)
      do 55 l=1,n
      temp=a(icolum,l)*tempr-b(icolum,l)*tempi
      b(icolum,l)=a(icolum,l)*tempi+b(icolum,l)*tempr
  55  a(icolum,l)=temp
      do 70 l1=1,n
      if(l1-icolum) 60, 70, 60
  60  tempa=a(l1,icolum)
      tempb=b(l1,icolum)
      a(l1,icolum)=0.0
      b(l1,icolum)=0.0
      do 65 l=1,n
      b(l1,l)=b(l1,l)-a(icolum,l)*tempb-b(icolum,l)*tempa
      a(l1,l)=a(l1,l)-a(icolum,l)*tempa+b(icolum,l)*tempb
  65  continue
  70  continue
      do 85 i=1,n
      l=n+1-i
      if(indexa(l)-indexb(l)) 75, 85, 75
  75  jrow=indexa(l)
      jcolum=indexb(l)
      do 80 k=1,n
      swap=a(k,jrow)
      swapi=b(k,jrow)
      a(k,jrow)=a(k,jcolum)
      b(k,jrow)=b(k,jcolum)
      a(k,jcolum)=swap
      b(k,jcolum)=swapi
  80  continue
  85  continue
  90  return
      end
      subroutine cscal(n,ca,cx,incx)
c***begin prologue  cscal
c***revision date  811015   (yymmdd)
c***category no.  f1a, m2
c***keywords  complex,blas,vector,scale
c***date written  october 1979
c
c***changed to double precision june 1985 by m.s.
c
c***author lawson c. (jpl),hanson r. (sla),
c                            kincaid d. (u texas), krogh f. (jpl)
c***purpose
c    complex vector scale x = a*x
c***description
c                b l a s  subprogram
c    description of parameters
c
c     --input--
c        n  number of elements in input vector(s)
c       ca  complex scale factor
c       cx  complex vector with n elements
c     incx  storage spacing between elements of cx
c
c     --output--
c    cscal  complex result (unchanged if n.le.0)
c
c     replace complex cx by complex ca*cx.
c     for i = 0 to n-1, replace cx(1+i*incx) with  ca * cx(1+i*incx)
c
c
c***references
c  lawson c.l., hanson r.j., kincaid d.r., krogh f.t.,
c   *basic linear algebra subprograms for fortran usage*,
c  algorithm no. 539, transactions on mathematical software,
c  volume 5, number 3, september 1979, 308-323
c***routines called  (none)
c***end prologue  cscal
c
      implicit real*8 (a-h,o-z)
      complex*16 ca,cx(1)
c***first executable statement  cscal
      if(n .le. 0) return
      ns = n*incx
	  do 10 i = 1,ns,incx
	  cx(i) = ca*cx(i)
   10     continue
      return
      end
      subroutine cswap(n,cx,incx,cy,incy)
c***begin prologue  cswap
c***revision date  811015   (yymmdd)
c***category no.  f1a
c***keywords  complex,blas,vector,interchange
c***date written  october 1979
c
c***changed to double precision june 1985 by m.s.
c
c***author lawson c. (jpl),hanson r. (sla),
c                            kincaid d. (u texas), krogh f. (jpl)
c***purpose
c     interchange complex vectors
c***description
c                b l a s  subprogram
c    description of parameters
c
c     --input--
c        n  number of elements in input vector(s)
c       cx  complex vector with n elements
c     incx  storage spacing between elements of cx
c       cy  complex vector with n elements
c     incy  storage spacing between elements of cy
c
c     --output--
c       cx  input vector cy (unchanged if n.le.0)
c       cy  input vector cx (unchanged if n.le.0)
c
c     interchange complex cx and complex cy
c     for i = 0 to n-1, interchange  cx(lx+i*incx) and cy(ly+i*incy),
c     where lx = 1 if incx .gt. 0, else lx = (-incx)*n, and ly is
c     defined in a similar way using incy.
c
c
c***references
c  lawson c.l., hanson r.j., kincaid d.r., krogh f.t.,
c   *basic linear algebra subprograms for fortran usage*,
c  algorithm no. 539, transactions on mathematical software,
c  volume 5, number 3, september 1979, 308-323
c***routines called  (none)
c***end prologue  cswap
c
      implicit real*8 (a-h,o-z)
      complex*16 cx(1),cy(1),ctemp
c***first executable statement  cswap
      if(n .le. 0)return
      if(incx.eq.incy.and.incx.gt.0) go to 20
      kx = 1
      ky = 1
      if(incx.lt.0) kx = 1+(1-n)*incx
      if(incy.lt.0) ky = 1+(1-n)*incy
	  do 10 i = 1,n
	  ctemp = cx(kx)
	  cx(kx) = cy(ky)
	  cy(ky) = ctemp
	  kx = kx + incx
	  ky = ky + incy
   10 continue
      return
   20 continue
      ns = n*incx
	  do 30 i=1,ns,incx
	  ctemp = cx(i)
	  cx(i) = cy(i)
	  cy(i) = ctemp
   30     continue
      return
      end
c/* @(#)dcalc.f 1.14      latest revision 4/15/87 13:03:13 */
c/* %w%      latest revision %g% %u% */
      subroutine dcalc(iwave,zp,zgkap,zgp)
c     subroutine dcalc******************************
c     this sub calc the prin value part of denominator
      implicit real*8 (a-h,k,m,o-y),complex*16(z)
      real*8 dks( 48),dwts( 48),dk2( 48),dwt2( 48)
      dimension kap(5,2),gam(5,2),zgkap(2),zgk(2),zgp(2),zregp(2)
      data zpold/(0.d00,0.d00)/
      zp2= zp**2
      if (zp.eq.zpold)return
      if (zp.eq.0.)zp=1.
      zpold=zp
      call gpbara(iwave,1,zp,kap,gam2,zgp(1))
      call gpbara(iwave,2,zp,kap,gam2,zgp(2))
c     calc const part of integrnd
      zgkap(1)=0.
      zregp(1)=zp2*zgp(1)**2
      zgkap(2)=0.
      zregp(2)=zp2*zgp(2)**2
c     do 100 ni=1,2
c     zgkap(ni)=0.
c 100 zregp(ni)=zp2*zgp(ni)**2
c     calc integral -map have to use more pts
      npt=24
      npt1=npt
      npt2=npt
      npt3=npt
      kmid=zp
      kmax=1800.d00
      call gauss (npt1,012,0.0d00 ,kmid,dks,dwts)
      call gauss (npt2,012,kmid,kmax,dk2,dwt2)                                 c
      npts=npt1+npt2+npt3
      do 10 n=1,npts
      if( n .gt. npt1) go to 9
      ki=dks(n)
      wti=dwts(n)
      go to 11
    9 continue
c     integrtal from kmax to infinity, using remainder not as good
      if(n .eq. (1+npt1+npt2)) call gauss (npt3,042,45.d03,kmax,dk2,
     1 dwt2)
      if( n.gt.(npt1+npt2)) go to 8
      nn=n-npt1
      ki=dk2(nn)
      wti=dwt2(nn)
      go to 11
    8 continue
      nn=n-npt1-npt2
      ki=dk2(nn)
      wti=dwt2(nn)
   11 continue
      k2=ki*ki
      zki =ki
c     find g at grid pts
      call gpbara(iwave,1,zki,kap,gam2,zgk(1))
      call gpbara(iwave,2,zki,kap,gam2,zgk(2))
      do 5 ni=1,2
      zgkap(ni)=zgkap(ni)+wti*(k2*zgk(ni)**2-zregp(ni))/(zp2-k2)
 5    continue
   10 continue
      if(nold.eq.o)nold=1
      return
      end
c/* @(#)detcal.f        1.46      latest revision 4/30/87 16:27:34 */
      subroutine detcal(nrei,xeri,detri,iflag)
c       version modified sept85 for antiproton bs,  8/86 for reltv 
      implicit real*8  (a-h,o-y)
      integer mypvt(192)
      complex*16 detin(2),work(192)
      complex*16 h, x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,
     1 x13,x14,x15,x16,x17
      complex*16 zi,zw(10), za(192,192),zdet
     1 ,den(192,3),sum1,sum2,sum3,zk22,zk33,zk11,zk2,zk1,zk3
     2 ,ze,zs,ar11,ar22,ar33
      real*8 kapg,mn,mn2,k0,ko,ko2,k,k2,kk,mpi,mpi2,mr,k02
     1, imtij,bnuc,bn,retij,xbarc,xi,xpi,xn,xgam,zcm,zch,theta,
     1dsig,bnucf,bnf,dflip,dnof,wsp,wsn,achp,achn,acmp,acmn,ymin1,
     2ymin2
c      NB NEED MATCH DIMEN AND NO ENTRIES IN DATA STM FOR IBM
      dimension m1(100),m2(100),y(50),gin(14,100),h(15,20),kapg(100
     1),  x1(15),x2(1),x3(4),
     2 x4(3),x5(10),x6(8),x7(6),x8(6),x9(2),x10(9),x11(1),
     3 x12(1),x13(1),x14(1),x15(4),x16(4),x17(2),intger(192),mm1(192),
     4 mm2(192),xeri(2),detri(2)
       dimension kk(130),u(36864,2),f(36864),gp(66),wt(66)
     1,ptptw(66),wa(130)
      dimension theta(95), yy(50), ff(4), dflip(95), dnof(95)
      dimension nifty(20),bnuc(20),bn(16),retij(14),imtij(14),bnucf
     1(20), dsig(95),chart(61,61),polar(61),denom(14,25),ecmden(25),
     2 bnf(16),  
     3 reul(25,2), aimul(25,2)
C      NEW DIMENSIONS AND EQUIVL FOR BS,BE CAREFUL WITH U NEED
C       F CAN BE USED.nb order u(11) u(21)...u(n1) u(12) u(22) ..u(n2)
c
      common /nlspfl/ gin,kapg,denom,ecmden,fmn,rew(14,100),aimw(14,100)
     1,akapd(100),nkapgs,necms,dummy(1)/foptp/f
      common /sec2/ bnuc,bn,bnucf,bnf,retij,imtij,xbarc,xi,xpi,xn,nz,nes
     1,nwaves,nifty,na
      common /sec4/ achp,acmp,wsp,achn,acmn,wsn/sec6/rcoul,rcut
      common /absp/ b0r,b0i,c0r,c0i
      common/det/mr,alfanz
      equivalence(u(1,1),y(1)),(u(51,1),yy(1)),(m2(1),dsig(1)),(chart(1,
     11),u(1,1)), (gin(1,1),u(1,1)), (aimw(12,34),dnof(1)), (aimw(14,41)
     2,dflip(1)),(aimw(2,49),polar(1)),(gin(1,1),reul(1,1)),(gin(1,5),
     3aimul(1,1)),(x1(1),h(1,1)),(x2(1),h(1,2)),(x3(1),h(1,3)),
     4(x4(1),h(1,4)),(x5(1),h(1,5)),(x6(1),h(1,6)),(x7(1),h(1,7)),
     5(x8(1),h(1,8)),(x9(1),h(1,9)),(x10(1),h(1,10)),(x11(1),h(1,11)),
     6(x12(1),h(1,12)),(x13(1),h(1,13)),(x14(1),h(1,14)),(x15(1),h(1,15)
     7),(x16(1),h(1,16)),(x17(1),h(1,17)),(ko,k0),(ko2,k02),
     9(m1(1),f(800)),(m2(1),f(1000)),(intger(1),f(1200)),(wa(1),mm1(1))
     a,(m1(1),mm1(1)),(m2(1),mm2(1)),(za(1,1),u(1,1)),
     b(h(1,1),u(1,1)),(zw(1),f(1))
c-1     ibm version with (u,h) equivalence removed to avoid block data
c
c      revised constants and masses Oct. 2000 for bound s problem
      data alfa /7.29735d-3/,an,ap,akmn,ak0,asig,ndim/
     1 939.56533d+0,938.271998d+0,493.677d+0,497.672d+0,1193.154d+0,
     2 192/,ix/1/
      data hbarc,api,pi/197.3269602d+0,139.57018d+0,3.14159265359d+0/
     1 ,x1/
     1 10h(p,n)     ,10h(k+,ko)   ,10hpi+       ,10hpi-       ,
     2 10hpi+,pi-   ,10hpi cex    ,10hpio       ,10hproton    ,
     3 10hcex(+pio) ,10hk+        ,10hk0        ,10hneutron    ,
     3 10hk-        ,10hk0bar     ,10hp bar      /,  
     4x2     /10h*5=be         /,x3   /10hbs not res ,10hres not bs,
     5 10hna         ,10hna        /,x4    /10htpin=const,
     6 10hnlsp-g(p)  ,10hlocal lapp  /,x5    /10hso,kapo    ,
     7  10hsin,kap0   ,10hso,no ang  ,2*10hundef     ,10he3b,no aay
     8,10hundef      ,10he3b, aay   ,10hundefined ,10hso,aay    /,x6/
     910hno spin ms,10hss+ds      ,10hss+ds+ts   ,10hspin,ms    ,
     a 10hus=uc=0   ,10hss+unitary,10hus=0,ucne0,
     a 10huim = 0    /
      data x7/
     b 10hno abs/s-f,10hno abs/s-s,10hno abs/s-p  ,10hvabs-old  ,
     c10hvabs-yool  ,10habs-old,sp/,x8    /10hno del/fld,10hfld,no del ,
     d10hdel,no fld,10hdel in,fld,10hk+ martin ,10hmartin,fld /,x9/
     e 10hno pauli  ,10hpauli tpin /,x10/10hno coul    ,10h+fc+phase ,
     f10h+fc        ,10hexact coul ,10hexact+f(q) ,10hbs,pt coul,
     f 10hbs,sph col,10hbs,real cl,10hbs,no coul/,x11/
     g10hpi0,pi-,k+/,x12    /10hpi+,ko,k0b/,x13    /
     h 10hpi-channel  /,x14    /10hpi+channel/,x15    /10h1 chanel=p,
     h 10h2 chan=pip,10h3 chn=ppin,
     i10h2 chan=pn /,x16    /10hrhog,rho2g,10hrhows,rh2g,10hrhows,r2ws,
     j10hhe/d      /,x17    /10hkmt       ,10hwatson    /
      if(ix.ne.1) go to 150
 10   continue
C     MN HERE REFERS TO NUCLEUS MASS,XN (IN COMMON) TO NUCLEOn
C     READ  MOMENTA  IN  CENTRE _OF_MASS
C      NR.LE.0  IS  NONRELATIVISTIC  NR .GT. 0 >     RELATIVISTIC  CASe
      read(5,*,end=650) nr,lxmax
      write (6,740) nr,lxmax
      read (5,*) ngp,kode,b,nang,ymin1,ymin2
      if (nang.eq.0) nang = 3
      write (6,760) ngp,kode,b,nang,ymin1,ymin2
      if (ngp.le.0) go to 650
      write (6,790) ngp
C     -----READ IN DATA
      read (5,*) achp,acmp,wsp,achn,acmn,wsn,rcoul,rcut
      write (6,880) achp,acmp,wsp,achn,acmn,wsn,rcoul,rcut
      read (5,*) nz,na,(nifty(n),n=1,20)
      write (6,830) nz,na,(nifty(n),n=1,20)
      nifty(15)=abs(nifty(15))
C     PRINT OUT NIFTY
      do 25 nif=1,17
      jj=0
	 j = nifty(nif)+1
	 if (nif.eq.1) j = nifty(nif)+3
	 if (nif.eq.2.or.nif.eq.11.or.nif.eq.12.or.nif.eq.13.or.nif.eq.
     1   14 )jj=1
	 if (jj.eq.1) j=1
	 write (6,710) nif,nifty(nif),h(j,nif)
	 if (jj.eq.1) write(6,711)
 25   continue
      read (5,*) nes,nwaves,b0r,b0i,c0r,c0i
      write (6,910) nes,nwaves,b0r,b0i,c0r,c0i
      if (b0i.gt.1.d-03.or.c0i.gt.1.d-03) go to 30
c      use old values if nothing is read in
      b0r = -0.04d+00
      b0i = 0.04d+00
      c0r = 0.d+00
      c0i = 0.08d+00
      write (6,910) nes,nwaves,b0r,b0i,c0r,c0i
	if (nwaves.lt.0)write(6,884)nwaves
 884    format(' Nwaves lt 0 for % partial absorption, Nwaves=',i4)
 30   if (na.ne.3) go to 40
c      square input size params for he3(prevoius read in squaed)
      achp = achp*achp
      acmp = acmp*acmp
      wsp = wsp*wsp
      write(6,880)achp,acmp,wsp
 40   continue
      mpi=api
      if ((nifty(1).eq.7).or.(nifty(1).eq.8).or.(nifty(1).eq.-1)) mpi =
     1 akmn
      if(nifty(1).eq.10.or.nifty(1).eq.11)mpi=akmn
	if (nifty(1).eq.12)mpi=ap
      zcm = acmp
      amass=na
      zch = achp
      zws = wsp
      nifty1 = nifty(1)
c     for pi 0(ko) , use   no coul if not b.s.
      if(nifty(10).gt.4) go to 45
      if (nifty1.eq.4.or.nifty1.eq.8.or.nifty1.eq.11) nifty(10) = 0
C      FOR B S ALWAYS USE WATSON POTEMNTIAl
 45    if(nifty(10).ge.5)nifty(17)=1
c       
c       DEFINE MASSES WHICH DONT CHANGE WITH ENERGY
c
c         mn here refers to nucleus mass,xn is nucleon
      mpi2 = mpi*mpi
      zi=(0.d+00,1.d+00)
      xn = (nz*ap+(na-nz)*an)/amass
      mn=xn*amass
      mn2=mn*mn
      xbarc = hbarc
      xi = pi
      xpi = mpi
      mr = mn*mpi/(mn+mpi)
      ar11=mr
      ar112=ar11*ar11
      alfanz=nz*alfa
      ar2=mr*mr
c           need readjust mass for other nuclei
C            SPECIAL COUPLED CHANNELS SECTION
      e12=akmn+ap-api-asig
      e13 = akmn + ap -ak0 - an
	if (nifty(1).eq.12)then
	e12=0.
	e13=ap +ap -an -an
	ak0=an
	akmn=ap
	end if
C       PBAR-P CHANNELS
C       1--PBAR P    2-- NOTHING OPEN !!  3- NBAR N
c
c      for n.r. case use haw's reltivis  sig pi channel reduced mass
	am1=akmn+ap
	am2=asig+api
	am3=ak0+an
	am12=am1**2
	am22=am2**2
	am32= am3**2
	am1m=(akmn-ap)**2
	am2m=(asig-api)**2
	am3m =(ak0 -an )**2
	a21=asig**2
	a22=api**2
	a31=ak0**2
	a32=an**2
	ar22=(am1+am2)*(am12-am2m)/8.d+00/am12
c       set booby trap on ar22 for pbar
	if(nifty(1).eq.12)ar22=0.
      ar222=ar22*2.d+00
      api2=api*api
      ar33=ak0*an/(ak0+an)
      ar332=ar33*2.d+00
      ak02=ak0**2.d+00
      an2=an**2.d+00
      asig2=asig**2.d+00
 65   continue
       if(nr.le.0.) write(6,1000)
       if(nr.gt.0.) write(6,980)
      write (6,830) nz,na,(nifty(n),n=1,20)
      write (6,780) mpi,amass,mn
      write (6,781)ar11,ar22,ar33
      rhl2 =4.d+00*mr/pi
      rhl22=4.d+00*ar22/pi
      rhl23=4.d+00*ar33/pi
      rhl3 = 2.d00/pi
      rhl4 = 2.*rhl3
      wron1=2.d+00*mn
      wron2=2.d+00*mpi
      wron3=mn/10.d+00
      wron4=mpi/10.d+00
      ldmax=lxmax
      if(lxmax.eq.0)ldmax=1
      lborn=ldmax
      xgam=0.d+00
c     include coulomb, set uop
      zp = 1
      nzp = nifty(1)+1
      go to (150,130,140,140,140,150,140,150,140,140,130,140,130), nzp
 130  zp = -1
      go to 150
 140  zp = 0.
c-----------come here directly on all but 1st call-----------------
  150  continue
C      bs ENERGY IS ke IN com FOR BOTH rEL AND nR
c      need redefine for scatt where KE  is lab KE
      e=xeri(1)
      eimin =xeri(2)
      zdet=zi*eimin
	ze=e+zdet
      e2=e+e12
      e3=e+e13
	if(nr.gt.0)goto 156
c             calc complex*16 e s and momenta for chan 2
      zk11= (e+zdet)*2.*mr
      zk22=(e2+zdet)*ar222
      zk33=(e3+zdet)*ar332
	go to 157
c       RELATIVISTIC ENERGY
c       take sqrt(s) = e + m1 +m2 as def of E, same in all channels
c       zk11 is now k2 in channel 1, NOT energy, ar is complex mr
  156   ze = e + zdet +mn +mpi
	zs = ze**2
	zk11 = (zs-am12)*(zs-am1m)/4./zs
	zk22 = (zs-am22)*(zs-am2m)/4./zs
	zk33 = (zs-am32)*(zs-am3m)/4./zs
	ar11 = sqrt((zk11 + mn2)*(zk11+mpi2))/ze
	ar22 = sqrt((zk22 + a21)*(zk22+a22))/ze
	ar33 = sqrt((zk33 + a31)*(zk33+a32))/ze
	if(ix.eq.1)write(6,781) ar11,ar22,ar33
 157   continue
c     place complex momentum in quadrant II for BS analytic continuation
c     place complex momentum in quadrant IV for Resonance (nif(3)=1)
       x = zk11
       yi = -zi*zk11
       th = atan2(yi,x) + 2.*pi
       if(nifty(3).eq.1) th = atan2(yi,x) 
       zk1 = sqrt(sqrt(x**2 + yi**2))*(cos(th/2.) + zi*sin(th/2.))
       x = zk22
       yi = -zi*zk22
       th = atan2(yi,x) +2.*pi
       if(nifty(3).eq.1) th = atan2(yi,x) 
       zk2 = sqrt(sqrt(x**2 + yi**2))*(cos(th/2.) + zi*sin(th/2.))
       x = zk33
       yi = -zi*zk33
       th = atan2(yi,x) +2.*pi
       if(nifty(3).eq.1) th = atan2(yi,x) 
       zk3 = sqrt(sqrt(x**2 + yi**2))*(cos(th/2.) + zi*sin(th/2.))
       if(ix.eq.1)
     1  write(6,108) e,eimin,e2,eimin,e3,eimin,zk1,zk2,zk3
c        theta(1),(2) used in call to optpbs if whant complex momentum
	 theta(1)=zk2
	 theta(2)=-zi*zk2
c
C-----------SET UP GRID POINTS---------------------------------------
c
      n1 = ngp+1
      n2 = n1+n1
      a = zk2
      a=abs(a)
c       special check if b=0,use different distrb of points
      if(ix.gt.1)go to 89
      if (b.eq.0.) a = 200.
      if(kode .lt. 0) go to 84
      if(ngp.gt.48) go to 82
      call gauss (ngp,kode,a,b,gp,wt)
      go to 83
c      ngp gt 48,redistrb pts, 0-ko ,k0-inf(half up to 10*ko)
 82   npt1=24
      npt2=ngp-npt1
      b=a
      call gauss(npt1,011,0.d+00,b,gp,wt)
      if(npt2.gt.48)write(6,107)
      b=8.*a
      call gauss(npt2,041,b,a,dflip,dnof)
      ngp=npt1+npt2
      write(6,730)ngp,npt1,npt2
      do 81 i=1,npt2
      ii=i + npt1
      gp(ii)=dflip(i)
   81 wt(ii)=dnof(i)
      go to 83
 107  format(20h  NPT2 GT 48 IN MAIN  )
c      special bs point distb for kode lt 0-
 84   cpmax=10.d+00**(-kode)
      catom=b
      write(6,883) kode,catom
 883  format(' kode < 0, use bs points,cpmax=10**-kode, '
     1 ,'catom <,>0 :gausbs-rhl mod of kt',
     2 'gausbs-kt version ,kode,catom=',i4,e12.4)
      nnuc=nang
      if(nang.eq.3)nnuc=0
      cnucl=ymin2
      csize=ymin1
      write(6,104) ngp,nnuc,catom,csize,cnucl,cpmax
 104  format(40h Ngp,Nnuc,Catom,Csize,Cnucl,Cpmax(MeV)=   ,2i4,4e12.4/)
      if(catom.ge.0.)
     1call gausb2(ngp,nnuc,catom,csize,cnucl,cpmax,gp,wt,ptptw)
      if(catom.lt.0.)
     1call gausbs(ngp,nnuc,-catom,csize,cnucl,cpmax,gp,wt,ptptw)
 83   continue
      write(6,105)(gp(i),i=1,ngp)
      write(6,106)(wt(i),i=1,ngp)
 105  format(8h points=  ,10e12.4)
 106  format(8h weight=  ,10e12.4)
 89   continue
      sum1=0.
      sum2=0.
      sum3=0.
      if (nr.gt.0) go to 100
c      -----------------------------------------------------------------
c     NONRELATIVISTIC      DENOMINATOR
c      -----------------------------------------------------------------
      do 90 i1=1,ngp
	k = gp(i1)
	kk(i1) = k
	k2 = k*k
	ak2w=k2*wt(i1)
	den(i1,1)=ak2w/(k2-zk11)*rhl2
	if(nifty(15).eq.0) go to 90
	den(i1,2)=ak2w/(k2-zk22)*rhl22
	den(i1,3)=ak2w/(k2-zk33)*rhl23
c            2 closed channels, replace den2 with den3 and dont use channel 3
c                    closed channels: (K-p,Kobar-n) or (p-pbar,n-nbar)
	if(nifty(15).eq.3)den(i1,2)=den(i1,3)
	if(nifty(15).eq.3)den(i1,3)=0.
c            include Prin Value sum in denom for open channel
	if(e      .gt.0.)sum1 = sum1 +wt(i1)/(k2-zk11)
	if(e2.gt.0.)sum2 = sum2 +wt(i1)/(k2-zk22)
	if(e3.gt.0.)sum3 = sum3 +wt(i1)/(k2-zk33)
 90     continue
c      SET DEN(N1) WITH DEN2 CONTAINING BOTH I EPS AND p SUBTRN
c      can use this procedure for any real/complex energy, here use reE>0 switch
      if(e      .lt.0.)den(n1,1)=0.
      if(e2.lt.0.)den(n1,2)=0.
      if(e3.lt.0.)den(n1,3)=0.
      kk(n1)=zk1
      if(nifty(15).eq.0.and.e.lt.0.) go to 120
c     channel 2 (usually the only open one)
      if(e2.gt.0.)then
c        channel 2 open
	 kk(n1)=zk2
	 zw(1)= - sum2*rhl22*zk22
	 zw(2)=  + zi*ar222*zk2
	 den(n1,2)= zw(1) +zw(2)
	 endif
      if(e3.gt.0.)then
c         channel 3 open
	  zw(1)= - sum3*rhl23*zk33
	  zw(2)=  + zi*ar332*zk3
	  den(n1,3)= zw(1) +zw(2)
	  endif
      if(e.gt.0.)then
c         channel 1 open
	  zw(1)= - sum1*rhl2*zk11
	  zw(2)=  + zi*ar2*zk1
	  den(n1,1)= zw(1) +zw(2)
	  endif
c     for closed channels set den(n1,channel)=0
c     for 2 channel case (nbar-n or kobar-n) switch channel 3 to 2
      if(nifty(15).eq.3)then
		  den(n1,2)=den(n1,3)
		  den(n1,3)=0.
		  endif
      go to 120
c        ---------------------------------------------------------------
c                    RELATIVISTIC    DENOMINATOR
c        ---------------------------------------------------------------
 100  continue
	do 110 i1=1,ngp
	k = gp(i1)
	kk(i1) = k
	k2 = k*k
	ak2w=k2*wt(i1)
	emn = sqrt(k2+mn2)
	epi = sqrt(k2+mpi2)
c       special low momentum expansion, maybe do for other channels ?
	if(k.lt.wron3)emn=(k2/wron1)*(1.-k2/wron1/wron1)+mn
	if(k.lt.wron4)epi=(k2/wron2)*(1.-k2/wron2/wron2)+mpi
	den(i1,1)=ak2w* rhl3/(emn+epi-ze)
	if(nifty(15).eq.0)goto 110
	den(i1,2)=ak2w*rhl3/(sqrt(k2+a21)+sqrt(k2+a22)-ze)
	den(i1,3)=ak2w*rhl3/(sqrt(k2+a31)+sqrt(k2+a32)-ze)
c            2 closed channels, replace den2 with den3 and dont use channel 3
c            closed channels: (K-p,Kobar-n) or (p-pbar,n-nbar)
	if(nifty(15).eq.3)den(i1,2)=den(i1,3)
	if(nifty(15).eq.3)den(i1,3)=0.
c       include Prin Value sum in denom for open channel
	if(e      .gt.0.)sum1 = sum1 +wt(i1)/(k2-zk11)
	if(e2.gt.0.)sum2 = sum2 +wt(i1)/(k2-zk22)
	if(e3.gt.0.)sum3 = sum3 +wt(i1)/(k2-zk33)
 110    continue
c      set den(n1) with den2 containing both i eps and P subtrn
      if(e      .lt.0.)den(n1,1)=0.
      if(e2.lt.0.)den(n1,2)=0.
      if(e3.lt.0.)den(n1,3)=0.
      kk(n1)=zk1
      if(nifty(15).eq.0.and.e.lt.0.) go to 120
      if(e2.gt.0.)then
c        channel 2 open (the usual case for k-)
	 kk(n1)=zk2
	 den(n1,2)= (-sum2*rhl4*zk22 +zi*2.*zk2)*ar22
	 endif
      if(e3.gt.0.)then
c         channel 3 open
	  den(n1,3)= (-sum3*rhl4*zk33 +zi*2.*zk3)*ar33
	  endif
      if(e.gt.0.)then
c         channel 1 open
	  den(n1,1)= (-sum1*rhl4*zk11 +zi*2.*zk1)*ar11
	  endif
c     for closed channels set den(n1,channel)=0
c     for 2 channel case (nbar-n or kobar-n) switch channel 3 to 2
      if(nifty(15).eq.3)then
		  den(n1,2)=den(n1,3)
		  den(n1,3)=0.
		  endif
c---------------end relativistic denom-------------------------------
  120 if(ix.eq.1.or.e.gt.0..or.e2.gt.0..or.e3.gt.0.)
     1  write(6,91) (den(n1,i),i=1,3)
      if(ix.ne.1)go to 203
      ii=3
      if(nifty(15).eq.0)ii=1
      do 202 i=1,ii
  202 write(6,102) (den(i1,i),i1=1,n1)
  102 format(' Zden(3/1,i)=',10e12.4)
c
C----- DO LOOP OVER L FOR THIS E AND  STATE
c
 203     l = ldmax-1
	 if(ix.eq.1)write (6,1070) l
c
C      BOUND STATE OR EIGENENERGY CALCULATION ( A* RESONANCES )
c      calc bound state e's only for ldmax
c
       theta(3)=10.
      call optpbs(kk,n1,wt,ldmax,ldmax,lborn,theta,e,eimin)
c      for three coupled channels possible, need 3*maxngp as dimen
c      nunit=ngp or n1 is size of block matrix
      nunit=ngp
      if(nifty(15).eq.1.or.nifty(15).eq.2)nunit=ngp+1
      nn=nunit
      ng2=nunit*2
      if(nifty(15).eq.2)nn=3*nunit
      if(nifty(15).eq.1.or.nifty(15).eq.3)nn=ng2
c      WRITE OUT MATRICs
      if(ix.ne.1)go to 208
      write(6,229)
      i=nn-1
      if (nn .lt.6) then
	 do 9995 ii = 1,nn
 9995    write(6,227) ii,(za(ii,ij),ij=1,nn)
      else
	 do 217 ii=1,nn,i
 217     write(6,227) ii,(za(ii,ij),ij=1,nn)
      endif
      ii=nn
 228  format(1h0,'  1-G(complex*16)*V matrix before  determinant'
     1 ,' found')
  229 format(1h0,' complex*16 potential matrix before 1-GV formed')
 227  format(1h0,i3,8e15.6/200(3x,8e15.6/))
 208  continue
C       DETERMINE CHANNEL FOR DEN + ROW OF SUBMATRX
      irow=1
      do 219 ii=1,nn
      if(ii.gt.nunit)irow=2
      if(ii.gt.ng2)irow=3
      do 216 ij=1,nn
      jjj=ij
      if(ij.gt.nunit)jjj=ij-nunit
      if(ij.gt.ng2)jjj=ij-ng2
      icol=1
c       revised version with g1 and g2 switched
      if(ij.gt.nunit)icol=2
      if(ij.gt.ng2)icol=3
      za(ii,ij)=za(ii,ij)*den(jjj,icol)
2140  continue
      if(ii.eq.ij)za(ii,ii)=za(ii,ii)+1.
  216 continue
  219 continue
C----------DO LOOP OVER II,IJ ENDs
c      write out matrics
      if(ix.gt.1)go to 298
      write(6,228)
      i=nn-1
      if (nn .lt. 6)then
	   do 9996 ii=1,nn
 9996        write(6,227) ii,(za(ii,ij),ij=1,nn)
      else
	   do 297 ii=1,nn,i
 297         write(6,227) ii,(za(ii,ij),ij=1,nn)
      endif
      if(nifty(6).eq.6.or.nifty(6).eq.7)go to 298
 298  continue
 212   continue
C    SEARCH FOR DET=0 WITH LANDE SUBRTN IN VCOUL
      call cgefa(za,192,nn,mypvt,info)
      call cgedi(za,192,nn,mypvt,detin,work,10)
c     calculate deter
      zdet=(1.,0.)
      do 223 i=1,nn
      ipvt=wa(i)
      if(ipvt.ne.i)zdet=-zdet
      zdet=zdet*za(i,i)
 223  continue
      if(ix.eq.1)ix=2
      rpow = detin(2)
      zdet = detin(1)*(10.d00**rpow)
      detri(1)=zdet
      detri(2)=-zi*zdet
      write(6,224)zk1,e,eimin,zdet
  224 format(' zk1=',2(1pE11.3),' Er,Ei=',
     1 2(1pe18.8),'  zdet=',2(e15.3)/)
      return
C     _______________     FORMATS
 650  write(6,841)
      stop
c
  108 format(' zE1,  zE2,   zE3           =',3(2e12.4,4x)/
     1       ' zk1,  zk2,   zk3           =',3(2e12.4,4x))
   91 format(' den(n1,1-3) (P, Del parts) =',3(2e12.4,4x))
 710   format(7h Nifty(,i2,2h)=,i3,10x,a8,a2)
  711 format(1h+,36x,5hshift )
 730  format (4x,4i4)
 740  format (' Nr Lmax =',5i4)
 760  format (2i5,f10.2,i5,2f10.0) 
 780  format(1h ,3x,' Mproj,  At,   MA  =',
     1       5x,3(f8.3,11x))
 781  format(1h ,3x,' ar11,   ar22, ar33=',
     1       5x,3(2f8.3,3x)/)
 790  format (1h ,18hNo of grid points=,i4)
 830  format (1h ,2i5,1x,i2,i2,8i1,10i3)
 841  format(12h EOF IN MAIN ,i3)
 850  format (i3,f15.5)
 851  format(29h g(p)-unit4/d(e)-unit2 print ,/i3,f15.5)
 860  format (f10.4/4(6f13.6,/),4f13.6)
 861  format(1x ,150(f10.4/6f13.6/6f13.6/2f13.6/))
 870  format (e10.4/8e10.3/6e10.3)
 871  format(1x ,150(f10.4/8e10.3/6e10.3/))
 880  format (8f10.5)
 910  format (28h Nes,Nwaves,b0r,b0i,c0r,c0i=,2i5,4f10.5)
 940  format (1h0,10x,8h  Energy,10x,15h    Momentum   ,/)
 950  format (1h ,7x,f14.4,10x,f14.4/)
 960  format (1h ,5f14.4)
 980  format (1h0,31h      RELATIVISTIC  CALCULATION)
 1000 format (1h0,' NONRELATIVSTC CALC')
 1070 format (1h0,30h----- Angular   Momentum  =   ,i4,5h-----)
      end
      subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
c***begin prologue  dogleg
c
c***changed to double precision
c
c***refer to  snsqe,snsq
c     **********
c
c     subroutine dogleg
c
c     given an m by n matrix a, an n by n nonsingular diagonal
c     matrix d, an m-vector b, and a positive number delta, the
c     problem is to determine the convex combination x of the
c     gauss-newton and scaled gradient directions that minimizes
c     (a*x - b) in the least squares sense, subject to the
c     restriction that the euclidean norm of d*x be at most delta.
c
c     this subroutine completes the solution of the problem
c     if it is provided with the necessary information from the
c     qr factorization of a. that is, if a = q*r, where q has
c     orthogonal columns and r is an upper triangular matrix,
c     then dogleg expects the full upper triangle of r and
c     the first n components of (q transpose)*b.
c
c     the subroutine statement is
c
c       subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
c
c     where
c
c       n is a positive integer input variable set to the order of r.
c
c       r is an input array of length lr which must contain the upper
c         triangular matrix r stored by rows.
c
c       lr is a positive integer input variable not less than
c         (n*(n+1))/2.
c
c       diag is an input array of length n which must contain the
c         diagonal elements of the matrix d.
c
c       qtb is an input array of length n which must contain the first
c         n elements of the vector (q transpose)*b.
c
c       delta is a positive input variable which specifies an upper
c         bound on the euclidean norm of d*x.
c
c       x is an output array of length n which contains the desired
c         convex combination of the gauss-newton direction and the
c         scaled gradient direction.
c
c       wa1 and wa2 are work arrays of length n.
c
c     subprograms called
c
c       minpack-supplied ... r1mach,enorm
c
c       fortran-supplied ... abs,dmax1,dmin1,sqrt
c
c     minpack. version of july 1978.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
c
c     epsmch is the machine precision.
c
c***routines called  r1mach,enorm
c***end prologue  dogleg
      implicit real*8 (a-h,o-z)
      integer n,lr
      real*8 delta
      real*8 r(lr),diag(n),qtb(n),x(n),wa1(n),wa2(n)
      integer i,j,jj,jp1,k,l
      real*8 alpha,bnorm,epsmch,gnorm,one,qnorm,sgnorm,sum,temp,zero
      real*8 r1mach,enorm
      data one,zero /1.0d0,0.0d0/
c***first executable statement  dogleg
      epsmch = r1mach(4)
c
c     first, calculate the gauss-newton direction.
c
      jj = (n*(n + 1))/2 + 1
      do 50 k = 1, n
	 j = n - k + 1
	 jp1 = j + 1
	 jj = jj - k
	 l = jj + 1
	 sum = zero
	 if (n .lt. jp1) go to 20
	 do 10 i = jp1, n
	    sum = sum + r(l)*x(i)
	    l = l + 1
   10       continue
   20    continue
	 temp = r(jj)
	 if (temp .ne. zero) go to 40
	 l = j
	 do 30 i = 1, j
	    temp = dmax1(temp,abs(r(l)))
	    l = l + n - i
   30       continue
	 temp = epsmch*temp
	 if (temp .eq. zero) temp = epsmch
   40    continue
	 x(j) = (qtb(j) - sum)/temp
   50    continue
c
c     test whether the gauss-newton direction is acceptable.
c
      do 60 j = 1, n
	 wa1(j) = zero
	 wa2(j) = diag(j)*x(j)
   60    continue
      qnorm = enorm(n,wa2)
      if (qnorm .le. delta) go to 140
c
c     the gauss-newton direction is not acceptable.
c     next, calculate the scaled gradient direction.
c
      l = 1
      do 80 j = 1, n
	 temp = qtb(j)
	 do 70 i = j, n
	    wa1(i) = wa1(i) + r(l)*temp
	    l = l + 1
   70       continue
	 wa1(j) = wa1(j)/diag(j)
   80    continue
c
c     calculate the norm of the scaled gradient direction,
c     normalize, and rescale the gradient.
c
      gnorm = enorm(n,wa1)
      sgnorm = zero
      alpha = delta/qnorm
      if (gnorm .eq. zero) go to 120
      do 90 j = 1, n
	 wa1(j) = (wa1(j)/gnorm)/diag(j)
   90    continue
c
c     calculate the point along the scaled gradient
c     at which the quadratic is minimized.
c
      l = 1
      do 110 j = 1, n
	 sum = zero
	 do 100 i = j, n
	    sum = sum + r(l)*wa1(i)
	    l = l + 1
  100       continue
	 wa2(j) = sum
  110    continue
      temp = enorm(n,wa2)
      sgnorm = (gnorm/temp)/temp
c
c     test whether the scaled gradient direction is acceptable.
c
      alpha = zero
      if (sgnorm .ge. delta) go to 120
c
c     the scaled gradient direction is not acceptable.
c     finally, calculate the point along the dogleg
c     at which the quadratic is minimized.
c
      bnorm = enorm(n,qtb)
      temp = (bnorm/gnorm)*(bnorm/qnorm)*(sgnorm/delta)
      temp = temp - (delta/qnorm)*(sgnorm/delta)**2
     *       + sqrt((temp-(delta/qnorm))**2
     *              +(one-(delta/qnorm)**2)*(one-(sgnorm/delta)**2))
      alpha = ((delta/qnorm)*(one - (sgnorm/delta)**2))/temp
  120 continue
c
c     form appropriate convex combination of the gauss-newton
c     direction and the scaled gradient direction.
c
      temp = (one - alpha)*dmin1(sgnorm,delta)
      do 130 j = 1, n
	 x(j) = temp*wa1(j) + alpha*x(j)
  130    continue
  140 continue
      return
c
c     last card of subroutine dogleg.
c
      end

c/* @(#)eigenc.f        1.1      latest revision 7/30/86 14:19:46 */
c/* %w%      latest revision %g% %u% */
      subroutine eigenc(n,hr,hi,imag,iw,p0,q0,ia,ib,ip,uu,vv,ndim)
c
c    rhl slight revision of
c     inverse iteration method
c     code provided by norman bardsley to kwon+tabakin s code bopit
c      n=no pts used (and in sub dimension of matrix)
c      (p0,q0) = intial complex*16 energy guess
c      this prog does only 1 eigenvalue per call
c      imag = 0 no imag potent used     =1  both real*8 + imag used
c      iw = -1 for intermediate print out    = 0 else
c      uu = complex*16 wavefuntion in p space (eigenvector) vv=dummy
c      other variables are just dummy passed to save space
c      nb hr hi matrices are destroyed by call to eigenc, need redefine
c
      implicit real*8 (a-h,o-z)
      complex*16 uu,vv,zero,one,uniti,temp,e0,ee
      dimension  hr(ndim,ndim),hi(ndim,ndim),ia(n),ib(n),ip(n),
     1 uu(ndim),vv(ndim)
      data  ivec,iterx,conv/1,20,1.e-6/
      data zero,one,uniti/(0.0,0.0),(1.0,0.0),(0.0,1.0)/
      itern=0
      temp=zero
      e0=p0*one+q0*uniti*imag
c
c      rhl fix of problem with nused(n)=ndim
      if(n.gt.ndim)write(6,900)
      if(n.gt.ndim)stop
 900  format(' *****n gt ndim in eigenc **********')
c      set up intial values for eigenvectors and h-eguess matrix
      do 40 i=1,n
      hr(i,i)=hr(i,i)-p0
       uu(i)=zero
      if(imag.ne.0) hi(i,i)=hi(i,i)-q0
 40    continue
      do 39 j=1,n
      if(imag.eq.0)uu(j)=one
      if(imag.ne.0)uu(j)=(1.0,1.0)
 39   continue
	if(imag.ne.0.and.imag.ne.1)return
c
c     finds the inverse of h matrix.
c
      call cmatin(hr,hi,ia,ib,ip,n,ndim)
c
c      mult psi by inverse of h-eguess
  50  imax=0
      amax=0.0
      do 70 i=1,n
      vv(i)=zero
      do 60 j=1,n
  60  vv(i)=vv(i)+(hr(i,j)*one+hi(i,j)*imag*uniti)*uu(j)
      if (abs(vv(i)).lt.amax) go to 70
      imax=i
      amax=abs(vv(i))
  70  continue
c      renorm psi, determine new e = ee (complex*16 guess e=e0)
      do 80 i=1,n
  80  uu(i)=vv(i)/vv(imax)
c
      ee=e0+one/vv(imax)
      if(iw.eq.-1) write(6,105) itern,ee,uu(ivec)
      change=abs((ee-temp)/ee)
      if(itern.ge.5.and.change.lt.conv) go to 85
      if(itern.eq.iterx) write(6,103)
      if(itern.eq.iterx) go to 84
      temp=ee
      itern=itern+1
      go to 50
c
c85   write(6,520)
85    continue
c520   format( ' momentum space wave function')
c     write(6,530)
c530   format( ' mom. in inv. fermis',' real*8 part   ',' imagi
c  a  nary part')
c     do 500 ij=1,n
c     aq=uu(ij)
c     bq=aimag(uu(ij))
c     write(6,510) p(ij),aq,bq
c510   format(3e15.5)
 500   continue
  84  p0=ee
      q0=-uniti*ee
  103 format(/15h non-convergent)
  105 format(19h itern,ee,uu(ivec)=,i4,5x,4(1pe20.10))
      return
      end
      double precision function enorm(n,x)
c***begin prologue  enorm
c
c***changed to double precision
c
c***refer to  snsqe,snsq,snlse,snls
c     **********
c
c     function enorm
c
c     given an n-vector x, this function calculates the
c     euclidean norm of x.
c
c     the euclidean norm is computed by accumulating the sum of
c     squares in three different sums. the sums of squares for the
c     small and large components are scaled so that no overflows
c     occur. non-destructive underflows are permitted. underflows
c     and overflows do not occur in the computation of the unscaled
c     sum of squares for the intermediate components.
c     the definitions of small, intermediate and large components
c     depend on two constants, rdwarf and rgiant. the main
c     restrictions on these constants are that rdwarf**2 not
c     underflow and rgiant**2 not overflow. the constants
c     given here are suitable for every known computer.
c
c     the function statement is
c
c       real function enorm(n,x)
c
c     where
c
c       n is a positive integer input variable.
c
c       x is an input array of length n.
c
c     subprograms called
c
c       fortran-supplied ... abs,sqrt
c
c     minpack. version of october 1979.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
c***routines called  (none)
c***end prologue  enorm
      implicit real*8 (a-h,o-z)
      integer n
      real*8 x(n)
      integer i
      real*8 agiant,floatn,one,rdwarf,rgiant,s1,s2,s3,xabs,x1max,x3max,
     *     zero
      data one,zero,rdwarf,rgiant /1.0d0,0.0d0,3.834d-20,1.304d19/
c***first executable statement  enorm
      s1 = zero
      s2 = zero
      s3 = zero
      x1max = zero
      x3max = zero
      floatn = n
      agiant = rgiant/floatn
      do 90 i = 1, n
	 xabs = abs(x(i))
	 if (xabs .gt. rdwarf .and. xabs .lt. agiant) go to 70
	    if (xabs .le. rdwarf) go to 30
c
c              sum for large components.
c
	       if (xabs .le. x1max) go to 10
		  s1 = one + s1*(x1max/xabs)**2
		  x1max = xabs
		  go to 20
   10          continue
		  s1 = s1 + (xabs/x1max)**2
   20          continue
	       go to 60
   30       continue
c
c              sum for small components.
c
	       if (xabs .le. x3max) go to 40
		  s3 = one + s3*(x3max/xabs)**2
		  x3max = xabs
		  go to 50
   40          continue
		  if (xabs .ne. zero) s3 = s3 + (xabs/x3max)**2
   50          continue
   60       continue
	    go to 80
   70    continue
c
c           sum for intermediate components.
c
	    s2 = s2 + xabs**2
   80    continue
   90    continue
c
c     calculation of norm.
c
      if (s1 .eq. zero) go to 100
	 enorm = x1max*sqrt(s1+(s2/x1max)/x1max)
	 go to 130
  100 continue
	 if (s2 .eq. zero) go to 110
	    if (s2 .ge. x3max)
     *         enorm = sqrt(s2*(one+(x3max/s2)*(x3max*s3)))
	    if (s2 .lt. x3max)
     *         enorm = sqrt(x3max*((s2/x3max)+(x3max*s3)))
	    go to 120
  110    continue
	    enorm = x3max*sqrt(s3)
  120    continue
  130 continue
      return
c
c     last card of function enorm.
c
      end
      subroutine fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,
     *                  wa1,wa2)
c***begin prologue  fdjac1
c
c***changed to double precision
c
c***refer to snsqe,snsq
c
c     subroutine fdjac1
c
c     this subroutine computes a forward-difference approximation
c     to the n by n jacobian matrix associated with a specified
c     problem of n functions in n variables. if the jacobian has
c     a banded form, then function evaluations are saved by only
c     approximating the nonzero terms.
c
c     the subroutine statement is
c
c       subroutine fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,
c                         wa1,wa2)
c
c     where
c
c       fcn is the name of the user-supplied subroutine which
c         calculates the functions. fcn must be declared
c         in an external statement in the user calling
c         program, and should be written as follows.
c
c         subroutine fcn(n,x,fvec,iflag)
c         integer n,iflag
c         real x(n),fvec(n)
c         ----------
c         calculate the functions at x and
c         return this vector in fvec.
c         ----------
c         return
c         end
c
c         the value of iflag should not be changed by fcn unless
c         the user wants to terminate execution of fdjac1.
c         in this case set iflag to a negative integer.
c
c       n is a positive integer input variable set to the number
c         of functions and variables.
c
c       x is an input array of length n.
c
c       fvec is an input array of length n which must contain the
c         functions evaluated at x.
c
c       fjac is an output n by n array which contains the
c         approximation to the jacobian matrix evaluated at x.
c
c       ldfjac is a positive integer input variable not less than n
c         which specifies the leading dimension of the array fjac.
c
c       iflag is an integer variable which can be used to terminate
c         the execution of fdjac1. see description of fcn.
c
c       ml is a nonnegative integer input variable which specifies
c         the number of subdiagonals within the band of the
c         jacobian matrix. if the jacobian is not banded, set
c         ml to at least n - 1.
c
c       epsfcn is an input variable used in determining a suitable
c         step length for the forward-difference approximation. this
c         approximation assumes that the relative errors in the
c         functions are of the order of epsfcn. if epsfcn is less
c         than the machine precision, it is assumed that the relative
c         errors in the functions are of the order of the machine
c         precision.
c
c       mu is a nonnegative integer input variable which specifies
c         the number of superdiagonals within the band of the
c         jacobian matrix. if the jacobian is not banded, set
c         mu to at least n - 1.
c
c       wa1 and wa2 are work arrays of length n. if ml + mu + 1 is at
c         least n, then the jacobian is considered dense, and wa2 is
c         not referenced.
c
c     subprograms called
c
c       minpack-supplied ... r1mach
c
c       fortran-supplied ... abs,dmax1,sqrt
c
c     minpack. version of june 1979.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
c
c     epsmch is the machine precision.
c
c***routines called  r1mach
c***end prologue  fdjac1
      implicit real*8 (a-h,o-z)
      integer n,ldfjac,iflag,ml,mu
      real*8 epsfcn
      real*8 x(n),fvec(n),fjac(ldfjac,n),wa1(n),wa2(n)
      integer i,j,k,msum
      real*8 eps,epsmch,h,temp,zero
      real*8 r1mach
      data zero /0.0d0/
c***first executable statement  fdjac1
      epsmch = r1mach(4)
c
      eps = sqrt(dmax1(epsfcn,epsmch))
      msum = ml + mu + 1
      if (msum .lt. n) go to 40
c
c        computation of dense approximate jacobian.
c
	 do 20 j = 1, n
	    temp = x(j)
	    h = eps*abs(temp)
	    if (h .eq. zero) h = eps
	    x(j) = temp + h
	    call fcn(n,x,wa1,iflag)
	    if (iflag .lt. 0) go to 30
	    x(j) = temp
	    do 10 i = 1, n
	       fjac(i,j) = (wa1(i) - fvec(i))/h
   10          continue
   20       continue
   30    continue
	 go to 110
   40 continue
c
c        computation of banded approximate jacobian.
c
	 do 90 k = 1, msum
	    do 60 j = k, n, msum
	       wa2(j) = x(j)
	       h = eps*abs(wa2(j))
	       if (h .eq. zero) h = eps
	       x(j) = wa2(j) + h
   60          continue
	    call fcn(n,x,wa1,iflag)
	    if (iflag .lt. 0) go to 100
	    do 80 j = k, n, msum
	       x(j) = wa2(j)
	       h = eps*abs(wa2(j))
	       if (h .eq. zero) h = eps
	       do 70 i = 1, n
		  fjac(i,j) = zero
		  if (i .ge. j - mu .and. i .le. j + ml)
     *               fjac(i,j) = (wa1(i) - fvec(i))/h
   70             continue
   80          continue
   90       continue
  100    continue
  110 continue
      return
c
c     last card of subroutine fdjac1.
c
      end
      subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)
c***begin prologue  fdjac2
c***refer to snlse,snls
c     **********
c
c     subroutine fdjac2
c
c     this subroutine computes a forward-difference approximation
c     to the m by n jacobian matrix associated with a specified
c     problem of m functions in n variables.
c
c     the subroutine statement is
c
c       subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)
c
c     where
c
c       fcn is the name of the user-supplied subroutine which
c         calculates the functions. fcn must be declared
c         in an external statement in the user calling
c         program, and should be written as follows.
c
c         subroutine fcn(m,n,x,fvec,iflag)
c         integer m,n,iflag
c         double precision x(n),fvec(m)
c         ----------
c         calculate the functions at x and
c         return this vector in fvec.
c         ----------
c         return
c         end
c
c         the value of iflag should not be changed by fcn unless
c         the user wants to terminate execution of fdjac2.
c         in this case set iflag to a negative integer.
c
c       m is a positive integer input variable set to the number
c         of functions.
c
c       n is a positive integer input variable set to the number
c         of variables. n must not exceed m.
c
c       x is an input array of length n.
c
c       fvec is an input array of length m which must contain the
c         functions evaluated at x.
c
c       fjac is an output m by n array which contains the
c         approximation to the jacobian matrix evaluated at x.
c
c       ldfjac is a positive integer input variable not less than m
c         which specifies the leading dimension of the array fjac.
c
c       iflag is an integer variable which can be used to terminate
c         the execution of fdjac2. see description of fcn.
c
c       epsfcn is an input variable used in determining a suitable
c         step length for the forward-difference approximation. this
c         approximation assumes that the relative errors in the
c         functions are of the order of epsfcn. if epsfcn is less
c         than the machine precision, it is assumed that the relative
c         errors in the functions are of the order of the machine
c         precision.
c
c       wa is a work array of length m.
c
c     subprograms called
c
c       user-supplied ...... fcn
c
c       minpack-supplied ... r1mach
c
c       fortran-supplied ... abs,dmax1,sqrt
c
c     minpack. version of june 1979.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
c
c     epsmch is the machine precision.
c
c***routines called  r1mach
c***end prologue  fdjac2
      integer m,n,ldfjac,iflag
      double precision epsfcn
      double precision x(n),fvec(m),fjac(ldfjac,n),wa(m)
      integer i,j
      double precision eps,epsmch,h,temp,zero
      double precision r1mach
      data zero /0.0d0/
c***first executable statement  fdjac2
      epsmch = r1mach(4)
c
      eps = sqrt(dmax1(epsfcn,epsmch))
      do 20 j = 1, n
	 temp = x(j)
	 h = eps*abs(temp)
	 if (h .eq. zero) h = eps
	 x(j) = temp + h
	 call fcn(m,n,x,wa,iflag)
	 if (iflag .lt. 0) go to 30
	 x(j) = temp
	 do 10 i = 1, m
	    fjac(i,j) = (wa(i) - fvec(i))/h
   10       continue
   20    continue
   30 continue
      return
c
c     last card of subroutine fdjac2.
c
      end
      subroutine fdump
c***begin prologue  fdump
c
c***changed to double precision
c
c***revision date  811015   (yymmdd)
c***category no.  q2,n2
c***keywords  symbolic dump,dump
c***date written  8/79
c***author  jones r.e. (sla)
c***purpose
c  symbolic dump (should be locally written)
c***description
c        ***note*** machine dependent routine
c        fdump is intended to be replaced by a locally written
c        version which produces a symbolic dump.  failing this,
c        it should be replaced by a version which prints the
c        subprogram nesting list.  note that this dump must be
c        printed on each of up to five files, as indicated by the
c        xgetua routine.  see xsetua and xgetua for details.
c
c     written by ron jones, with slatec common math library subcommittee
c     latest revision ---  23 may 1979
c
c***references
c   jones r.e., *slatec common mathematical library error handling
c    package*, sand78-1189, sandia laboratories, 1978.
c***routines called  (none)
c***end prologue  fdump
c
c***first executable statement  fdump
      return
      end
c/* @(#)gausb2.f        1.1      latest revision 7/30/86 14:20:03 */
c/* %w%      latest revision %g% %u% */
      subroutine gausb2(npt,nnuc,catom,csize,cnucl,cpmax,p,wp,ppw)
c     rhl mod of /tabakin + kwon  trig scaled grid
c     bound state problem
      implicit real*8 (a-h,o-z)
      dimension p(66),wp(66),ppw(66),x(64),wx(64)
      if(nnuc.ne.0) rata=catom/csize
      if(nnuc.eq.0) rata=0.
c     rata = catom/csize
      natom=npt-nnuc
c     atomic region
      if (natom.eq.0) go to 12
      call gauss(natom,010,-1.d+00,1.d+00,x,wx)
      do 10 i=1,natom
      z1=x(i)
      p(i)=catom*(1.+z1)/(1.-(1.-2.*rata)*z1)
      dpdx=2.*catom*(1.-rata)/((1.-(1.-2.*rata)
     1 *z1)**2)
      wp(i)=dpdx*wx(i)
 10   ppw(i)=p(i)*p(i)*wp(i)
      if(nnuc.eq.0) go to 50
 12   continue
c     nuclearr region grid
      nuc=nnuc/2
      if(nuc.ne.natom) call gauss(nuc,010,-1.d+00,1.d+00,x,wx)
      n1=natom+nuc
      n2=natom+1
      do 20 j=n2,n1
      p(j)=(csize+cnucl+(cnucl-csize)*x(j-natom))/2.
      dpdx=(cnucl-csize)/2.
      wp(j)=dpdx*wx(j-natom)
      ppw(j)=p(j)*p(j)*wp(j)
 20   continue
      n3=natom+2*nuc
      if(n3.ne.npt) write(6,900)npt,nnuc,nuc
 900  format(' nnuc not div by 2 in gausb2,npt,nnuc,nuc=',3i3)
      do 30 j=n1+1,n3
      p(j)=(cnucl+cpmax+(cpmax-cnucl)*x(j-n1))/2.
      dpdp=(cpmax-cnucl)/2.
      wp(j)=dpdp*wx(j-n1)
      ppw(j)=p(j)*p(j)*wp(j)
 30   continue
 50   continue
      return
      end
c/* @(#)gausbs.f        1.1      latest revision 7/30/86 14:20:15 */
c/* %w%      latest revision %g% %u% */
      subroutine gausbs(npt,nnuc,catom,csize,cnucl,cpmax,p,wp,ppw)
c     tabakin + kwon  trig scaled grid
c     bound state problem
      implicit real*8 (a-h,o-z)
      dimension p(66),wp(66),ppw(66),x(64),wx(64)
      pid4=3.14159265/4.
      if(nnuc.ne.0) rata=catom/csize
      if(nnuc.eq.0) rata=0.
c     rata = catom/csize
      ratb=cnucl/cpmax
      natom=npt-nnuc
c     atomic region
      if (natom.eq.0) go to 12
      call gauss(natom,010,-1.d+00,1.d+00,x,wx)
      do 10 i=1,natom
      z1=pid4*(1.+x(i))
      p(i)=catom*sin(z1)/(cos(z1) + rata*sin(z1))
      dpdx=catom*pid4/(cos(z1)+rata*sin(z1))**2
      wp(i)= dpdx*wx(i)
 10   ppw(i)=p(i)*p(i)*wp(i)
      if(nnuc.eq.0) go to 50
 12   continue
c     nuclearr region grid
      if(nnuc.ne.natom) call gauss(nnuc,010,-1.d+00,1.d+00,x,wx)
      do 20 j=natom+1,npt
      z2=pid4*(1.+x(j-natom))
      p(j)=csize+cnucl*sin(z2)/(cos(z2)+ratb*sin(z2))
      dpdx=cnucl*pid4/(cos(z2)+ratb*sin(z2))**2
      wp(j)=dpdx*wx(j-natom)
      ppw(j)=p(j)*p(j)*wp(j)
 20   continue
 50   continue
      return
      end
c/* @(#)gauss.f 1.1      latest revision 7/30/86 14:20:30 */
c/* %w%      latest revision %g% %u% */
*      subroutine gauss(npts, kode, a, b, xs, wts)
c     head.h Version 1.1
c***********************************************************************
c
c       Oregon State University Nuclear Theory Group
c       Multi-Channel Anti Kaon -- Nucleon Code
c       Guangliang He & Rubin H. Landau
c
c***********************************************************************
c     fortran.h Version 1.4
c***********************************************************************
c     gauss.F Version 1.3
c     Latest Revision Date 10:35:54 8/16/91
c***********************************************************************
c
      subroutine gauss(npts,job,a,b,x,w)
c
c***********************************************************************
c     rescale rescals the gauss-legendre grid points and weights
c
c     npts     number of points
c     job = 0  rescalling uniformly between (a,b)
c           1  for integral (0,b) with 50% points inside (0, ab/(a+b))
c           2  for integral (a,inf) with 50% inside (a,b+2a), with
c              rational scaling
c     job = 3  for integral (a, inf) with 50% inside (a, b+2a), with
c              atan() scaling. (a+b) has to great than zero.
c
c     x, w     output grid points and weights.
c***********************************************************************
c
      integer npts,job,m,i,j,kode
      double precision x(npts),w(npts),a,b,xi
      double precision t,t1,pp,p1,p2,p3,aj
      double precision eps,pi,zero,two,one,half,quarter,pio4
      parameter (pi = 3.14159265358979323846264338328, eps = 3.0d-15)
      parameter (zero=0.0d0,one=1.0d0,two=2.0d0)
      parameter (half=0.5d0,quarter=0.25d0)
c
c *** FIRST EXECTUABLE *************************************************
c
      kode=job
      job=0
      m=(npts+1)/2
      do 10020 i=1,m
	 t=cos(pi*(i-quarter)/(npts+half))
10000 continue
	 p1=one
	 p2=zero
	 aj=zero
	 do 10010 j=1,npts
	    p3=p2
	    p2=p1
	    aj=aj+one
	    p1=((two*aj-one)*t*p2-(aj-one)*p3)/aj
10010    continue
	 pp=npts*(t*p1-p2)/(t*t-one)
	 t1=t
	 t=t1-p1/pp
c
	 if(abs(t-t1).gt.eps) goto 10000
c
	 x(i)=-t
	 x(npts+1-i)=t
	 w(i)=two/((one-t*t)*pp*pp)
	 w(npts+1-i)=w(i)
10020 continue
c
c***********************************************************************
c     rescale the grid points
c***********************************************************************
c
      if (job.eq.0) then
c     scale to (a,b) uniformly
	 do 10030 i=1,npts
	    x(i)=x(i)*(b-a)/two+(b+a)/two
	    w(i)=w(i)*(b-a)/two
10030    continue
c
      elseif (job.eq.1) then
c     scale to (0,b) with 50% points inside (0,ab/(a+b))
	 do 10040 i=1,npts
	    xi=x(i)
	    x(i)=a*b*(one+xi)/(b+a-(b-a)*xi)
	    w(i)=w(i)*two*a*b*b/((b+a-(b-a)*xi)*(b+a-(b-a)*xi))
10040    continue
c
      elseif (job.eq.2) then
c     scale to (a,inf) with 50% points inside (a,b+2a)
	 do 10050 i=1,npts
	    xi=x(i)
	    x(i)=(b*xi+b+a+a)/(one-xi)
	    w(i)=w(i)*two*(a+b)/((one-xi)*(one-xi))
10050    continue
      else if (job .eq. 3) then
c     scale to (a, inf) with 50% points inside (a, b+2a)
	 pio4 = atan(one)
	 if (a+b .le. zero) stop 'gauss: a+b less than zero'
	 do 10060 i = 1, npts
	    xi = x(i)
	    x(i) = a+(a+b)*tan(pio4*(one+xi))
	    w(i) = w(i)*(a+b)*pio4
     &           /(cos(pio4*(one+xi))*cos(pio4*(one+xi)))
10060    continue
      else
	 stop 'gauss: Wrong value of job'
      endif
      job=kode
c
      return
      end

c/* @(#)gpbara.f        1.7      latest revision 4/15/87 12:59:41 */
c/* %w%      latest revision %g% %u% */
	   subroutine gpbara(iwave,ispin1,zp,akap,gam2,zgofp)
c ***********************************************************************
c
c       the annihilation form factor for the graz nbar-n potentials
c       the p|p and n|n channels are truly coupled, yet no interaction
c       within the annihilation channel. the model is therefore a coupled 
c       channels optical potential.
c       s and p waves included, coupled channels-spin not 
c
c       iwave code      1- 1s0  2- 3s1  3- 1p1  4- 3p0  5- 3p1
c       the l and j refer to intial nnbar channels
c       in boson-boson channel have power p**jnn and bb state
c                       1- 1s0  2- 1p1  3- 1p1  4- 1s0  5- 1p1
c
c       storage         kap(iwave,ispin)        gam(iwave,ispin)        
c
c                       u(11) u(21)...(n1) u(12) (22)....(n2).....
c
	   implicit real*8 (a-h,m,k,o-x), complex*16 (z)
	   dimension kap(5,2), gam(5,2)
	dimension nifty(20),bnuc(20),bnucf(20),bn(16),bnf(16),retij(14),
     1  aimtij(14)
	 common/sec2/bnuc,bn,bnucf,bnf,retij,aimtij,hbarc,pi,mpi,mn,nz,
     1  nes,nwaves,nifty,na
	   data kap/2953.,44744.,1148461.,2930.,120271., 
     1           3151.,40452.,891622.,598.,174495./,
     2       gam/2.6697,3.0032,4.6462,1.2277,3.2559,
     3           2.6697,3.9662,4.6462,1.2277,3.2559/,nint/0/
	   if (nint.ne.0) go to 100
c
c       intialize, convert to mev, square gamma
c
	   nint=1
	   write(6,900) (kap(iwave,ispin),gam(iwave,ispin),ispin=1,2)
 900    format(' in g(p)-anninilation  k,gam=',100(5f20.8/))
 901    format(' in g(p)-anninilation  k,gam=',100(5d20.8/))
	   do 10 ispin=1,2
	   do 10 isett=1,5
	   gam(isett,ispin)=(gam(isett,ispin)*hbarc)**2
	   n=1
	   if(isett.eq.2.or.isett.eq.4)n=2
	   if(isett.eq.3.or.isett.eq.5)n=3
c     use nwaves as percent cutoff on annihilation
c
 10     kap(isett,ispin)=  kap(isett,ispin)*(hbarc**n)
	   write(6,901) (kap(iwave,ispin),gam(iwave,ispin),ispin=1,2)
c
c       intialization over
c
 100    zp2=zp*zp
	   zgofp=1./(zp2+gam(iwave,ispin1))
	   if(iwave.lt.3)go to 110
c       pwave
	   zgofp=zp*(zgofp**2)
 110    akap=kap(iwave,ispin1)
	   gam2=gam(iwave,ispin1)
	   return
	   end
c/* @(#)gpbarp.f        1.1      latest revision 7/30/86 14:21:15 */
c/* %w%      latest revision %g% %u% */
	   subroutine gpbarp(ispin1,p,alamb,gofp)
c ***********************************************************************
c
c       the separable potential form factor for the graz nbar-n potentials
c       the p|p and n|n channels are truly coupled, yet no interaction
c       within the annihilation channel. the model is therefore a coupled 
c       channels optical potential.
c       s and p waves included, coupled channels-spin not 
c       nes in /sec2/ is iset, the switch for spin-isospin state
c       ispin1 = isospin(0or1)+1        n=(1,5)=term # in form factor
c
c       iset code       1- 1s0  2- 3s1  3- 1p1  4- 3p0  5- 3p1
c
c       storage         al(iset,ispin)  b(n,iset,ispin) c(n,iset,ispin)
c
c                       u(11) u(21)...(n1) u(12) (22)....(n2).....
c
	   implicit real*8 (a-h,m,k,o-z)
	   dimension al(5,2), b(5,5,2), c(5,5,2)
	dimension nifty(20),bnuc(20),bnucf(20),bn(16),bnf(16),retij(14),
     1  aimtij(14)
	common/sec2/bnuc,bn,bnucf,bnf,retij,aimtij,hbarc,pi,mpi,mn,nz,
     1  nes,nwaves,nifty,na
	   equivalence (iset,nes)
	   data al/-1.,-1.,-1.,-1.,+1.,-1.,-1.,-1.,+1.,-1./,
     1  b/1.,1.866066,2.6878754,3.4822022,4.2566996, 1.,2.,3.,4.,5.,
     2    1.2,2.2392792,3.2254504,4.1786426,5.1080396,1.1,2.2,3.3,
     3    4.4,5.5,1.4,3.0009657,4.6877173,6.4327107,8.2223325,
     4    1.,1.866066,2.6878754,3.4822022,4.2566996,1.2,2.4,3.6,4.8,6.
     5    , 1.2,2.2392792,3.2254504,4.1786426,5.1080396, 0.8,1.6,2.4,
     6    3.2,4.,1.,1.866066,2.6878754,3.4822022,4.2566996/
       data c
     a /-16.947223,19.613007,704.60953,-1620.0131,929.20364,7.715196,
     1 -308.80525,1745.465,-2823.1814,1366.6833, 2.2254985,-188.92457,
     2 2764.8436,-12098.478,12273.396, 3.975071,-413.1827,-125.86102,
     3 9082.9976,-12059.309, 47.040483,377.98701,-7419.3042,19707.559,
     4 -13200.458, 9.4015746,-354.45965,2489.8515,-4728.6714,2647.6517,
     5 -2.0706683,-29.268615,987.27114,-2414.2214,1484.4254,-1.2249177,
     6 116.21112,-822.3519,-1397.174,3540.5651, -0.6954015,110.12122,
     7-613.65913,1716.1131,-1535.9172, 0.53770179,-44.214116,876.85139,
     8 -4972.9828,5427.2397/,nint/0/
	   if (nint .ne. 0)go to 200
c
c       intialize, convert to mev, square beta(called b)
c       convert l to rhl conventn.. vrhl=vgraz
c
	   write (6,903) iset
	   write (6,900)  (al(iset,iso),iso=1,2)
	   write (6,901)  ((b(n,iset,iso),n=1,5),iso=1,2)
	   write (6,902)  ((c(n,iset,iso),n=1,5),iso=1,2)
 900    format(' al',/100(5f20.8/))
 901    format(' b ',/100(5f20.8/))
 902    format(' c ',/100(5f20.8/))
 903    format('0in vpbarp, nes=',i3,
     1  ': 1-1s0  2-3s1  3-1p1  4-3p0  5-3p1')
	   do 100 iso=1,2
	   do 100 isett=1,5
	   if(isett.le.2) al(isett,iso)=al(isett,iso)*hbarc*pi/2.
	   if (isett.gt.2) al(isett,iso)= al(isett,iso)*(hbarc**3)*pi/2.
	   do 100 n=1,5
 100    b(n,isett,iso)= (b(n,isett,iso)*hbarc)**2
c       write(6,903)
	   write (6,900) (al(iset,iso),iso=1,2)
	   write (6,901) ((b(n,iset,iso),n=1,5),iso=1,2)
	   write (6,902) ((c(n,iset,iso),n=1,5),iso=1,2)
	   nint=1
c
c       intialization over
c
 200    gofp=0.
	   p2=p*p
c       sum over 5 terms in form factor
	   do 250 n=1,5
	   if (iset.lt.3)term=c(n,iset,ispin1)/(p2+b(n,iset,ispin1))
	if (iset.ge.3)term=c(n,iset,ispin1)*p/((p2+b(n,iset,ispin1))**2)
 250    gofp=gofp+term
	   alamb=al(iset,ispin1)
	   return
	   end
      integer function i1mach(i)
c***begin prologue  i1mach
c***revision date  850615   (yymmdd)
c
c***changed to double precision
c
c***category no.  q
c***keywords  machine constants,integer
c***date written   1975
c***author fox p.a.,hall a.d.,schryer n.l. (bell labs)
c***purpose
c returns integer machine dependent constants
c***description
c * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c   these machine constant routines must be activated for
c   a particular environment.
c * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     i1mach can be used to obtain machine-dependent parameters
c     for the local machine environment.  it is a function
c     subroutine with one (input) argument, and can be called
c     as follows, for example
c
c          k = i1mach(i)
c
c     where i=1,...,16.  the (output) value of k above is
c     determined by the (input) value of i.  the results for
c     various values of i are discussed below.
c
c  i/o unit numbers.
c    i1mach( 1) = the standard input unit.
c    i1mach( 2) = the standard output unit.
c    i1mach( 3) = the standard punch unit.
c    i1mach( 4) = the standard error message unit.
c
c  words.
c    i1mach( 5) = the number of bits per integer storage unit.
c    i1mach( 6) = the number of characters per integer storage unit.
c
c  integers.
c    assume integers are represented in the s-digit, base-a form
c
c               sign ( x(s-1)*a**(s-1) + ... + x(1)*a + x(0) )
c
c               where 0 .le. x(i) .lt. a for i=0,...,s-1.
c    i1mach( 7) = a, the base.
c    i1mach( 8) = s, the number of base-a digits.
c    i1mach( 9) = a**s - 1, the largest magnitude.
c
c  floating-point numbers.
c    assume floating-point numbers are represented in the t-digit,
c    base-b form
c               sign (b**e)*( (x(1)/b) + ... + (x(t)/b**t) )
c
c               where 0 .le. x(i) .lt. b for i=1,...,t,
c               0 .lt. x(1), and emin .le. e .le. emax.
c    i1mach(10) = b, the base.
c
c  single-precision
c    i1mach(11) = t, the number of base-b digits.
c    i1mach(12) = emin, the smallest exponent e.
c    i1mach(13) = emax, the largest exponent e.
c
c  double-precision
c    i1mach(14) = t, the number of base-b digits.
c    i1mach(15) = emin, the smallest exponent e.
c    i1mach(16) = emax, the largest exponent e.
c
c  to alter this function for a particular environment,
c  the desired set of data statements should be activated by
c  removing the c from column 1.  also, the values of
c  i1mach(1) - i1mach(4) should be checked for consistency
c  with the local operating system.
c
c***references
c  fox p.a., hall a.d., schryer n.l.,*framework for a portable library*,
c  acm transaction on mathematical software, vol. 4, no. 2,
c  june 1978, pp. 177-188.
c***routines called  (none)
c***end prologue  i1mach
c
      implicit real*8 (a-h,o-z)
      integer imach(16),output
c
c     equivalence (imach(4),output)
c
c     machine constants for the burroughs 1700 system.
c
c     data imach( 1) /    7 /
c     data imach( 2) /    2 /
c     data imach( 3) /    2 /
c     data imach( 4) /    2 /
c     data imach( 5) /   36 /
c     data imach( 6) /    4 /
c     data imach( 7) /    2 /
c     data imach( 8) /   33 /
c     data imach( 9) / z1ffffffff /
c     data imach(10) /    2 /
c     data imach(11) /   24 /
c     data imach(12) / -256 /
c     data imach(13) /  255 /
c     data imach(14) /   60 /
c     data imach(15) / -256 /
c     data imach(16) /  255 /
c
c     machine constants for the burroughs 5700 system.
c
c     data imach( 1) /   5 /
c     data imach( 2) /   6 /
c     data imach( 3) /   7 /
c     data imach( 4) /   6 /
c     data imach( 5) /  48 /
c     data imach( 6) /   6 /
c     data imach( 7) /   2 /
c     data imach( 8) /  39 /
c     data imach( 9) / o0007777777777777 /
c     data imach(10) /   8 /
c     data imach(11) /  13 /
c     data imach(12) / -50 /
c     data imach(13) /  76 /
c     data imach(14) /  26 /
c     data imach(15) / -50 /
c     data imach(16) /  76 /
c
c     machine constants for the burroughs 6700/7700 systems.
c
c     data imach( 1) /   5 /
c     data imach( 2) /   6 /
c     data imach( 3) /   7 /
c     data imach( 4) /   6 /
c     data imach( 5) /  48 /
c     data imach( 6) /   6 /
c     data imach( 7) /   2 /
c     data imach( 8) /  39 /
c     data imach( 9) / o0007777777777777 /
c     data imach(10) /   8 /
c     data imach(11) /  13 /
c     data imach(12) / -50 /
c     data imach(13) /  76 /
c     data imach(14) /  26 /
c     data imach(15) / -32754 /
c     data imach(16) /  32780 /
c
c     machine constants for the cdc 6000/7000 series.
c
c     data imach( 1) /    5 /
c     data imach( 2) /    6 /
c     data imach( 3) /    7 /
c     data imach( 4) /6loutput/
c     data imach( 5) /   60 /
c     data imach( 6) /   10 /
c     data imach( 7) /    2 /
c     data imach( 8) /   48 /
c     data imach( 9) / 00007777777777777777b /
c     data imach(10) /    2 /
c     data imach(11) /   48 /
c     data imach(12) / -974 /
c     data imach(13) / 1070 /
c     data imach(14) /   96 /
c     data imach(15) / -927 /
c     data imach(16) / 1070 /
c
c     machine constants for the cray 1
c
c     data imach( 1) /   100 /
c     data imach( 2) /   101 /
c     data imach( 3) /   102 /
c     data imach( 4) /   101 /
c     data imach( 5) /    64 /
c     data imach( 6) /     8 /
c     data imach( 7) /     2 /
c     data imach( 8) /    63 /
c     data imach( 9) /  777777777777777777777b /
c     data imach(10) /     2 /
c     data imach(11) /    48 /
c     data imach(12) / -8192 /
c     data imach(13) /  8191 /
c     data imach(14) /    96 /
c     data imach(15) / -8192 /
c     data imach(16) /  8191 /
c
c     machine constants for the data general eclipse s/200
c
c     data imach( 1) /   11 /
c     data imach( 2) /   12 /
c     data imach( 3) /    8 /
c     data imach( 4) /   10 /
c     data imach( 5) /   16 /
c     data imach( 6) /    2 /
c     data imach( 7) /    2 /
c     data imach( 8) /   15 /
c     data imach( 9) /32767 /
c     data imach(10) /   16 /
c     data imach(11) /    6 /
c     data imach(12) /  -64 /
c     data imach(13) /   63 /
c     data imach(14) /   14 /
c     data imach(15) /  -64 /
c     data imach(16) /   63 /
c
c     machine constants for the harris 220
c
c     data imach( 1) /       5 /
c     data imach( 2) /       6 /
c     data imach( 3) /       0 /
c     data imach( 4) /       6 /
c     data imach( 5) /      24 /
c     data imach( 6) /       3 /
c     data imach( 7) /       2 /
c     data imach( 8) /      23 /
c     data imach( 9) / 8388607 /
c     data imach(10) /       2 /
c     data imach(11) /      23 /
c     data imach(12) /    -127 /
c     data imach(13) /     127 /
c     data imach(14) /      38 /
c     data imach(15) /    -127 /
c     data imach(16) /     127 /
c
c     machine constants for the honeywell 600/6000 series.
c
c     data imach( 1) /    5 /
c     data imach( 2) /    6 /
c     data imach( 3) /   43 /
c     data imach( 4) /    6 /
c     data imach( 5) /   36 /
c     data imach( 6) /    6 /
c     data imach( 7) /    2 /
c     data imach( 8) /   35 /
c     data imach( 9) / o377777777777 /
c     data imach(10) /    2 /
c     data imach(11) /   27 /
c     data imach(12) / -127 /
c     data imach(13) /  127 /
c     data imach(14) /   63 /
c     data imach(15) / -127 /
c     data imach(16) /  127 /
c
c     machine constants for the hp 2100
c     3 word double precision option with ftn4
c
c     data imach(1) /      5/
c     data imach(2) /      6 /
c     data imach(3) /      4 /
c     data imach(4) /      1 /
c     data imach(5) /     16 /
c     data imach(6) /      2 /
c     data imach(7) /      2 /
c     data imach(8) /     15 /
c     data imach(9) /  32767 /
c     data imach(10)/      2 /
c     data imach(11)/     23 /
c     data imach(12)/   -128 /
c     data imach(13)/    127 /
c     data imach(14)/     39 /
c     data imach(15)/   -128 /
c     data imach(16)/    127 /
c
c     machine constants for the hp 2100
c     4 word double precision option with ftn4
c
c     data imach(1) /      5 /
c     data imach(2) /      6 /
c     data imach(3) /      4 /
c     data imach(4) /      1 /
c     data imach(5) /     16 /
c     data imach(6) /      2 /
c     data imach(7) /      2 /
c     data imach(8) /     15 /
c     data imach(9) /  32767 /
c     data imach(10)/      2 /
c     data imach(11)/     23 /
c     data imach(12)/   -128 /
c     data imach(13)/    127 /
c     data imach(14)/     55 /
c     data imach(15)/   -128 /
c     data imach(16)/    127 /
c
c     machine constants for the ibm 360/370 series,
c     the xerox sigma 5/7/9, the sel systems 85/86, and
c     the perkin elmer (interdata) 7/32.
c
c     data imach( 1) /   5 /
c     data imach( 2) /   6 /
c     data imach( 3) /   7 /
c     data imach( 4) /   6 /
c     data imach( 5) /  32 /
c     data imach( 6) /   4 /
c     data imach( 7) /   2 /
c     data imach( 8) /  31 /
c     data imach( 9) / z7fffffff /
c     data imach(10) /  16 /
c     data imach(11) /   6 /
c     data imach(12) / -64 /
c     data imach(13) /  63 /
c     data imach(14) /  14 /
c     data imach(15) / -64 /
c     data imach(16) /  63 /
c
c     machine constants for the pdp-10 (ka processor).
c
c     data imach( 1) /    5 /
c     data imach( 2) /    6 /
c     data imach( 3) /    5 /
c     data imach( 4) /    6 /
c     data imach( 5) /   36 /
c     data imach( 6) /    5 /
c     data imach( 7) /    2 /
c     data imach( 8) /   35 /
c     data imach( 9) / "377777777777 /
c     data imach(10) /    2 /
c     data imach(11) /   27 /
c     data imach(12) / -128 /
c     data imach(13) /  127 /
c     data imach(14) /   54 /
c     data imach(15) / -101 /
c     data imach(16) /  127 /
c
c     machine constants for the pdp-10 (ki processor).
c
c     data imach( 1) /    5 /
c     data imach( 2) /    6 /
c     data imach( 3) /    5 /
c     data imach( 4) /    6 /
c     data imach( 5) /   36 /
c     data imach( 6) /    5 /
c     data imach( 7) /    2 /
c     data imach( 8) /   35 /
c     data imach( 9) / "377777777777 /
c     data imach(10) /    2 /
c     data imach(11) /   27 /
c     data imach(12) / -128 /
c     data imach(13) /  127 /
c     data imach(14) /   62 /
c     data imach(15) / -128 /
c     data imach(16) /  127 /
c
c     machine constants for pdp-11 fortran[s supporting
c     32-bit integer arithmetic.
c
c     data imach( 1) /    5 /
c     data imach( 2) /    6 /
c     data imach( 3) /    5 /
c     data imach( 4) /    6 /
c     data imach( 5) /   32 /
c     data imach( 6) /    4 /
c     data imach( 7) /    2 /
c     data imach( 8) /   31 /
c     data imach( 9) / 2147483647 /
c     data imach(10) /    2 /
c     data imach(11) /   24 /
c     data imach(12) / -127 /
c     data imach(13) /  127 /
c     data imach(14) /   56 /
c     data imach(15) / -127 /
c     data imach(16) /  127 /
c
c     machine constants for pdp-11 fortran[s supporting
c     16-bit integer arithmetic.
c
c     data imach( 1) /    5 /
c     data imach( 2) /    6 /
c     data imach( 3) /    5 /
c     data imach( 4) /    6 /
c     data imach( 5) /   16 /
c     data imach( 6) /    2 /
c     data imach( 7) /    2 /
c     data imach( 8) /   15 /
c     data imach( 9) / 32767 /
c     data imach(10) /    2 /
c     data imach(11) /   24 /
c     data imach(12) / -127 /
c     data imach(13) /  127 /
c     data imach(14) /   56 /
c     data imach(15) / -127 /
c     data imach(16) /  127 /
c
c     machine constants for the univac 1100 series. ftn compiler
c
c
c     data imach( 1) /    5 /
c     data imach( 2) /    6 /
c     data imach( 3) /    1 /
c     data imach( 4) /    6 /
c     data imach( 5) /   36 /
c     data imach( 6) /    4 /
c     data imach( 7) /    2 /
c     data imach( 8) /   35 /
c     data imach( 9) / o377777777777 /
c     data imach(10) /    2 /
c     data imach(11) /   27 /
c     data imach(12) / -128 /
c     data imach(13) /  127 /
c     data imach(14) /   60 /
c     data imach(15) /-1024 /
c     data imach(16) / 1023 /
c
c     machine constants for the univac 1100 series. for compiler
c
c     data imach( 1) /    5 /
c     data imach( 2) /    6 /
c     data imach( 3) /    7 /
c     data imach( 4) /    6 /
c     data imach( 5) /   36 /
c     data imach( 6) /    6 /
c     data imach( 7) /    2 /
c     data imach( 8) /   35 /
c     data imach( 9) / o377777777777 /
c     data imach(10) /    2 /
c     data imach(11) /   27 /
c     data imach(12) / -128 /
c     data imach(13) /  127 /
c     data imach(14) /   60 /
c     data imach(15) /-1024/
c     data imach(16) / 1023 /
c
c
c     machine constants for the vax 11/780
c
c     data imach(1) /    5 /
c     data imach(2) /    6 /
c     data imach(3) /    5 /
c     data imach(4) /    6 /
c     data imach(5) /   32 /
c     data imach(6) /    4 /
c     data imach(7) /    2 /
c     data imach(8) /   31 /
c     data imach(9) /2147483647 /
c     data imach(10)/    2 /
c     data imach(11)/   24 /
c     data imach(12)/ -127 /
c     data imach(13)/  127 /
c     data imach(14)/   56 /
c     data imach(15)/ -127 /
c     data imach(16)/  127 /
c
c     machine constants for the ridge 32

      data imach(1)  /      5 /
      data imach(2)  /      6 /
      data imach(3)  /      6 /
      data imach(4)  /      6 /
      data imach(5)  /     32 /
      data imach(6)  /      4 /
      data imach(7)  /      2 /
      data imach(8)  /     31 /
      data imach(9)  / 2147483647 /
      data imach(10) /      2 /
      data imach(11) /     23 /
      data imach(12) /   -127 /
      data imach(13) /    128 /
      data imach(14) /     53 /
      data imach(15) /  -1023 /
      data imach(16) /   1024 /
c
c***first executable statement  i1mach
c
      if (i .lt. 1  .or.  i .gt. 16) go to 10
c
      i1mach=imach(i)
      return
c
   10 continue
      write(output,9000)
9000  format(39h1error    1 in i1mach - i out of bounds )
c
c     call fdump
c
c
      stop
      end
      integer function icamax(n,cx,incx)
c***begin prologue  icamax
c***revision date  811015   (yymmdd)
c***category no.  f1a
c***keywords  blas,vector,complex,component,largest component
c***date written  october 1979
c
c***changed to double precision june 1985 by m.s.
c
c***author lawson c. (jpl),hanson r. (sla),
c                            kincaid d. (u texas), krogh f. (jpl)
c***purpose
c    find largest component of complex vector
c***description
c                b l a s  subprogram
c    description of parameters
c
c     --input--
c        n  number of elements in input vector(s)
c       cx  complex vector with n elements
c     incx  storage spacing between elements of cx
c
c     --output--
c   icamax  smallest index (zero if n.le.0)
c
c      returns the index of the component of cx having the
c      largest sum of magnitudes of real and imaginary parts.
c     icamax = first i, i = 1 to n, to minimize
c        abs(real(cx(1-incx+i*incx))) + abs(imag(cx(1-incx+i*incx)))
c
c
c***references
c  lawson c.l., hanson r.j., kincaid d.r., krogh f.t.,
c   *basic linear algebra subprograms for fortran usage*,
c  algorithm no. 539, transactions on mathematical software,
c  volume 5, number 3, september 1979, 308-323
c***routines called  (none)
c***end prologue  icamax
c
      implicit real*8 (a-h,o-z)
      complex*16 cx(1)
c***first executable statement  icamax
      icamax = 0
      if(n.le.0) return
      icamax = 1
      if(n .le. 1) return
      ns = n*incx
      ii = 1
      summax = abs(dble(cx(1))) + abs(dimag(cx(1)))
	  do 20 i=1,ns,incx
	  sumri = abs(dble(cx(i))) + abs(dimag(cx(i)))
	  if(summax.ge.sumri) go to 10
	  summax = sumri
	  icamax = ii
   10     ii = ii + 1
   20     continue
      return
      end
      function j4save(iwhich,ivalue,iset)
c***begin prologue  j4save
c
c***changed to double precision
c
c***refer to  xerror
c     abstract
c        j4save saves and recalls several global variables needed
c        by the library error handling routines.
c
c     description of parameters
c      --input--
c        iwhich - index of item desired.
c                 = 1 refers to current error number.
c                 = 2 refers to current error control flag.
c                 = 3 refers to current unit number to which error
c                     messages are to be sent.  (0 means use standard.)
c                 = 4 refers to the maximum number of times any
c                     message is to be printed (as set by xermax).
c                 = 5 refers to the total number of units to which
c                     each error message is to be written.
c                 = 6 refers to the 2nd unit for error messages
c                 = 7 refers to the 3rd unit for error messages
c                 = 8 refers to the 4th unit for error messages
c                 = 9 refers to the 5th unit for error messages
c        ivalue - the value to be set for the iwhich-th parameter,
c                 if iset is .true. .
c        iset   - if iset=.true., the iwhich-th parameter will be
c                 given the value, ivalue.  if iset=.false., the
c                 iwhich-th parameter will be unchanged, and ivalue
c                 is a dummy parameter.
c      --output--
c        the (old) value of the iwhich-th parameter will be returned
c        in the function value, j4save.
c
c     written by ron jones, with slatec common math library subcommittee
c     adapted from bell laboratories port library error handler
c     latest revision ---  23 may 1979
c
c***references
c   jones r.e., *slatec common mathematical library error handling
c    package*, sand78-1189, sandia laboratories, 1978.
c***routines called  (none)
c***end prologue  j4save
c
      logical iset
      integer iparam(9)
      data iparam(1),iparam(2),iparam(3),iparam(4)/0,1,0,10/
      data iparam(5)/1/
      data iparam(6),iparam(7),iparam(8),iparam(9)/0,0,0,0/
c***first executable statement  j4save
      j4save = iparam(iwhich)
      if (iset) iparam(iwhich) = ivalue
      return
      end
      subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sigma,wa1,
     *                 wa2)
c***begin prologue  lmpar
c***refer to snlse,snls
c     **********
c
c     subroutine lmpar
c
c     given an m by n matrix a, an n by n nonsingular diagonal
c     matrix d, an m-vector b, and a positive number delta,
c     the problem is to determine a value for the parameter
c     par such that if x solves the system
c
c           a*x = b ,     sqrt(par)*d*x = 0 ,
c
c     in the least squares sense, and dxnorm is the euclidean
c     norm of d*x, then either par is zero and
c
c           (dxnorm-delta) .le. 0.1*delta ,
c
c     or par is positive and
c
c           abs(dxnorm-delta) .le. 0.1*delta .
c
c     this subroutine completes the solution of the problem
c     if it is provided with the necessary information from the
c     qr factorization, with column pivoting, of a. that is, if
c     a*p = q*r, where p is a permutation matrix, q has orthogonal
c     columns, and r is an upper triangular matrix with diagonal
c     elements of nonincreasing magnitude, then lmpar expects
c     the full upper triangle of r, the permutation matrix p,
c     and the first n components of (q transpose)*b. on output
c     lmpar also provides an upper triangular matrix s such that
c
c            t   t                   t
c           p *(a *a + par*d*d)*p = s *s .
c
c     s is employed within lmpar and may be of separate interest.
c
c     only a few iterations are generally needed for convergence
c     of the algorithm. if, however, the limit of 10 iterations
c     is reached, then the output par will contain the best
c     value obtained so far.
c
c     the subroutine statement is
c
c       subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sigma,
c                        wa1,wa2)
c
c     where
c
c       n is a positive integer input variable set to the order of r.
c
c       r is an n by n array. on input the full upper triangle
c         must contain the full upper triangle of the matrix r.
c         on output the full upper triangle is unaltered, and the
c         strict lower triangle contains the strict upper triangle
c         (transposed) of the upper triangular matrix s.
c
c       ldr is a positive integer input variable not less than n
c         which specifies the leading dimension of the array r.
c
c       ipvt is an integer input array of length n which defines the
c         permutation matrix p such that a*p = q*r. column j of p
c         is column ipvt(j) of the identity matrix.
c
c       diag is an input array of length n which must contain the
c         diagonal elements of the matrix d.
c
c       qtb is an input array of length n which must contain the first
c         n elements of the vector (q transpose)*b.
c
c       delta is a positive input variable which specifies an upper
c         bound on the euclidean norm of d*x.
c
c       par is a nonnegative variable. on input par contains an
c         initial estimate of the levenberg-marquardt parameter.
c         on output par contains the final estimate.
c
c       x is an output array of length n which contains the least
c         squares solution of the system a*x = b, sqrt(par)*d*x = 0,
c         for the output par.
c
c       sigma is an output array of length n which contains the
c         diagonal elements of the upper triangular matrix s.
c
c       wa1 and wa2 are work arrays of length n.
c
c     subprograms called
c
c       minpack-supplied ... r1mach,enorm,qrsolv
c
c       fortran-supplied ... abs,dmax1,dmin1,sqrt
c
c     minpack. version of june 1979.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
c
c     dwarf is the smallest positive magnitude.
c
c***routines called  r1mach,enorm,qrsolv
c***end prologue  lmpar
      integer n,ldr
      integer ipvt(n)
      double precision delta,par
      double precision r(ldr,n),diag(n),qtb(n),x(n),
     &                 sigma(n),wa1(n),wa2(n)
      integer i,iter,j,jm1,jp1,k,l,nsing
      double precision dxnorm,dwarf,fp,gnorm,parc,parl,
     &                 paru,p1,p001,sum,temp,zero
      double precision r1mach,enorm
      data p1,p001,zero /1.0d-1,1.0d-3,0.0d0/
c***first executable statement  lmpar
      dwarf = r1mach(1)
c
c     compute and store in x the gauss-newton direction. if the
c     jacobian is rank-deficient, obtain a least squares solution.
c
      nsing = n
      do 10 j = 1, n
	 wa1(j) = qtb(j)
	 if (r(j,j) .eq. zero .and. nsing .eq. n) nsing = j - 1
	 if (nsing .lt. n) wa1(j) = zero
   10    continue
      if (nsing .lt. 1) go to 50
      do 40 k = 1, nsing
	 j = nsing - k + 1
	 wa1(j) = wa1(j)/r(j,j)
	 temp = wa1(j)
	 jm1 = j - 1
	 if (jm1 .lt. 1) go to 30
	 do 20 i = 1, jm1
	    wa1(i) = wa1(i) - r(i,j)*temp
   20       continue
   30    continue
   40    continue
   50 continue
      do 60 j = 1, n
	 l = ipvt(j)
	 x(l) = wa1(j)
   60    continue
c
c     initialize the iteration counter.
c     evaluate the function at the origin, and test
c     for acceptance of the gauss-newton direction.
c
      iter = 0
      do 70 j = 1, n
	 wa2(j) = diag(j)*x(j)
   70    continue
      dxnorm = enorm(n,wa2)
      fp = dxnorm - delta
      if (fp .le. p1*delta) go to 220
c
c     if the jacobian is not rank deficient, the newton
c     step provides a lower bound, parl, for the zero of
c     the function. otherwise set this bound to zero.
c
      parl = zero
      if (nsing .lt. n) go to 120
      do 80 j = 1, n
	 l = ipvt(j)
	 wa1(j) = diag(l)*(wa2(l)/dxnorm)
   80    continue
      do 110 j = 1, n
	 sum = zero
	 jm1 = j - 1
	 if (jm1 .lt. 1) go to 100
	 do 90 i = 1, jm1
	    sum = sum + r(i,j)*wa1(i)
   90       continue
  100    continue
	 wa1(j) = (wa1(j) - sum)/r(j,j)
  110    continue
      temp = enorm(n,wa1)
      parl = ((fp/delta)/temp)/temp
  120 continue
c
c     calculate an upper bound, paru, for the zero of the function.
c
      do 140 j = 1, n
	 sum = zero
	 do 130 i = 1, j
	    sum = sum + r(i,j)*qtb(i)
  130       continue
	 l = ipvt(j)
	 wa1(j) = sum/diag(l)
  140    continue
      gnorm = enorm(n,wa1)
      paru = gnorm/delta
      if (paru .eq. zero) paru = dwarf/dmin1(delta,p1)
c
c     if the input par lies outside of the interval (parl,paru),
c     set par to the closer endpoint.
c
      par = dmax1(par,parl)
      par = dmin1(par,paru)
      if (par .eq. zero) par = gnorm/dxnorm
c
c     beginning of an iteration.
c
  150 continue
	 iter = iter + 1
c
c        evaluate the function at the current value of par.
c
	 if (par .eq. zero) par = dmax1(dwarf,p001*paru)
	 temp = sqrt(par)
	 do 160 j = 1, n
	    wa1(j) = temp*diag(j)
  160       continue
	 call qrsolv(n,r,ldr,ipvt,wa1,qtb,x,sigma,wa2)
	 do 170 j = 1, n
	    wa2(j) = diag(j)*x(j)
  170       continue
	 dxnorm = enorm(n,wa2)
	 temp = fp
	 fp = dxnorm - delta
c
c        if the function is small enough, accept the current value
c        of par. also test for the exceptional cases where parl
c        is zero or the number of iterations has reached 10.
c
	 if (abs(fp) .le. p1*delta
     *       .or. parl .eq. zero .and. fp .le. temp
     *            .and. temp .lt. zero .or. iter .eq. 10) go to 220
c
c        compute the newton correction.
c
	 do 180 j = 1, n
	    l = ipvt(j)
	    wa1(j) = diag(l)*(wa2(l)/dxnorm)
  180       continue
	 do 210 j = 1, n
	    wa1(j) = wa1(j)/sigma(j)
	    temp = wa1(j)
	    jp1 = j + 1
	    if (n .lt. jp1) go to 200
	    do 190 i = jp1, n
	       wa1(i) = wa1(i) - r(i,j)*temp
  190          continue
  200       continue
  210       continue
	 temp = enorm(n,wa1)
	 parc = ((fp/delta)/temp)/temp
c
c        depending on the sign of the function, update parl or paru.
c
	 if (fp .gt. zero) parl = dmax1(parl,par)
	 if (fp .lt. zero) paru = dmin1(paru,par)
c
c        compute an improved estimate for par.
c
	 par = dmax1(parl,par+parc)
c
c        end of an iteration.
c
	 go to 150
  220 continue
c
c     termination.
c
      if (iter .eq. 0) par = zero
      return
c
c     last card of subroutine lmpar.
c
      end
c/* @(#)optpbs.f        1.7      latest revision 10/31/86 11:52:23 */
c/* %w%      latest revision %g% %u% */
      subroutine optpbs (kk,n1,wt,lpia,lmaxin,lborn,fx,er,eim)
c      coupled channels version...see dimens + equiv
c       modified sept85 for anti protons, n.b. ze now passed
c *** see r h landau, program lpott, 1981
      implicit real*8 (a-h,k,o-z)
      real*8 imtij,mpi,mn,k,kp,kk,imvabs,k1
      complex*16 vkp(3,3),zi,zk,zkp,zk2,zu(192,192),ze
      dimension rho(4,50), bi(100), vcoull(30), revabs(30,2), imvabs(30,
     12) ,wt(64)
      dimension nifty(20), bnuc(20), bnucf(20), bn(16), bnf(16), retij
     1(14),imtij(14),kk(64),u(36864,2),f(18432,2),fx(30),
     2 ui(192,192),ur(192,192),vr(6),vi(6)
       dimension ff(715,6)
      dimension np(5)
      common /sec2/ bnuc,bn,bnucf,bnf,retij,imtij,hbarc,pi,mpi,mn,nz,nes
     1,nwaves,nifty,na
      common/nlspfl/u/foptp/f
      common /sec4/ achp,acmp,wsp,achn,acmn,wsn
      common /sec5/ drho(4,50),rho,xkp,xk,imax/sec6/rcoul,rcut
      equivalence (vkp(1,1),bnuc(1)),(iset,nes) ,(ur(1,1),u(1,1)),
     1(ui(1,1),u(1,2)),(f(1,1),vr(1)),(f(7,1),vi(1)),(f(1,1),ff(1,1))
       equivalence(zu(1,1),u(1,1))
      data vcoull,revabs,imvabs/150*0.0e+00/,ikp,zi/0,(0.e+00,1.e+00)/
     1 ,nint,eri,eimi/0,0.,0./
c
c        need reset dimensns and equivs as max no grid increase
c      should use n=ngp +1 for scattering
c-----1st  equiv----------------------------------------------------------
c      ur+ui(3n,3n) , u( (3n)**2,2)
c-----2nd  equivalence----------------------------------------------------
c      f(n(n+1)*2)-in main      f(n(n+1),2)-in optpbs
c      ff(ndimf/6, 6) also needed ,and equivalenced to f(1)
c     for equiv need length+1 to overlay
c     also see intialize statements which now change with dimensions
c     calculate partial wave decomp of rho- the form factor
c     rho(nff,i) is the i-1 th patial wave ampltd of the form factor
c     for a value of k a nd kp. nff determines which f.f. is used
c      nff=1 proton matter     2 neutron matter
c      3 protone spin  4 neutron spin
c     do all lvalues on 1stcall,then store on file7, ft is for temp
c
c       comple k2 needed for kp channel 2-open,fx(1),(2)=re,im k
c     fx(3) lt 1 define ur,ui       ge 1    define complex*16 zu
c     zk2=fx(1)+zi*fx(2)
      nx=fx(3)
      nint=nint+1
	if(er.eq.eri.and.eim.eq.eimi)goto 9
c       energy has changed, recalculate u for e-dependent potential
c
	eri=er
	eimi=eim
	go to 8
    9 if (lpia.ne.1.or.nint.ne.1) go to 250
c     1st call to optpbs, calculate all u s
c      intialize (cant use data for this due to commoms)
    8 do 7 i=1,2
      do 7 j=1,18432
    7 f(j,i)=0.
      do 6 j=1,36864
      u(j,1)=0.
    6 u(j,2)=0.
      if(kk(n1).gt.0.)k1=kk(n1)
      if(nifty(10).gt.4)k1=kk(1)
      if (nint.eq.1) write(6,410) k1
      k = -200.
c     ikp = k- p switch for special sub for opt pot
c       also for pbar-p
      if(nifty(1).ge.10.and.na.eq.1)ikp=1
      if(ikp.eq.1) go to 5
c     if(nifty(6).ne.4.and.nifty(6).ne.6)call tncm (k)
  5   n2 = 2*n1
      nn = na-nz
c      make sure lmax=1 case does not run into indef on intial runs
      ldmax = lmaxin
      if (ldmax.eq.1) ldmax = 2
      rhl3 = (achp/hbarc)**2/2.
      rhl4 = (achn/hbarc)**2/2.
c     addeded terms for c13 spin...may want to change size  here
      rhl5 = (wsn/hbarc)**2/2.
      as = 1./6.
      if(nz.ne.0)ap = (nz-2.)/6./nz*(acmp/achp)**2
      if (nz.le.2) ap = 0.
      if(nn.ne.0)an = (nn-2.)/6./nn*(acmn/achn)**2
      if (nn.le.2) an = 0.
c     special c13 , add in p1/2 n and rescale
      if (na.eq.13) an = (nn-3.)/6./(nn-1.)*(acmn/achn)**2
c     do loop over all grid points
c     size of u must be 4*(ngp+1)**2,2
      rewind 7
      if (nifty(6).eq.6) nifty(17) = 1
      aovera = (na-1.)/na
      if (nifty(17).eq.1) aovera = 1.
      if (na.le.2) aovera = 1.
c      set u(strong)=0
      if ( nifty(6).eq.4)  aovera = 0.
      xhl = aovera*2.*pi*pi
      aparam = 1.
      i = nifty(1)+1
      go to (30,10,20,20,20,30,20,30,20,20,10,20,10), i
 10   aparam = -1.
      go to 30
 20   aparam = 0.
 30   if(rcoul.lt.1.e-10 )rcoul=1.e-6
      aparam = (1.5*nz*aparam/rcoul)*hbarc*7.29735e-3
      xhl1 = rcoul/hbarc
      xhl2 = rcut/hbarc
      nif = nifty(10)
      do 240 i1=1,n1
	 do 240 i2=1,i1
	    k = kk(i2)
	    kp = kk(i1)
	    lmax = ldmax+5
	    if(nifty(6).eq.6) go to 110
	    if (nifty(16).ne.0) go to 90
      if(ikp.gt.0)go to 110
c      fo shell model nuclesu(c,0) use analytic expression with il*s
c     different neutron and protone densiteies possible
	    zp = k*kp*rhl3
	    xp = (k*k+kp*kp)*rhl3/2.
c           call beslia (zp,bi)
	    b = 0.
	    if ((zp-xp).gt.-150.) b = exp(zp-xp)
	    do 40 i=1,lmax
	       rho(3,i) = 0.
	       rho(1,i) = b*(2.*i-1)*((1.-4.*ap*(xp-i+1))*bi(i)+4.*zp*ap
     1         *bi(i+1))
 40         continue
	    if (na.eq.13) go to 60
	    if ((achp.ne.achn).or.(acmp.ne.acmn)) go to 60
c     equal n and p densities
	    do 50 i=1,lmax
	       rho(4,i) = 0.
	       rho(2,i) = rho(1,i)
 50         continue
	    nff = 1
	    go to 110
c     unequal n and p densites in h.o. model
 60         nff = 2
	    zn = k*kp*rhl4
	    xn = (k*k+kp*kp)*rhl4/2.
c           call beslia (zn,bi)
	    b = 0.
	    if ((zn-xn).gt.-150.) b = exp(zn-xn)
	    do 70 i=1,lmax
	       rho(2,i) = b*(2.*i-1)*((1.-4.*an*(xn-i+1))*bi(i)+4.*zn*an
     1         *bi(i+1))
	       rho(4,i) = 0.
 70         continue
	    if (na.ne.13) go to 110
c     spin for c13
	    nff = 4
	    zs = k*kp*rhl5
	    xs = (k*k+kp*kp)*rhl5/2.
c           call beslia (zs,bi)
	    b = 0.
	    if ((zs-xs).gt.-150.) b = exp(zs-xs)
c     normalize to 1/nn so can mult by nn and have 1 neutron
	    b = b/nn
	    do 80 i=1,lmax
	       rho(4,i) = b*(2.*i-1)*((1.-4.*as*(xs-i+1))*bi(i)+4.*zs*as
     1         *bi(i+1))
	       rho(2,i) = ((nn-1.)/nn)*rho(2,i)+rho(4,i)
 80         continue
	    if ((i1.eq.i2).and.(i1.eq.n1).and.nint.eq.1) 
     1       write (6,350) k,(rho(4,i),i=1,lmax)
	    go to 110
c      for wood saxon typr(or arbitrary) nuclera n density,calc rho nume
 90         xkp = kp
	    xk = k
	    imax = ldmax+5
	    lmax = imax
c     deteermine number of f.f.*s needed,if nifty(6)=3, need spin in
	    nff = 2
	    if ((achp.eq.achn).and.(acmp.eq.acmn).and.(wsp.eq.wsn)) nff
     1      = 1
	    if (nifty(6).eq.3) nff = 4
c      for rho**2 in ca
	    if (nifty(16).eq.2) nff = 4
c           call rhokff (nff)
	    if (nff.ne.1) go to 110
c     nff=1, seet n and p matter densites=, all spin dens=0
	    do 100 i=1,lmax
	       rho(2,i) = rho(1,i)
	       rho(3,i) = 0.
	       rho(4,i) = 0.
 100        continue
 110        continue
	    if (nifty(10).lt.3.or.nifty(10).eq.8) go to 120
c     include cutoff coulomb potential in momentum space
c      nb vl contains scale factor ,i.e. it is v(j=l+1/2)
	    ak = k
	    akp = kp
c             dec83 add,no coul for open chan
	  if(i1.eq.n1.or.i2.eq.n1)goto120
	  if(ak.le.0..or.akp.le.0.)go to 120
	  if(i1.ne.i2) go to 118
c         k = kp, need special vcoul for bs so remove singularity
	  if(nifty(10).lt.5) go to 118
	  call vcdiag(ak,i2,aparam,xhl1,xhl2,vcoull,ldmax,
     1 nif,n1,kk,wt)
	  go to 120
 118        call vcoulb (ak,akp,aparam,xhl1,xhl2,vcoull,ldmax,nif)
 120  if (nifty(7) .lt. 3 .and. nifty(7).ne.5) go to 130
c      include absorption on 2 nucleons
c           call vabs (k,kp,kk(n1),revabs,imvabs,ldmax)
 130        continue
c     new version....outter do loopn over kp,k,,innerover all l s
	    xmn = mn
c        if(nifty(6).ne.6.and.ikp.ne.1)call tnuc (kp,k,kk(n1),xmn)
c     do loop over l                                                 opt
c
	    do 230 ldum=1,ldmax
	       if ((ldum.gt.lborn).and.((i1+i2).ne.(2*n1))) go to 230
	       j1 = ldum
	       urc = 0.
	     vc=0
	       uic = 0.
	       urs = 0.
	       uis = 0.
	       if (nifty(6).eq.3) rhl1 = j1*(j1-1.)
	    if(nifty(6).eq.6.or.ikp.eq.1) go to 205
c     do over pi-n angular momentum;s,p,d,f,g,
	       do 200 j2=1,5
		  j4 = iabs(j1-j2)+1
		  j5 = j1+j2-1
		  if (nifty(6).eq.3) rhl2 = j2*(j2-1.)
		  do 190 j3=j4,j5
c******************* rhl = cgc2(j3,j2,j1)
c     spin independt potential
		     if (nifty(1).eq.1.or.nifty(1).eq.8) go to 140
c     pi plus scatterinf
		     jpr = 4*j2-3
		     jpi = jpr+1
		     jnr = jpr+2
		     jni = jpr+3
		     go to 150
 140                 continue
c     pi minus scatterin, rearrange code
		     jnr = 4*j2-3
		     jni = jnr+1
		     jpr = jnr+2
		     jpi = jnr+3
 150                 continue
		     if (nifty(1).ne.4) go to 160
c     pi zero scattering,  take average of p and n amplitudes
		     bnuc(jpr) = (bnuc(jpr)+bnuc(jnr))/2.
		     bnuc(jnr) = bnuc(jpr)
		     bnuc(jpi) = (bnuc(jpi)+bnuc(jni))/2.
		     bnuc(jni) = bnuc(jpi)
		     if (nifty(6).ne.3) go to 160
		     bnucf(jpr) = (bnucf(jpr)+bnucf(jnr))/2.
		     bnucf(jnr) = bnucf(jpr)
		     bnucf(jpi) = (bnucf(jpi)+bnucf(jni))/2.
		     bnucf(jni) = bnucf(jpi)
 160                 continue
		     urc = urc+rhl*(nz*rho(1,j3)*bnuc(jpr)+nn*rho(2,j3)*
     1               bnuc(jnr))
		     uic = uic+rhl*(nz*rho(1,j3)*bnuc(jpi)+nn*rho(2,j3)*
     1               bnuc(jni))
 180                 if (nifty(6).ne.3) go to 190
c     ---- ---spin dependnet potential
		     if (j1.eq.1) go to 190
		     rhl = rhl*(rhl2+rhl1-j3*(j3-1.))/(2.*rhl1)
		     urs = urs+rhl*(nz*rho(3,j3)*bnucf(jpr)+nn*rho(4,j3)
     1               *bnucf(jnr))
		     uis = uis+rhl*(nz*rho(3,j3)*bnucf(jpi)+nn*rho(4,j3)
     1               *bnucf(jni))
 190              continue
c          rhl=(2pi**2/2l+1)(a-1/a)
 200           continue
 205             continue
c     form potentials for state of definite j
	       rhl = xhl/(2.*j1-1.)
	       vc = 0.
	       if (nifty(10).lt.3.or.nifty(10).eq.8) go to 210
	       vc = vcoull(j1)*aovera
	       if (nifty(6).eq.6) vc = vcoull(j1)
 210           if (nifty(7).ne.3.and.nifty(7).ne.4) go to 220
c     add in vabs to central potential
c     aovera scaling already done for vabs in sub vabs
	       urc = urc+revabs(j1,1)/rhl
	       uic = uic+imvabs(j1,1)/rhl
	       if (nifty(7).ne.4.or.nifty(6).ne.3) go to 220
c     include spin-dependent absorption
	       urs = urs+revabs(j1,2)/rhl
	       uis = uis+imvabs(j1,2)/rhl
 220           continue
c     set up potential needed for this l in integral eqtn
c     nb the f is - imv as stored, latter reversed
	       npot1 = 2*ldum
	       f(npot1-1,1) = rhl*(urc+(j1-1)*urs)+vc
	       f(npot1,1) = -rhl*(uic+(j1-1)*uis)
       ff1=f(npot1-1,1)
       ff2=f(npot1,1)
c      pure coulomb or non k-p
	     if(nifty(6).eq.6 .or. ikp. ne .1) go to 226
c
c      special k-p/or pbar p/section coupled or uncoupled 
c      nes=iset=switch for model (k-p) or channel(pbar) 
c      only l=0 for k-, yet l=1 for pbar. user beware!
 225  zkp=kp
      zk= k
  227   if(nifty(1).eq.12)then
c       pbar-p section, need use complex on-shell energy
c       can change to complex momentum too
	ze=er+zi*eim
	call vpbar(kp,k,ze,ldum)
	else
	call vkbarp(zkp,zk)
	endif
      f(npot1-1,1)=vkp(1,1)+f(npot1-1,1)
      ff1=f(npot1-1,1)
c    nb reverse sign of v here for storage as f
      f(npot1  ,1)=+zi*vkp(1,1)+f(npot1  ,1)
	ff2=f(npot1,1)
      if(nifty(15).eq.0) go to 226
	ff(npot1-1,1)=ff1
	ff(npot1,1)= ff2
c      coupled channels kp
c      store matrices with new, nonlinear forms
c      not -imv any longer  so need add 2 imv to correct the above f
c      nifty(15)=0 ..uncoupled    1...k-p+sigma-pi
c           2...k-p,sig pi,ko/ n  3...k-p,  k0/ n
c      nspin for coupled kp sigma-pi kobarn
c         1---v11  2---v22  3---v33
c         4---v12  5---v13  6---v23
c
      ff2=ff2-2.*zi*vkp(1,1)
	ff(npot1,1)=ff2
      i=nifty(15)
      go to (231,231,232),i
c        k-p plus sigma pi (coul only for k-p)
 231  ff(npot1-1,2)=vkp(2,2)
      ff(npot1,2)=-zi*vkp(2,2)
      ff(npot1-1,4)=vkp(1,2)
      ff(npot1,4)=-zi*vkp(1,2)
      if(nifty(15).eq.1) go to 226
c      kip, sigma-pi, kobar n
 232  ff(npot1-1,3)=vkp(3,3)
      ff(npot1,3)=-zi*vkp(3,3)
      ff(npot1-1,5)=vkp(1,3)
      ff(npot1,5)=-zi*vkp(1,3)
      if (nifty(15).eq.3) go to 226
      ff(npot1-1,6)=vkp(2,3)
      ff(npot1,6)=-zi*vkp(2,3)
 226           nnn = npot1-1
		 if(kp.eq.k1.and.kp.eq.k.and.ldum.eq.1.and.nint.eq.1)
     1           write(6,410)k1
	if(k.eq.k1.and.ikp.eq.1.and.lpia.eq.2.and.nint.eq.1
     1     .and.ldum.eq.2)goto228
	if(k.eq.k1.and.ikp.eq.1.and.lpia.eq.2.and.ldum.eq.1
     1     .and.ldum.eq.1)goto402
       if(k.ne.k1.or.ldum.ne.1.or.nint.ne.1) go to 402
c      write out potential matrices on printer (half offshell)
 228            write (6,400)kp,ff1
     1         ,vc,revabs(j1,1),revabs(j1,2),ff2,imvabs(j1,1),
     2         imvabs(j1,2)
      n=nifty(15)
      if(i1.ne.1.and.i1.ne.8) go to 402
      if(n.eq.1.or.n.eq.2)write(6,401)(i,ff(nnn,i),ff(npot1,i),i=2,4,2)
      if(n.eq.2.or.n.eq.3)write(6,401)(i,ff(nnn,i),ff(npot1,i),i=3,5,2)
      if(n.eq.2          )write(6,401)(i,ff(nnn,i),ff(npot1,i),i=6,6,1)
 401  format(1x,i2,19x,e19.4,41x,e19.4)
 402           if (nifty(6).ne.3) go to 230
	       f(npot1-1,2) = rhl*(urc-j1*urs)+vc
	       f(npot1,2) = -rhl*(uic-j1*uis)
c     do loop over ldum ends
 230        continue
	    npot1m = 2*lborn
	    if (ldmax.lt.lborn) npot1m = 2*ldmax
	    if ((i1.eq.i2).and.(i1.eq.n1)) npot1m = 2*ldmax
c      new unformatted read for time saving
	    if (nifty(6).ne.3.and.nifty(15).eq.0) write (7)
     1 (f(npot1,1),npot1=1,npot1m)
	    if(nifty(6).eq.3)write (7)((f(npot2-1,i),f(npot2,i),i=1,2
     1      ),npot2=2,npot1m,2)
      if(nifty(15).ne.0)write(7) ((ff(npot2-1,i),ff(npot2,i),i=1,6
     1 ),npot2=2,npot1m,2)
c     do loop over kp and k ends
 240  continue
c
c----------------read u from unit7 on nonintial calls
c
 250  rewind 7
c     pots all stored,read in the values for cuurrent lpia
c     npts =? of i1,i2 combos(followed by u s for all l s)
      npts = (n1*(1+n1))/2
c       mod 9/85 so dont need lpia=1 to set up
c     n.b....the values of vl(i,j) for all l but fixed
c     i,j are on one record via unformatted write
c      need skip to new record to get new ij,
      nspin = 1
      npot = 2
      if(nifty(15).ne.0)nspin=6
      if (nifty(6).ne.3) go to 260
c     spin-dependent case, now 4 v s adjacent for each l vs 2 if nf
      nspin = 2
      npot = 2*npot
c      coupled channels kn, 12 v s adjacent,non imv sign reversal
 260  ldum = 2*(lpia-1)*nspin
       if(nifty(15).ne.0) go to 345
      if (lpia.gt.lborn) go to 290
c     do loop over npts (i.e. records)
c       nb ldum contains both no of l's and nspin(pot types) to skip
      do 280 npt=1,npts
	 if (lpia.ne.1) go to 270
	 read (7) (f(npt,i),i=1,npot)
	 go to 280
c        space to l value needed
 270     read(7)(rhl,i=1,ldum),(f(npt,i),i=1,npot)
 280  continue
      go to 310
c     lpia .gt. lborn, read only last on-shell point
 290  nskip = npts-1
      do 300 i=1,nskip
	 read (7)
 300  continue
      read (7) (rhl,i=1,ldum),(f(npts,i),i=1,npot)
c     set up potential matrix
 310  istart = 1
      if (lpia.gt.lborn) istart = n1
      do 340 i1=istart,n1
	 do 340 i2=istart,i1
	    iflag = 0
	    npt = ((i1-1)*i1)/2+i2
      if(nifty(10).ge.5)go to 331
	    npot1 = (i2-1)*n2+i1
 320        npot2 = npot1+n1*(n2+1)
c      imv:npot4    -imv:npot3 (the one stored for uncoupled channels)
	    npot3 = npot2-n1
	    npot4 = npot1+n1
	    u(npot1,1) = f(npt,1)
	    u(npot2,1) = u(npot1,1)
c      rhl change feb82 for bs
      if(nifty(6).eq.7)f(npt,2)=0.
	    u(npot3,1) = f(npt,2)
	    u(npot4,1) = -f(npt,2)
c     spin-dependnt potentials
	    if (nifty(6).ne.3) go to 330
c      special turn off of spin for bs to avoid mstrix def
       write(6,999)
 999   format(' ********** spin section in opt no good at 999 ')
c           u(npot1,2) = f(npt,3)
c           u(npot2,2) = u(npot1,2)
c           u(npot3,2) = f(npt,4)
	     if(nifty(6).eq.7) u(npot3,2)=0.
	    u(npot4,2) = -u(npot3,2)
c     set up potntial with u(kp,k)=u(k,kp)..symmetric u
 330        if ((i1.eq.i2).or.(iflag.eq.1)) go to 340
	    npot1 = (i1-1)*n2+i2
	    iflag = 1
	    go to 320
c      bound state case, set up seperate ur,ui
 331  if(nx.lt.1)ur(i1,i2)=f(npt,1)
      if(nx.lt.1)ui(i1,i2)=-f(npt,2)
      if(nx.gt.1)zu(i1,i2)=f(npt,1)-zi*f(npt,2)
      if(i1.eq.i2) go to 340
c      off diagnol symmetrty
      if(nx.lt.1)ur(i2,i1)=ur(i1,i2)
      if(nx.lt.1)ui(i2,i1)=ui(i1,i2)
      if(nx.gt.1)zu(i2,i1)=zu(i1,i2)
 340  continue
      return
c
c     special coulpled channels calc....now for k-p,kobar-n, sigma-pi
c      1st fill lower half of full matrix,then latter upper
c then set ur,ui to super matrix (ur,ui equiv to the above u)
c                 v1     v4     v5
c                 v4     v2     v6
c                 v5     v6     v3
c      npt is order written out/read in from storage
c       here only for b.s. with on shell point now included as n1
 345  ngp=n1
      if(nifty(15).eq.3)ngp=n1-1
      ng2=2*ngp
      ng3=3*ngp
c      npt is locatn in triangle matrix stored linearly
c             rearrang with (i,j) index and above symmetry
      if(lpia.gt.lborn) write(6,420) lpia,lborn
 420  format(' lpia gt lborn in optpbs bs********',2i3)
      do 347 i1=1,ngp
c     skip on shell point if ko k-
       if(i1.eq.n1.and.nifty(15).eq.3) go to 347
      npot1=i1+ngp
      npot3=i1+ng2
      do 346 i2=1,i1
      npot2=i2+ngp
      npot4=i2+ng2
      if(lpia.eq.1)read(7)(vr(i),vi(i),i=1,6)
      if(lpia.ne.1)read(7)(rhl,i=1,ldum),(vr(i),vi(i),i=1,6)
	if(nifty(6).eq.7)then
c       set imv = 0
	do 425 i=1,6
 425    vi(i)=0.
	endif
c      off diagnol block matrices,build in internal symmetry
      if(nifty(15).eq.3) go to 421
      if(nx.lt.1)ur(npot1,i2)=vr(4)
      if(nx.lt.1)ui(npot1,i2)=vi(4)
      if(nx.lt.1)ur(npot3,i2)=vr(5)
      if(nx.lt.1)ui(npot3,i2)=vi(5)
      if(nx.lt.1)ur(npot3,npot2)=vr(6)
      if(nx.lt.1)ui(npot3,npot2)=vi(6)
      if(nx.gt.1)zu(npot1,i2)=vr(4)+zi*vi(4)
      if(nx.gt.1)zu(npot3,i2)=vr(5)+zi*vi(5)
      if(nx.gt.1)zu(npot3,npot2)=vr(6)+zi*vi(6)
      go to 422
c     n(15)=3 st up 2 block matrix
 421  zu(npot1,i2)=vr(5)+zi*vi(5)
 422  continue
      if(i1.eq.i2) go to 348
c      off diagnol block submatrix in lower lh side,fill its
c       non diagnol elements via symmetry
      if(nifty(15).eq.3) go to 423
      if(nx.lt.1)ur(npot2,i1)=vr(4)
      if(nx.lt.1)ui(npot2,i1)=vi(4)
      if(nx.lt.1)ur(npot4,i1)=vr(5)
      if(nx.lt.1)ui(npot4,i1)=vi(5)
      if(nx.lt.1)ur(npot4,npot1)=vr(6)
      if(nx.lt.1)ui(npot4,npot1)=vi(6)
      if(nx.gt.1)zu(npot2,i1)=vr(4)+zi*vi(4)
      if(nx.gt.1)zu(npot4,i1)=vr(5)+zi*vi(5)
      if(nx.gt.1)zu(npot4,npot1)=vr(6)+zi*vi(6)
      go to 348
 423  if(nx.gt.1)zu(npot2,i1)=vr(5)+zi*vi(5)
c      block diagnol matrices   (only lower part as rest to follow)
 348  if(nx.lt.1)ur(i1,i2)=vr(1)
      if(nx.lt.1)ui(i1,i2)=vi(1)
      if(nx.lt.1)ur(npot1,npot2)=vr(2)
      if(nx.lt.1)ui(npot1,npot2)=vi(2)
      if(nx.lt.1)ur(npot3,npot4)=vr(3)
      if(nx.lt.1)ui(npot3,npot4)=vi(3)
      if(nx.gt.1)zu(i1,i2)=vr(1)+zi*vi(1)
      if(nifty(15).eq.3) go to 424
      if(nx.gt.1)zu(npot1,npot2)=vr(2)+zi*vi(2)
      if(nx.gt.1)zu(npot3,npot4)=vr(3)+zi*vi(3)
      go to 346
 424  if(nx.gt.1)zu(npot1,npot2)=vr(3)+zi*vi(3)
 346  continue
 347      continue
      istart=ng3
      if(nifty(15).eq.1.or.nifty(15).eq.3) istart=ng2
c      build up upper rh part of super matrix via symmetry

      do 349 i2=2,istart
       i22=i2-1
       do 349 i1=1,i22
      if(nx.gt.1)zu(i1,i2)=zu(i2,i1)
      if(nx.lt.1)ur(i1,i2)=ur(i2,i1)
 349  if(nx.lt.1)ui(i1,i2)=ui(i2,i1)
c      redefine on shell v for channel 2 (others mult by 0 in main)
c      only second row of block matrices changed
      if(nint.eq.1.or.nifty(19).eq.0)go to 245
      if(nifty(15).eq.3) go to 245
      write(6,244)
 244  format(' define on shell v for channnel2 ala e')
      istart=n1+1
      iend=ng2
      iend2=1
      zk=kk(n1)
      if(nifty(19).eq.2)zk=zk2
      np(1)=n1
      np(2)=ng2
      np(3)=ng3
      do 242 ii=istart,iend
      idum=ii-n1
      zkp=kk(idum)
      if(ii.eq.iend) iend2=n1
      do 242 ij=1,iend2
      if(ii.ne.iend) go to 243
c       last lint, fill in entire n1 th row for all block 2
      zk=kk(ij)
      if(ij.eq.n1.and.nifty(19).eq.2)zk=zk2
      np(1)=ij
      np(2)=ij+n1
      np(3)=ij+ng2
  243 call vkbarp(zkp,zk)
      do 242 i=1,3
      jjj=np(i)
      if(nx.lt.1)ur(ii,jjj)=vkp(2,i)
      if(nx.lt.1)ui(ii,jjj)=-zi*vkp(2,i)
      if(nx.gt.1)zu(ii,jjj)=vkp(2,i)
      if(nint.eq.2.and.nx.lt.1)
     1write(6,246)ii,jjj,zkp,zk,ur(ii,jjj),ui(ii,jjj)
  246 format(' i,j,zkp,zk,u=',2i5,6e15.6)
  242 continue
 245  continue
      return
c
 350  format (30h c13 neutron spin ff for k=kp=,e15.5/10(10x,6e15.6/))
 360  format (2i5,2x,i1,i2,8i1,10i3,f10.3,e13.6)
 370  format (8f10.4)
 380  format (5e16.7)
 390  format (2i5/6e16.7)
 400  format (3x,3e19.4,2e11.3,e19.4,2e11.3)
 410  format(28h0opt poten for k =k0(k1-bs)=,e18.5,13hand l=0(1-bs)/
     1     20x,2hkp,
     2 13x,3hrev,16x,2hvc,15h rabs(n0 flip) ,11h rabs(flip),17x,
     3 4h-imv,27h imvabs(noflip),imabs(flip)/)
      end
c/* @(#)psibs.f 1.1      latest revision 7/30/86 14:21:41 */
c/* %w%      latest revision %g% %u% */
      subroutine psibs(zsip,npt,nn,ldum,irmax,deltar,ptptw,pt,wt)
c      normalize p space wf and calculate r space wf via kwon-tabakin
      implicit real*8 (a-h,o-y)
      complex*16 zsip,psir,sump,sumr,cnorm
      dimension ptptw(npt),sbl(50),pt(npt),wt(npt),zsip(nn)
      data pi,hbarc/3.14159265,197.3289/
c     most of this sub follows kwon-tabakin 's sub wavfn
c
c      this portion normalizes the p-space w.f.
c
      sump =0.
c      nb psi ... are complex*16 numbers in this sub
      do 10 i=1,npt
      ptptw(i)=wt(i)*(pt(i)**2)
      sump=sump+zsip(i)*zsip(i)*ptptw(i)
c      nb this renom mixes real*8 and im parts of wf and can give nodes in
 10   continue
c     cnorm=1./csqrt(sump)
      cnorm=1./sqrt(sump)
      write(6,100) cnorm
      do 20 i=1,npt
      zsip(i)=cnorm*zsip(i)
 20   continue
      write(6,110)
      write(6,120) (pt(i),zsip(i),i=1,npt)
      write(6,130)
c
c     r space w.f.
c     irmax= no rvalues wanted     deltar= increment in r
      fac=sqrt(2./pi)
      do 40 i=1,irmax
      sumr=0.
      rval=deltar*(i-1)
      if(rval.eq.0.)rval=1.e-8
      do 30 j=1,npt
      pr=pt(j)*rval/hbarc
c     call spbesl(ldum,pr,sbl)
      sbl(ldum)=sin(pr)/pr
c    ********special l=0 only
      sumr=sumr+sbl(ldum)*zsip(j)*ptptw(j)
 30   continue
      psir=fac*sumr
      write(6,140) rval,psir
 40   continue
 100  format(' cnorm <complex*16> =',1pe14.6,e14.6,/)
 110  format(6x,1hp,10x,8hre(zsip),12x,8him(psip))
 120  format(6(1pe14.6,2x))
 130  format(/,' r (fermi)     re(psir)         im(psir)')
 140  format(0pf8.2,2x,1pe14.6,2x,e14.6)
      return
      end
      subroutine qform(m,n,q,ldq,wa)
c***begin prologue  qform
c
c***changed to double precision
c
c***refer to  snsqe,snsq
c     **********
c
c     subroutine qform
c
c     this subroutine proceeds from the computed qr factorization of
c     an m by n matrix a to accumulate the m by m orthogonal matrix
c     q from its factored form.
c
c     the subroutine statement is
c
c       subroutine qform(m,n,q,ldq,wa)
c
c     where
c
c       m is a positive integer input variable set to the number
c         of rows of a and the order of q.
c
c       n is a positive integer input variable set to the number
c         of columns of a.
c
c       q is an m by m array. on input the full lower trapezoid in
c         the first min(m,n) columns of q contains the factored form.
c         on output q has been accumulated into a square matrix.
c
c       ldq is a positive integer input variable not less than m
c         which specifies the leading dimension of the array q.
c
c       wa is a work array of length m.
c
c     subprograms called
c
c       fortran-supplied ... min0
c
c     minpack. version of january 1979.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
c
c     zero out upper triangle of q in the first min(m,n) columns.
c
c***routines called  (none)
c***end prologue  qform
      implicit real*8 (a-h,o-z)
      integer m,n,ldq
      real*8 q(ldq,m),wa(m)
      integer i,j,jm1,k,l,minmn,np1
      real*8 one,sum,temp,zero
      data one,zero /1.0d0,0.0d0/
c***first executable statement  qform
      minmn = min0(m,n)
      if (minmn .lt. 2) go to 30
      do 20 j = 2, minmn
	 jm1 = j - 1
	 do 10 i = 1, jm1
	    q(i,j) = zero
   10       continue
   20    continue
   30 continue
c
c     initialize remaining columns to those of the identity matrix.
c
      np1 = n + 1
      if (m .lt. np1) go to 60
      do 50 j = np1, m
	 do 40 i = 1, m
	    q(i,j) = zero
   40       continue
	 q(j,j) = one
   50    continue
   60 continue
c
c     accumulate q from its factored form.
c
      do 120 l = 1, minmn
	 k = minmn - l + 1
	 do 70 i = k, m
	    wa(i) = q(i,k)
	    q(i,k) = zero
   70       continue
	 q(k,k) = one
	 if (wa(k) .eq. zero) go to 110
	 do 100 j = k, m
	    sum = zero
	    do 80 i = k, m
	       sum = sum + q(i,j)*wa(i)
   80          continue
	    temp = sum/wa(k)
	    do 90 i = k, m
	       q(i,j) = q(i,j) - temp*wa(i)
   90          continue
  100       continue
  110    continue
  120    continue
      return
c
c     last card of subroutine qform.
c
      end
c/* @(#)qlfun.f 1.1      latest revision 7/30/86 14:21:51 */
c/* %w%      latest revision %g% %u% */
      subroutine qlfun(l,z,ql)
c
c     legrendre functn of second kind for all l values
c     from tabakin s bopit
      implicit real*8 (a-h,o-z)
      dimension fac(15),ql(15)
      if(z.gt.1) go to 30
c     rhl mod use series for z gt 1 vs 10 as tabakin
      if(z.eq.1)write(6,900)
      if(z.eq.1)return
 900  format('qlfun called for z=1')
      ql(1)=0.5*dlog(abs((z+1)/(z-1)))
      ql(2)=z*ql(1)-1.
      if(l.le.2) go to 20
      do 10 n=2,l
      ql(n+1)=(2*n-1)*z*ql(n)-(n-1)*ql(n-1)
 10   ql(n+1)=ql(n+1)/n
 20   return
c
c     series expansion
c
 30   fac(1)=1.
      nmax=l+2
      zsqinv=1/z/z
      do 50 n=1,nmax
 50   fac(n+1)=fac(n)*n/(2.*n+1)
       do 70 lp1=1,l
      ll=lp1-1
      prefac=fac(ll+1)/(z**(ll+1))
      fk=1.
      alfa=0.5*ll
      beta=alfa-.5
      gama=ll+.5
      sum=1.
      do 60 k=1,100
      fk=fk*zsqinv*(alfa+k)*(beta+k)/(gama+k)/k
      if(fk.lt.1.e-17) go to 70
 60   sum=sum+fk
 70   ql(lp1)= prefac*sum
      return
      end
      subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,sigma,acnorm,wa)
c***begin prologue  qrfac
c
c***changed to double precision
c
c***refer to  snlse,snls,snsqe,snsq
c     **********
c
c     subroutine qrfac
c
c     this subroutine uses householder transformations with column
c     pivoting (optional) to compute a qr factorization of the
c     m by n matrix a. that is, qrfac determines an orthogonal
c     matrix q, a permutation matrix p, and an upper trapezoidal
c     matrix r with diagonal elements of nonincreasing magnitude,
c     such that a*p = q*r. the householder transformation for
c     column k, k = 1,2,...,min(m,n), is of the form
c
c                           t
c           i - (1/u(k))*u*u
c
c     where u has zeros in the first k-1 positions. the form of
c     this transformation and the method of pivoting first
c     appeared in the corresponding linpack subroutine.
c
c     the subroutine statement is
c
c       subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,sigma,acnorm,wa)
c
c     where
c
c       m is a positive integer input variable set to the number
c         of rows of a.
c
c       n is a positive integer input variable set to the number
c         of columns of a.
c
c       a is an m by n array. on input a contains the matrix for
c         which the qr factorization is to be computed. on output
c         the strict upper trapezoidal part of a contains the strict
c         upper trapezoidal part of r, and the lower trapezoidal
c         part of a contains a factored form of q (the non-trivial
c         elements of the u vectors described above).
c
c       lda is a positive integer input variable not less than m
c         which specifies the leading dimension of the array a.
c
c       pivot is a logical input variable. if pivot is set true,
c         then column pivoting is enforced. if pivot is set false,
c         then no column pivoting is done.
c
c       ipvt is an integer output array of length lipvt. ipvt
c         defines the permutation matrix p such that a*p = q*r.
c         column j of p is column ipvt(j) of the identity matrix.
c         if pivot is false, ipvt is not referenced.
c
c       lipvt is a positive integer input variable. if pivot is false,
c         then lipvt may be as small as 1. if pivot is true, then
c         lipvt must be at least n.
c
c       sigma is an output array of length n which contains the
c         diagonal elements of r.
c
c       acnorm is an output array of length n which contains the
c         norms of the corresponding columns of the input matrix a.
c         if this information is not needed, then acnorm can coincide
c         with sigma.
c
c       wa is a work array of length n. if pivot is false, then wa
c         can coincide with sigma.
c
c     subprograms called
c
c       minpack-supplied ... r1mach,enorm
c
c       fortran-supplied ... dmax1,sqrt,min0
c
c     minpack. version of december 1978.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
c
c     epsmch is the machine precision.
c
c***routines called  r1mach,enorm
c***end prologue  qrfac
      implicit real*8 (a-h,o-z)
      integer m,n,lda,lipvt
      integer ipvt(lipvt)
      logical pivot
      real*8 a(lda,n),sigma(n),acnorm(n),wa(n)
      integer i,j,jp1,k,kmax,minmn
      real*8 ajnorm,epsmch,one,p05,sum,temp,zero
      real*8 r1mach,enorm
      data one,p05,zero /1.0d0,5.0d-2,0.0d0/
c***first executable statement  qrfac
      epsmch = r1mach(4)
c
c     compute the initial column norms and initialize several arrays.
c
      do 10 j = 1, n
	 acnorm(j) = enorm(m,a(1,j))
	 sigma(j) = acnorm(j)
	 wa(j) = sigma(j)
	 if (pivot) ipvt(j) = j
   10    continue
c
c     reduce a to r with householder transformations.
c
      minmn = min0(m,n)
      do 110 j = 1, minmn
	 if (.not.pivot) go to 40
c
c        bring the column of largest norm into the pivot position.
c
	 kmax = j
	 do 20 k = j, n
	    if (sigma(k) .gt. sigma(kmax)) kmax = k
   20       continue
	 if (kmax .eq. j) go to 40
	 do 30 i = 1, m
	    temp = a(i,j)
	    a(i,j) = a(i,kmax)
	    a(i,kmax) = temp
   30       continue
	 sigma(kmax) = sigma(j)
	 wa(kmax) = wa(j)
	 k = ipvt(j)
	 ipvt(j) = ipvt(kmax)
	 ipvt(kmax) = k
   40    continue
c
c        compute the householder transformation to reduce the
c        j-th column of a to a multiple of the j-th unit vector.
c
	 ajnorm = enorm(m-j+1,a(j,j))
	 if (ajnorm .eq. zero) go to 100
	 if (a(j,j) .lt. zero) ajnorm = -ajnorm
	 do 50 i = j, m
	    a(i,j) = a(i,j)/ajnorm
   50       continue
	 a(j,j) = a(j,j) + one
c
c        apply the transformation to the remaining columns
c        and update the norms.
c
	 jp1 = j + 1
	 if (n .lt. jp1) go to 100
	 do 90 k = jp1, n
	    sum = zero
	    do 60 i = j, m
	       sum = sum + a(i,j)*a(i,k)
   60          continue
	    temp = sum/a(j,j)
	    do 70 i = j, m
	       a(i,k) = a(i,k) - temp*a(i,j)
   70          continue
	    if (.not.pivot .or. sigma(k) .eq. zero) go to 80
	    temp = a(j,k)/sigma(k)
	    sigma(k) = sigma(k)*sqrt(dmax1(zero,one-temp**2))
	    if (p05*(sigma(k)/wa(k))**2 .gt. epsmch) go to 80
	    sigma(k) = enorm(m-j,a(jp1,k))
	    wa(k) = sigma(k)
   80       continue
   90       continue
  100    continue
	 sigma(j) = -ajnorm
  110    continue
      return
c
c     last card of subroutine qrfac.
c
      end
      subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sigma,wa)
c***begin prologue  qrsolv
c***refer to  snlse,snls
c     **********
c
c     subroutine qrsolv
c
c     given an m by n matrix a, an n by n diagonal matrix d,
c     and an m-vector b, the problem is to determine an x which
c     solves the system
c
c           a*x = b ,     d*x = 0 ,
c
c     in the least squares sense.
c
c     this subroutine completes the solution of the problem
c     if it is provided with the necessary information from the
c     qr factorization, with column pivoting, of a. that is, if
c     a*p = q*r, where p is a permutation matrix, q has orthogonal
c     columns, and r is an upper triangular matrix with diagonal
c     elements of nonincreasing magnitude, then qrsolv expects
c     the full upper triangle of r, the permutation matrix p,
c     and the first n components of (q transpose)*b. the syste