c complex refractive index spheroid
c
c gsub.f: subroutines of g.f
c 13 jan 1992: separated subroutines to work on driver
c 12 oct: loop over size possible
c 14 apr: dimension of rr corrected in plegen
c 7 apr: n series convergence assessment
c 3 apr: iconv=-1 for nose on; removed redundant printout
c 28 mar: relaxed cbowkm convergence by 2 and altered printout at failure
c developed subroutine for single particle output (sinfin)
c introduced tabs, cleaning up presentation
c f.f
c 26feb developed subroutine sinsph (single particle results)
c 26feb change to lower case; better loop structure
c 25 feb incorporate long wavelength limit approximation (see lw.f)
c 24 feb new input format better suited to wavelength + refractive index loop
c
c change epslon and related items to common
c data epslon/.5d-15/
c also changed -i**(-j) to 1/(-i)**j where i is complex (mips machines)
c
c modified back for sparc and mips machines, feb 1992
c modified december 1987 to run on cray x-mp: f.cst
cm* are comments added in the process of adding complex refractive index
c 
      block data
      implicit real*8 (a-h,o-z)
c these are the named common blocks that do not appear in main program
      common/coefot/dn(40,100),rn(100),d(100),drho(100),ap(100),dneg
     $     (40),an(20),jupn(40),jdn,jup
      complex*16 uvo3,xyo3,uvi1,xyi1
      common/uvxy/uvo3(40,40),xyo3(40,40),uvo1(40,40),xyo1(40,40),
     $     uvi1(40,40),xyi1(40,40)
      complex*16 di,rni
      common/coefin/rni(100),di(100),jdni,jupi
      complex*16 a,b,alpha,beta
      common/matrix/a(80,80),b(80),alpha(40),beta(40),matmax
      complex*16 t11,t12,t22,t21
      common/ampfun/tsigma(40,10),tchi(40,10),t11(36,7),t12(36,7),t22
     $     (36,7),t21(36,7)
      complex*16 egvi,ri1,rip1,xlami
      common/rwaves/egv,r1,rp1,r2,rp2,egvi,ri1,rip1,xlam(40),xlami(40)
      complex*16 cexte1,cexte2,cextm3,cextm4,cscte2,csctm4
      common/xsectn/ cexte1,cexte2,cextm3,cextm4,cscte2,csctm4
      common/misc/ar2(40),f(40),g(40),chi(40),sigma(40),pie(40,40)
      common/limits/kmax,kmaxr,kmax1,kmaxmm,kmaxr4,m,n,nm
      logical test,testi,comp,easy,out
      common/logic/test,testi,comp,easy,out
      complex*16 cexte,cextm,cscate,cscatm
      common/outpu/cexte,cextm,cscate,cscatm,iconv
      end
c 
c
      subroutine sinfin(aoverb)
c write output from single particle calculation
      implicit real*8 (a-h,o-z)
      complex*16 cexte,cextm,cscate,cscatm
      common/outpu/cexte,cextm,cscate,cscatm,iconv
      complex*16 cfindx
      common/param/cfindx,rfindx,anglin,y,lsig
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
c
      pi=3.141592653589793d0
      cee=cexte
      cme=cextm
      cep=(0.0d0,-1.0d0)*cexte
      cmp=(0.0d0,-1.0d0)*cextm
      write(6,500)
 500  format(1x,t16,'area*k**2',2x, 'qexte',6x, 'qextm',6x,
     $     'tau',8x,'polar',6x, 'qpte',7x, 'qptm',7x,'birefr',5x,
     $     'qscte',6x, 'qsctm',6x,'conv=1')
c
c basic information before normalization
c
      sumext=(cextm+cexte)/2.d0
      pol=(cextm-cexte)/2.d0
      phasdf=dcmplx(0.0d0,-1.0d0)*(cextm-cexte)/2.d0
c
c scaled areas (have a k**2 factor, as do the cross-sections)
      ab=aoverb**2
      if(sgn.eq.-1.0d0) ab=1.0d0/ab
      a=sgn*pi*cout*cout/(ab-1.0d0)
      cos2=dcos(anglin)**2
      sin2=dsin(anglin)**2
      proj1=cos2+sin2/ab
      proj2=(sin2*cos2*(1.0d0/ab-1)**2)/(sin2+cos2/ab)
      pa=a/dsqrt(proj1)
      pga=a/dsqrt(proj1-proj2)
      pab=pi*cout*cout*aoverb/(aoverb**2-1.0d0)
c
c scaled volume
c need to input approariate parameters if writing this
c      if(sgn.eq.1.0d0) then
c          volume=4.0d0*pi*scalea*scaleb**2/3.0d0
c          else
c             volume=4.0d0*pi*scaleb*scalea**2/3.0d0
c          endif
c       write(6,505) volume
c505   format(' ',7hvolume=,d10.4,5h*w**3)
c
c choose edge-on (scaled) area
c
      qtepab=cexte/pab
      qtmpab=cextm/pab
      qpepab=dimag(cexte)/pab
      qpmpab=dimag(cextm)/pab
      qscate=cscate/pab
      qscatm=cscatm/pab
      extpab=sumext/pab
      polpab=pol/pab
      birpab=phasdf/pab
      brdcab=birpab*polpab
      write(6,504)pab,qtepab,qtmpab,extpab,polpab,
     $     qpepab,qpmpab,birpab,qscate,qscatm,iconv
 504  format(1x,'side-on',t14,10(1x,1pd10.3),3x,i2)
c     
      degpol=pol/sumext
      birdic=phasdf/pol
      albedo=(cscate+cscatm)/2.d0/sumext
      write(6,*) ' degree of polarization, birefr/polar, albedo'
      write(6,506) degpol,birdic,albedo
 506  format(1x,3(1x,d11.4))
      return
      end
c
      subroutine sinsph(needed,mold,nmax,mmax,mbot,mtop,kmaxp,cnvrgn)
c single spheroid subroutine
      implicit real*8 (a-h,o-z)
      complex*16 cexte,cextm,cscate,cscatm,cexold,cmxold,cesold,cmsold
      common/outpu/cexte,cextm,cscate,cscatm,iconv
      complex*16 cexte1,cexte2,cextm3,cextm4,cscte2,csctm4
      common/xsectn/ cexte1,cexte2,cextm3,cextm4,cscte2,csctm4
      real*8 jbs,nbs
      complex*16 jbsi,nbsi
      common/spfnct/p(120,22),q(120,22),pl(120),jbs(120),nbs(120),
     $     jbsi(120),nbsi(120),qneg(22),maxdeg,maxord
      dimension qh(120,22)
      equivalence (qh(1,1),q(1,1))
      common/misc/ar2(40),f(40),g(40),chi(40),sigma(40),pie(40,40)
      complex*16 dconjg
      complex*16 cfindx
      common/param/cfindx,rfindx,anglin,y,lsig
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
      common/limits/kmax,kmaxr,kmax1,kmaxmm,kmaxr4,m,n,nm
      logical made
      logical needed
      logical wanted
      logical test,testi,comp,easy,out
      common/logic/test,testi,comp,easy,out
      common/output/ wanted
c     
      pi=3.141592653589793d0
      if(mmax.gt.mold)made=.false.
      mold=max0(mold,mmax)
c
cm* outside incident wave
      if(zeta.ne.0.0d0)call spbess(nmax,jbs,nbs)
c inside
cm*
      if(zeta.ne.0.0d0)call cspbes(nmax,jbsi,nbsi)
c outside again
      c=cout
      if(sgn.eq.-1.0d0) goto 2
cm* outside prolate
      comp=.true.
      easy=.false.
      if(zeta.le.1.30d0) comp=.true.
      if(zeta.gt.1.050d0.and.zeta.le.1.30d0) easy=.true.
      if(zeta.gt.1.30d0.and.zeta.le.2.0d0.and.c.ge.10.0d0) comp=.false.
      if(zeta.gt.2.0d0.and.zeta.le.5.0d0.and.c.ge.1.0d0)comp=.false.
      if(zeta.gt.5.0d0.and.c.ge.1.0d-2) comp=.false.
      if(comp.and..not.made)call plegen(zeta,y,mmax,nmax,lsig,p,maxdeg,
     $     maxord)
      if(comp.and..not.made)call qlegen(zeta,y,mmax,nmax,lsig,q,maxdeg,
     $     maxord)
      if(comp)made=.true.
      goto 3
 2    continue
cm* outside oblate
      if(zeta.ne.0.0d0.and..not.made)call qhyper(zeta,mmax,nmax,qh,maxd
     $     eg,maxord)
      made=.true.
      comp=.false.
 3    continue
c
cm* main loop
c set for nose on
      iconv=-1
      cexte=(0.0d0,0.0d0)
      cextm=(0.0d0,0.0d0)
      cscate=(0.0d0,0.0d0)
      cscatm=(0.0d0,0.0d0)
      do 20 mm=mbot,mtop
         m=mm-1
         if(anglin.ne.0.0d0)call elegen(nmax,m,anglin,pl)
         if(m.ne.0.and.comp)call qlgneg(zeta,m,q(1,m+1),q(2,m+1),qneg)
c test,testi are set are in old version
c out is true if wavefunctions underflow which cuts down kmax
         out=.false.
c      series run from m,m+1,..,n such that m+1 and n are of the same
c      parity,define kmax=n-m which implies kmax is odd
         kmax=kmaxp
         kmaxr=(kmax+1)/2
         kmax1=kmax+1
         kmaxmm=kmax+mm
         kmaxr4=4*kmaxr
         cexte1=(0.0d0,0.0d0)
         cexte2=(0.0d0,0.0d0)
         cextm3=(0.0d0,0.0d0)
         cextm4=(0.0d0,0.0d0)
         cscte2=(0.0d0,0.0d0)
         csctm4=(0.0d0,0.0d0)
c make vectors of eigenvalues in,out-side
         call egnvec
         if(wanted) write(6,666)
 666     format(
     $        3x,4hm  n,t21,2hr1,t44,3hrp1,t67,2hr2,t91,3hrp2,
     $        t116,3hegv,t130,3hacc/
     $        t19,5hsigma,t44,3hchi,      t67,4hegvi
     $        ,t130,3hjup/
     $        t19,3hri1,                  t67,4hrip1)
cm* main inner loop
         do 30 nn=mm,kmaxmm
            if(out) goto 30
            n=nn-1
            call radwav
            call capfun
 30      continue
cm*
         call scater
         if(needed) call tangle
         call temod1
         call temod2
         if(needed) call teamp
         call tmmod3
         call tmmod4
         if(needed) call tmamp
c since (cos m*a,cos n*a)=2*pi if m=n=0 scattering x-section must be
c doubled for m=0
         if(m.eq.0) cscte2=2.0d0*cscte2
         if(m.eq.0) csctm4=2.0d0*csctm4
c exp(i*(k*r-w*t)) convention is used -im(cext) gives proper phases
         cexte=-4.0d0*pi*dconjg(cexte1+cexte2)+cexte
         cextm=-4.0d0*pi*dconjg(cextm3+cextm4)+cextm
         cscate=pi*cscte2+cscate
         cscatm=pi*csctm4+cscatm
c
c no convergence check for nose on case (m=1 so no intermediate output each m)
         if(m.eq.1.and.anglin.eq.0.0d0) goto 23
c
         if(wanted) write(6,900) cexte,cscate
         if(wanted) write(6,901) cextm,cscatm
 900     format(1x,20hte mode  cexte*k**2=,2d24.16,2x,
     $        12hcscate*k**2=,2d24.16)
 901     format(1x,20htm mode  cextm*k**2=,2d24.16,2x,
     $        12hcscatm*k**2=,2d24.16/)
c
c no convergence check possible when m=0
         if(m.eq.0) goto 19
c
cm* oblique incidence: check for convergence w.r.t. mm (given kmaxp)
c
cm* check for convergence
c convergence test for extinction and phases
         cxten=cexte
         cpten=(0.0d0,-1.0d0)*cexte
         cxtmn=cextm
         cptmn=(0.0d0,-1.0d0)*cextm
         cnvr1=(cexte-cexold)/cxten
         cnvr2=(cextm-cmxold)/cxtmn
         cnvr3=(0.0d0,-1.0d0)*(cexte-cexold)/cpten
         cnvr4=(0.0d0,-1.0d0)*(cextm-cmxold)/cptmn
c original m,n series convergence test is applicable only in the case of
c real refractive index, in the use of absorption = cscatt
cm* remove these two which are only valid for real refractive index
cm*       cnvrge=(cexte-cscate)/cxten
cm*       cnvrgm=(cextm-cscatm)/cxtmn
cm* introduce the old scattering comparison (real part)
         cscten = cscate
         csctmn = cscatm
         cnvr5 = (cscate-cesold)/cscten
         cnvr6 = (cscatm-cmsold)/csctmn
cm* note that cabs = cextinction - cscattering
cm* will not be converged to the same accuracy by these
cm* criteria, especially if small
         cnvrp=dmax1(dabs(cnvr1),dabs(cnvr2),dabs(cnvr3),
     $        dabs(cnvr4),dabs(cnvr5),dabs(cnvr6) )
         if(cnvrp.ge.cnvrgn) goto 19
         write(6,800)m
 800     format(1x,'given n series, m series has converged at m=',i2)
         write(6,900) cexte,cscate
         write(6,901) cextm,cscatm
         iconv=1
c converged; skip out of loop with flag iconv set
         goto 21
cm* prepare for next loop
 19      continue
c if not converged, do not update ...old
         if(mm.eq.mtop) goto 22
         cexold=cexte
         cmxold=cextm
         cesold=cscate
         cmsold=cscatm
 20   continue
c end of main loop
c
 22   continue
c not converged by end of loop; write last two iterations; set flag iconv
      iconv=0
      write(6,801)m
 801  format(1x,'given n series, m series not converged by m=',i2)
      write(6,900) cexold,cesold
      write(6,901) cmxold,cmsold
 23   continue
c nose on printout
      write(6,900) cexte,cscate
      write(6,901) cextm,cscatm
c
c converged hits here out of loop
 21   continue
c
c output available to main program in /outpu/
c
      if(needed) call trans
c transformation matrix has self-contained output
c
      return
      end
c
      subroutine spbess(nmax,jbs,nbs)
c     spherical bessel and neumann functions stores as jbs(m+1),nbs(m+1)
cm* jbs by gautschi's algorithm, except for czeta .gt. nmax
      implicit real*8 (a-h,o-z)
      common /accura/ epslon,epslo2
      real*8 jbs(nmax),nbs(nmax)
      logical wanted
      common/output/ wanted
      logical flag
      dimension f(120),rr(120)
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
cm* c is real
      nflag=nmax+1
      jflag=nflag
      c = cout
      nmax1=nmax-1
      pi=3.141592653589793d0
      czeta=c*zeta
      czeta2=czeta/2.0d0
      gamma=dsqrt(pi)
      p=dsqrt(pi/2.0d0/czeta)
 100  if(idint(czeta).gt.nmax) goto 300
c
cm* prepare for gautschi
      do 110 j=1,nmax
 110  f(j)=0.0d0
      cgamma=.50d0*gamma
      sum=dsqrt(czeta2)/cgamma
cm* get starting nu
      d1=dlog(10.0d0)*16.0d0+dlog(4.0d0)
      r=0.0d0
      if(nmax.gt.1)r=nmax*tinvrt(.50d0*d1/dfloat(nmax))
      s=1.35910d0*czeta*tinvrt(.735760d0*d1/czeta)
      nu=1+idint(dmax1(r,s))
cm* gautschi loop
 115  m=0
      el=1.0d0
      lim=nu/2
 120  m=m+1
      el=el*(dfloat(m)+.50d0)/dfloat(m+1)
      if(m.lt.lim)goto 120
      n=2*m
      r=0.0d0
      s=0.0d0
 125  r=1.0d0/(2.0d0*(.50d0+dfloat(n))/czeta-r)
      elam=0.0d0
      if(mod(n,2).ne.0) goto 130
      el=el*dfloat(n+2)/dfloat(n+1)
      elam=el*(dfloat(n)+.50d0)
 130  s=r*(elam+s)
      if(n.le.nmax1) rr(n)=r
      n=n-1
      if(n.ge.1) goto 125
      jbs(1)=sum/(1.0d0+s)*p
      flag=.true.
      do 136 j=1,nmax1
         if(jbs(j).eq.0.0d0) goto 135
         if(dlog(dabs(jbs(j)))+dlog(dabs(rr(j))).ge.-170.0d0)
     $        goto 135
         jbs(j+1)=0.0d0
         if(flag)jflag=j
cm* jbs(jflag+1), i.e. order jflag, first to underflow
         flag=.false.
         goto 136
 135     jbs(j+1)=rr(j)*jbs(j)
 136  continue
      do 140 j=1,nmax
         if(f(j).eq.0.0d0.and.jbs(j).eq.0.0d0) goto 140
         if(jbs(j).eq.0.0d0) goto 142
         if(dabs(1.0d0-f(j)/jbs(j)).le.epslon) goto 140
 142     do 145 k=1,nmax
 145     f(k)=jbs(k)
         nu=nu+5
         goto 115
 140  continue
      goto 400
cm* forwards recursion ok when large czeta
 300  jbs(1)=dsin(czeta)/czeta
      jbs(2)=jbs(1)/czeta-dcos(czeta)/czeta
      flag=.true.
      do 350 j=3,nmax
 350  jbs(j)=dfloat(2*(j-2)+1)/czeta*jbs(j-1)-jbs(j-2)
 400  continue
      if(.not.flag.and.wanted) write(6,410) jflag
 410  format(1x,45h***spherical bessel function underflow,order ,i3)
cm*
cm* forwards recursion
      nbs(1)=-dcos(czeta)/czeta
      nbs(2)=nbs(1)/czeta-jbs(1)
      do 450 j=3,nmax
         nflag=j-2
         if(dabs(nbs(j-1)).gt.1.0d70) goto 500
cm* nbs(nflag+1), i.e. order nflag, first to overflow
 450  nbs(j)=dfloat(2*(j-2)+1)/czeta*nbs(j-1)-nbs(j-2)
      nflag=nmax+1
      goto 600      
 500  if(wanted)write(6,510) nflag
 510  format(1x,45h***spherical neumann function overflow,order ,i3)
      do 520 jj=nflag+2,nmax
 520  nbs(jj)=1.0d70
 600  continue
cm* wronskian check
      nup = min0(nmax,jflag,nflag)
      w=1./(czeta*czeta)
      feps = 50.*epslon
      do 601 j=1,nup-1
         ww = jbs(j+1)*nbs(j)-jbs(j)*nbs(j+1)
         wd = dabs((w-ww)/w)
         if(wd.gt.feps) write(6,602) j,w,ww,wd
 601  continue
 602  format(1x,'b wr: spbess',1x,i3,6e10.3)
      if(.not.wanted) return
cm* printout
      write(6,605)
 605  format(1x,26hspherical bessel functions)
      do 604 j=1,min0(nmax,nup+1)
         write(6,603) j,jbs(j),nbs(j)
 604  continue
 603  format(1x,i3,6e11.4)
      return
      end
c
      function tinvrt(y)
c inverts y=t*ln(t) for spherical bessel function routine's start pt
      implicit real*8(a-h,o-z)
      if(y.gt.10.0d0) goto 50
      p=5.79410d-5*y-1.761480d-3
      p=y*p+2.086450d-2
      p=y*p-.1290130d0
      p=y*p+.857770d0
      tinvrt=y*p+1.01250d0
      return
 50   z=dlog(y)-.7750d0
      p=(.7750d0-dlog(z))/(1.0d0+z)
      p=1.0d0/(1.0d0+p)
      tinvrt=y*p/z
      return
      end
c
      subroutine cspbes(nmax,jbs,nbs)
c spherical bessel and neumann functions stores as jbs(m+1),nbs(m+1)
      implicit real*8 (a-h,o-z)
      common /accura/ epslon,epslo2
cm* jbs by gautschi's algorithm, except for czeta .gt. nmax
cm*
cm* for complex argument
cm* new complex typing
      complex*16 jbs(nmax),nbs(nmax)
      complex*16 c,czeta,czeta2,czero,f,nlast,p,r,rr,s,sum,w,ww
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
      complex*16 cfindx
      common/param/cfindx,rfindx,anglin,y,lsig
cm*
      logical wanted
      common/output/ wanted
      logical flag
      dimension f(120),rr(120)
cm* c is complex
cm* these flags are used as order indices, rather than fortran indices
      nflag=nmax+1
      jflag=nflag
      c = cin
      nmax1=nmax-1
      pi=3.141592653589793d0
      czeta=c*zeta
      czeta2=czeta/2.0d0
      gamma=dsqrt(pi)
      p=cdsqrt(pi/2.0d0/czeta)
cm* following split is a guess which should work for reasonable values 
cm* of imaginary part
      aczeta=cdabs(czeta)
 100  if(int(aczeta).gt.nmax) goto 300
cm* prepare for gautschi
      czero=(0.0,0.0)
      do 110 j=1,nmax
 110  f(j)=czero
      cgamma=.50d0*gamma
      sum=cdsqrt(czeta2)/cgamma
cm* get starting nu
      d1=dlog(10.0d0)*16.0d0+dlog(4.0d0)
      ar=0.0d0
      if(nmax.gt.1)ar=nmax*tinvrt(.50d0*d1/dfloat(nmax))
cm* following a rough guess too; note that nu iterates if not enough
      as=1.35910d0*aczeta*tinvrt(.735760d0*d1/aczeta)
      nu=1+idint(dmax1(ar,as))
cm* gautschi loop
 115  m=0
      el=1.0d0
      lim=nu/2
 120  m=m+1
      el=el*(dfloat(m)+.50d0)/dfloat(m+1)
      if(m.lt.lim)goto 120
      n=2*m
      r=czero
      s=czero
 125  r=1.0d0/(2.0d0*(.50d0+dfloat(n))/czeta-r)
      elam=0.0d0
      if(mod(n,2).ne.0) goto 130
      el=el*dfloat(n+2)/dfloat(n+1)
      elam=el*(dfloat(n)+.50d0)
 130  s=r*(elam+s)
      if(n.le.nmax1) rr(n)=r
      n=n-1
      if(n.ge.1) goto 125
      jbs(1)=sum/(1.0d0+s)*p
      flag=.true.
      do 136 j=1,nmax1
         if(jbs(j).eq.0.0d0) goto 135
cm* underflow check on the modulus
         if(dlog(cdabs(jbs(j)))+dlog(cdabs(rr(j))).ge.-170.d0)
     $        goto 135
         jbs(j+1)=czero
         if(flag)jflag=j
cm* jbs(jflag+1), i.e. order jflag, first to underflow
         flag=.false.
         goto 136
 135     jbs(j+1)=rr(j)*jbs(j)
 136  continue
      do 140 j=1,nmax
         if(f(j).eq.czero.and.jbs(j).eq.czero) goto 140
         if(jbs(j).eq.czero) goto 142
cm* check convergence on both real and imaginary
         rf = dreal(f(j))
         rj = dreal(jbs(j))
         if(rj.eq.0.d0) goto 142
         rfj = 1.d0-rf/rj
         af = dimag(f(j))
         aj = dimag(jbs(j))
cm? hopefuly this will not allow a false convergence, i.e.
cm* to get here rj is already non-zero, and so aj should be
cm* moving away from zero if it is going to be non-zero
         if(aj.eq.0.d0) afj= 0.d0
         if(aj.ne.0.d0) afj = 1.d0-af/aj
cm*   if(dabs(1.0d0-f(j)/jbs(j)).le.epslon) goto 140
         if(dabs(rfj).le.epslon .and. dabs(afj).le.epslon) goto 140
 142     do 145 k=1,nmax
 145     f(k)=jbs(k)
         nu=nu+5
         goto 115
 140  continue
      goto 400
cm* forwards recursion ok when large czeta
 300  jbs(1)=cdsin(czeta)/czeta
      jbs(2)=jbs(1)/czeta-cdcos(czeta)/czeta
      flag=.true.
      do 350 j=3,nmax
 350  jbs(j)=dfloat(2*(j-2)+1)/czeta*jbs(j-1)-jbs(j-2)
 400  continue
      if(.not.flag.and.wanted) write(6,410) jflag
 410  format(1x,45h***spherical bessel function underflow,order ,i3)
cm*
cm* forwards recursion
      nbs(1)=-cdcos(czeta)/czeta
      nbs(2)=nbs(1)/czeta-jbs(1)
      do 450 j=3,nmax
         nflag=j-2
cm* overflow check on the modulus
         if(cdabs(nbs(j-1)).gt.1.0d70) goto 500
cm* nbs(nflag+1), i.e. order nflag, first to overflow
 450  nbs(j)=dfloat(2*(j-2)+1)/czeta*nbs(j-1)-nbs(j-2)
      nflag=nmax+1
      goto 600      
 500  if(wanted)write(6,510) nflag
 510  format(1x,45h***spherical neumann function overflow,order ,i3)
cm* use nlast rather than (1.e70,1.e70)
      nlast = nbs(nflag+1)
      if(wanted) write(6,511) nlast
 511  format(1x,5hnlast,2x,2e10.3)
      do 520 jj=nflag+2,nmax
 520  nbs(jj)= nlast
 600  continue
cm* wronskian check
      nup = min0(nmax,jflag,nflag)
      w=1.d0/(czeta*czeta)
      rw = dreal(w)
      aw = dimag(w)
      if(aw.eq.0.d0) awww = 0.d0
      fepfac=50.d0+2000.d0*
     $     dsqrt(dabs(dimag(cfindx)/(dreal(cfindx)-1.d0)))
cm?? remove this 3.e3 limit only after testing
      fepf = dmin1(fepfac,3.d3)
      feps = epslon*fepf
      do 601 j=1,nup-1
         ww = jbs(j+1)*nbs(j)-jbs(j)*nbs(j+1)
cm* check real and imaginary parts of wronskian separately
         rww = dreal(ww)
         aww = dimag(ww)
         rwww = dabs(1.d0-rww/rw)
         if(aw.ne.0.) awww = dabs(1.d0-aww/aw)
         if(rwww.gt.feps .or. awww.gt.feps) 
     $        write(6,602) j,w,ww,rwww,awww
 601  continue
 602  format(1x,'b wr: cspbes',1x,i3,6e10.3)
      if(.not.wanted) return
cm* printout
      write(6,605)
 605  format(1x,33hinside spherical bessel functions)
      do 604 j=1,min0(nmax,nup+1)
         write(6,603) j,jbs(j),nbs(j)
 604  continue
 603  format(1x,i3,6e11.4)
      return
      end
c
      subroutine plegen(zeta,y,mmax,nmax,lsig,p,maxdeg,maxord)
c associated legendre polynomials p(n,m) stores as p(n+1,m+1)
      implicit real*8(a-h,o-z)
      logical wanted
      common/output/ wanted
      dimension p(maxdeg,maxord),rr(22)
      do 100 nn=1,nmax
         n=nn-1
         n2=n+2
         if(zeta.ne.1.0d0.and.n.ne.0) goto 20
         p(nn,1)=1.0d0
         do 10 j=2,mmax
 10      p(nn,j)=0.0d0
         goto 100
 20      if(lsig.eq.0) zeta1=dsqrt(zeta*zeta-1.0d0)
         if(lsig.eq.1) zeta1=dsqrt(y*(2.0d0+y))
         sum=(zeta+zeta1)**n
         zeta1=2.0d0*zeta/zeta1
         r=0.0d0
         s=0.0d0
         mmax1=mmax-1
         do 50 j=1,n
            k=n+1-j
            r=(n+1-k)/(k*zeta1+(k+n+1)*r)
            s=r*(2.0d0+s)
            if(k.le.mmax1) rr(k)=r
 50      continue
         p(nn,1)=sum/(1.0d0+s)
         lim=min0(mmax-1,n)
         do 60 j=1,mmax1
            rj=0.0d0
            if(j.le.lim) rj=rr(j)
            if(dabs(p(nn,j)).gt.1.0d70) goto 105
 60      p(nn,j+1)=(j+n)*rj*p(nn,j)
         goto 100
 105     if(wanted) write(6,110) j,n
 110     format(1x,' ***legendre first kind overflow, order ',
     $        i3,2x,7hdegree ,i3)
         do 115 k=j,mmax1
 115     p(nn,k)=1.0d70
 100  continue
      return
      end
c
      subroutine qlegen(zeta,y,mmax,nmax,lsig,q,maxdeg,maxord)
c associated legendre polynomials q(n,m) stores as q(n+1,m+1)
      implicit real*8(a-h,o-z)
      common /accura/ epslon,epslo2
      logical wanted
      common/output/ wanted
      dimension q(maxdeg,maxord),f(120),rr(120)
      if(lsig.eq.0) y=zeta-1.0d0
      if(lsig.eq.0)zeta1=dsqrt(zeta*zeta-1.0d0)
      if(lsig.eq.1) zeta1=dsqrt(y*(y+2.0d0))
      q(1,1)=.50d0*dlog((zeta+1.0d0)/y)
      q(1,2)=-1.0d0/zeta1
      do 100 mm=1,mmax
         m=mm-1
         if(m.le.1) goto 20
         ej=dfloat(m-1)
         q(1,mm)=-ej*2.0d0*zeta*q(1,mm-1)/
     $        zeta1-ej*(ej-1.0d0)*q(1,mm-2)
 20      do 30 j=1,nmax
 30      f(j)=0.0d0
         nmax1=nmax-1
         nu=20+int(1.25*nmax)
 40      r=0.0d0
         do 50 j=1,nu
            k=nu-j+1
            r=dfloat(k+m)/(dfloat(2*k+1)*zeta-dfloat(k-m+1)*r)
            if(k.le.nmax1) rr(k)=r
 50      continue
         do 60 j=1,nmax1
            if(dabs(q(j,mm)).le.1.0d-70) goto 61
 60      q(j+1,mm)=rr(j)*q(j,mm)
         goto 65
 61      if(wanted) write(6,62) m,j
 62      format(1x,40h***legendre second kind underflow order ,
     $        i3,2x,7hdegree ,i3)
         do 63 k=j,nmax1
 63      q(k,mm)=1.0d-70
 65      do 70 j=1,nmax
            if(dabs(1.0d0-f(j)/q(j,mm)).le.epslon) goto 70
            do 80 k=1,nmax
 80         f(k)=q(k,mm)
            nu=nu+10
            goto 40
 70      continue
 100  continue
      return
      end
c
      subroutine qhyper(zeta,mmax,nmax,qh,maxdeg,maxord)
c calculates q(n,m) of i*zeta by hypergeometric series and recursion
c stores as i**(n+2m+1)q(n,m)=qh(n+1,m+1)
      implicit real*8(a-h,o-z)
      common /accura/ epslon,epslo2
      real*8 jr
      logical flag
      dimension qh(maxdeg,maxord)
      nmax1=nmax-1
      nmax2=nmax-2
      zeta1=dsqrt(zeta*zeta+1.0d0)
      z=(zeta+zeta1)
      z2=-1.0d0/z/z
      e2n=2.0d0
c     dact=fact(n)/fact(2n)
      dact=1.0d0
      flag=.true.
      do 10 nn=1,nmax1
         n=nn-1
         f=1.0d0
         a1=.50d0
         b1=dfloat(n+1)
         c1=dfloat(2*n+3)/2.0d0
         f2=a1
         f4=z2*b1/c1
         f1=0.0d0
         jr=1
 15      continue
         f=f2*f4+f
         if(dabs(f-f1).le.epslon*dabs(f)) goto 20
         f2=f2*(a1+jr)/(jr+1.0d0)
         f4=f4*(b1+jr)*z2/(c1+jr)
         f1=f
         jr=jr+1
         goto 15
 20      continue
         if(flag) goto 21
         if(dlog(e2n*dact*f)-(n+1)*dlog(z).ge.-177) goto 21
         qh(nn,1)=0.0d0
         goto 22
 21      qh(nn,1)=e2n*dact*f*z**(-n-1)
 22      if(qh(nn,1).le.1.0d-70)flag=.false.
         if(qh(nn,1).gt.1.0d-70) flag=.true.
         e2n=e2n*2.0d0
         dact=dact*dfloat(nn)/dfloat(2*nn+1)
 10   continue
      mmax1=mmax-1
      do 30 mm=1,mmax1
         m=mm-1
         qh(nmax,mm)=(dfloat(nmax2+m)*qh(nmax1-1,mm)-
     $        dfloat(2*nmax2+1)*zeta*qh(nmax1,mm))/dfloat(nmax2-m+1)
         do 40 nn=1,nmax1
            n=nn-1
            qh(nn,mm+1)=(dfloat(n-m+1)*qh(nn+1,mm)+
     $           dfloat(n+m+1)*zeta*qh(nn,mm))/zeta1
 40      continue
 30   continue
      return
      end
c
      subroutine qlgneg(zeta,m,qq0,qq1,qneg)
c calculates q(-n,m) stores as qneg(n)
      implicit real*8(a-h,o-z)
      dimension qneg(m)
cm* these elements are already a function of zeta
      q0=qq0
      q1=qq1
      if(m.eq.0) return
      do 10 jj=1,m
         j=1-jj
         qneg(jj)=(dfloat(2*j+1)*zeta*q0-
     $        dfloat(j-m+1)*q1)/dfloat(j+m)
         q1=q0
         q0=qneg(jj)
 10   continue
      return
      end
c
      subroutine elegen(nmax,m,theta,pl)
c calculates,by hypergeometric series,associated legendre functions
c for angle of incidence for a given order(m)
c ferrer's definition is used
cm* used in fandg outside only
      implicit real*8 (a-h,o-z)
      common /accura/ epslon,epslo2
      real*8 mult,norm
      logical flag
      dimension pl(nmax)
      xsqr=dcos(theta)
      y=dsin(theta)
      pi=3.141592653589793d0
      if(dabs(theta-pi/2.0d0).le.epslon) xsqr=0.0d0
      if(dabs(theta-pi/2.0d0).le.epslon) y=1.0d0
      if(theta.le.epslon) xsqr=1.0d0
      if(theta.le.epslon) y=0.0d0
      x=xsqr*xsqr
      flag=.false.
      if(m.eq.0.and.y.eq.0.0d0) flag=.true.
      if(.not.flag) ym=y**m
      do 20 nn=1,nmax
         n=nn-1
         pl(nn)=0.0d0
         if(n.lt.m) goto 20
         if(flag)  pl(nn)=1.0d0
         if(flag) goto 20
 2       k=0
         if(mod(n-m,2).ne.0) k=1
         dfact=1.0d0
         lim=n+m-1+k
         if(lim.lt.1) goto 30
         limb=(n-m-k)/2
         do 27 j=1,lim
            dj1=1.0d0
            if(mod(j,2).ne.0) dj1=dfloat(j)
            dj2=1.0d0
            if(j.le.limb) dj2=dfloat(j)
 27      dfact=dfact*dj1/dj2
 30      norm=dfloat((-1)**((n-m-k)/2))*dfact/(2.0d0**((n-m-k)/2))
         if(x.ge.epslo2) goto 5
         if(k.eq.0) pl(nn)=norm
         goto 20
 5       f1=1.0d0
         a=dfloat(m-n+k)/2.0d0
         b=dfloat(m+n+1+k)/2.0d0
         c=.50d0+dfloat(k)
         r=1.0d0
         term=a*b*x/c
 15      f=f1+term
         if(term.eq.0.0d0) goto 35
         mult=(a+r)*(b+r)*x/(c+r)/(r+1.0d0)
         term=term*mult
         r=r+1.0d0
         f1=f
         goto15
 35      pl(nn)=norm*xsqr**k*ym*f
 20   continue
      return
      end
c
      subroutine radwav
      implicit real*8(a-h,o-z)
      real*8 jbs,nbs
      complex*16 jbsi,nbsi
      logical wanted
      logical test,testi,comp,easy,out
      common/coefot/dn(40,100),rn(100),d(100),drho(100),ap(100),dneg
     $     (40),an(20),jupn(40),jdn,jup
      complex*16 di,rni
      common/coefin/rni(100),di(100),jdni,jupi
      common/misc/ar2(40),f(40),g(40),chi(40),sigma(40),pie(40,40)
      common/spfnct/p(120,22),q(120,22),pl(120),jbs(120),nbs(120),
     $     jbsi(120),nbsi(120),qneg(22),maxdeg,maxord
      dimension qh(120,22)
      equivalence (qh(1,1),q(1,1))
      complex*16 cfindx
      common/param/cfindx,rfindx,anglin,y,lsig
      common/logic/test,testi,comp,easy,out
      complex*16 egvi,ri1,rip1,xlami
      common/rwaves/egv,r1,rp1,r2,rp2,egvi,ri1,rip1,xlam(40),xlami(40)
      common/limits/kmax,kmaxr,kmax1,kmaxmm,kmaxr4,m,n,nm
      common/output/ wanted
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
      logical select
      select=.false.
      if(out)return
      nm=n-m
cm* wronskian used to check which evaluation of the
cm* spheroidal wave functions on outside to adopt
      w=1.0d0/cout/(zeta*zeta-sgn)
c 
cm* outside incident
      jj=kmax
      egv=xlam(nm+1)
      call bowkmp(n,m,egv,rn,jj)
      call excoef(n,m,rn,jj,jdn,jup,d,ddlog,dsgn)
      jupn(nm+1)=jup
      do 4 jj=jdn,jup,2
 4    dn(nm+1,jj)=d(jj)
c 
c inside
      jji=kmax
cm*
      egvi=xlami(nm+1)
      call cbowkm(n,m,egvi,rni,jji)
      call cexcoe(n,m,rni,jji,jdni,jupi,di,ddlogi,dsgni)
      call cradia(n,m,jdni,jupi,di,jbsi,jupi+m,ri1,rip1)
c 
cm* outside scattered
      c=cout
      if(sgn.eq.1.0d0) goto 200
      if(zeta.ne.0.0d0) goto 100
cm* oblate zeta = 0, outside scattered
      call xpaneg(n,m,egv,d(2),d(3),dneg,ddlog,dsgn,ddlog2,dsgn2)
      call obzero(n,m,dneg,d,jup,r1,rp1,r2,rp2)
      goto 300
cm* oblate outside scattered
 100  call radial(n,m,jdn,jup,d,jbs,jup+m,r1,rp1)
      call baber(n,m,egv,an,ap,lima)
      call radob(n,m,r1,rp1,qh,maxdeg,maxord,an,ap,lima,r2,rp2)
      if(.not.(zeta.gt.1.0d0.or.c*zeta.gt.10.0d0))goto 300
      call radial(n,m,jdn,jup,d,nbs,jup+m,r2a,rp2a)
      wa=r1*rp2a-r2a*rp1
      wb=r1*rp2-r2*rp1
      if(dabs(w-wa).ge.dabs(w-wb))       goto 300
      r2=r2a
      rp2=rp2a
      goto 300
cm* prolate outside scattered
 200  call radial(n,m,jdn,jup,d,jbs,jup+m,r1,rp1)
      if(comp)call xpaneg(n,m,egv,d(2),d(3),dneg,ddlog,dsgn,ddlog2,dsgn
     $     2)
      if(m.ne.0.and.comp)call dsubro(n,m,egv,dneg(2*m),dneg(2*m-1),drho
     $     ,ldrho)
      if(m.eq.0.and.comp)call dsubro(n,m,egv,d(2),d(3),drho,ldrho)
      if(.not.(comp.and.easy)) goto 250
      call radpr(n,m,dneg,drho,ldrho,d,jup,qneg,q,p,maxdeg,maxord,lsig,
     $     y,r2a,rp2a,ddlog,dsgn,ddlog2,dsgn2)
      call radial(n,m,jdn,jup,d,nbs,jup+m,r2b,rp2b)
      wa=r1*rp2a-r2a*rp1
      wb=r1*rp2b-r2b*rp1
      if(dabs(w-wa).le.dabs(w-wb))goto 225
      r2=r2b
      rp2=rp2b
      goto 300
 225  r2=r2a
      rp2=rp2a
      select=.true.
      goto 300
 250  if(.not.comp)call radial(n,m,jdn,jup,d,nbs,jup+m,r2,rp2)
      if(comp)call radpr(n,m,dneg,drho,ldrho,d,jup,qneg,q,p,maxdeg,maxo
     $     rd,lsig,y,r2,rp2,ddlog,dsgn,ddlog2,dsgn2)
      if(comp)select=.true.
c 
cm* main continuation for outside oblate or prolate
 300  continue
      diff=1-(r1*rp2-r2*rp1)/w
      if(diff.ne.0.0d0)iacc=idint(dlog10(dabs(diff)))
      if(diff.eq.0.0d0)iacc=-16
      if(wanted) write(6,301)m,n,r1,rp1,r2,rp2,egv,iacc
 301  format(1x,2i3,5d24.16,i4)
      if(wanted.and.select) write(6,500)
 500  format(1x,11hselect true)
c    sigmas and chi's for cross sections
      call angl(m,jdn,jup,d,anglin,pl,nmax,chie,sigmae)
      chi(nm+1)=chie
      sigma(nm+1)=sigmae
      if(wanted) write(6,777) sigmae,chie,egvi,jup
      if(wanted) write(6,777) ri1,rip1
 777  format(7x,4d24.16,25x,i3)
c scaling factors for matrices and alphas,betas
      ar2(nm+1)=r2
      if(dabs(rp2).lt.dabs(r2)) ar2(nm+1)=rp2
c  f's and g's for constant vector
      call fandg (m,anglin,d,jdn,jup,pl,nmax,fe,ge)
      f(nm+1)=fe
      g(nm+1)=ge
      if(dabs(r1).lt.1.0d70.and.dabs(r1).gt.1.0d-70) return
      write(6,24) n
 24   format(1x,'wave function out of range 1.d+/-70, degree ',i2)
      out=.true.
      kmax=nm
      if(mod(nm,2).eq.0) kmax=nm-1
      kmax1=kmax+1
      kmaxr=kmax1/2
      kmaxr4=kmaxr*4
      return
      end
c
      subroutine egnvec
c vector of eigenvalues found from diagonalizing real symmetric tri-diag
c matrix(tql1).finds egv(m,m),egv(m+1,m)... as xlam(1),xlam(2)... out
c and in -side spheroid.accuracy is improved and d(r,mn) made by bowkmp
      implicit real*8(a-h,o-z)
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
      common/limits/kmax,kmaxr,kmax1,kmaxmm,kmaxr4,m,n,nm
      complex*16 cfindx
      common/param/cfindx,rfindx,anglin,y,lsig
      complex*16 egvi,ri1,rip1,xlami
      common/rwaves/egv,r1,rp1,r2,rp2,egvi,ri1,rip1,xlam(40),xlami(40)
      complex*16 cdiag,csubd,cc
      dimension diag(100),subd(100),cdiag(100),csubd(100)
      integer r,s
      maxdim=100
      msq=m**2
      m2=2*m
      nl=kmax1
      numax=nl-1
c find outside eigenvalues first
      c=cout
c find eigenvalues of even parity first
      s=0
c determine matrix size and check for overflow and other errors
      n=numax+3*idint(dsqrt(c))+2
      if((nl.le.0).or.(n.gt.maxdim)) goto 100
c set the diagonal elements.note that we write r=2*q+s
 5    r=s
      do 1 i=1,n
         xa=dfloat((m+r)*(m+r+1))
         diag(i)=(2.0d0*xa-dfloat(2*msq+1))*sgn*c*c/
     $        dfloat((m2+2*r-1)*(m2+2*r+3))+xa
 1    r=r+2
c now set the subdiagonal elements
      r=s+2
      do 2 i=2,n
         subd(i)=dsqrt(dfloat((m2+r)*(m2+r-1))*
     $        dfloat(r*(r-1))/dfloat((m2+2*r-3)*
     $        (m2+2*r+1)))*c*c/dfloat(m2+2*r-1)
 2    r=r+2
      call tql1(n,diag,subd,ia)
      if(ia.ne.0) goto 101
c tql1 returns eigenvalues of same parity in ascending order.
      j=0
      ls=s+1
      do 3 i=ls,nl,2
         xlam(i)=diag(j+1)
 3    j=j+1
      if(s.eq.1) goto 10
c find eigenvalues of odd parity
      s=1
      goto 5
cm*
 10   continue
cm* now do inside
      cc=cin
c find eigenvalues of even parity first
      s=0
c determine matrix size and check for overflow and other errors
      n=numax+3*idint(dsqrt(cdabs(cc)))+2
      if((nl.le.0).or.(n.gt.maxdim)) goto 100
c set the diagonal elements.note that we write r=2*q+s
 1005 r=s
      do 1001 i=1,n
         xa=dfloat((m+r)*(m+r+1))
         cdiag(i)=(2.0d0*xa-dfloat(2*msq+1))*sgn*cc*cc
     $        /dfloat((m2+2*r-1)*(m2+2*r+3))+xa
 1001 r=r+2
c now set the subdiagonal elements
      r=s+2
      do 1002 i=2,n
         csubd(i)=dsqrt(dfloat((m2+r)*(m2+r-1))*dfloat(r*(r-1))
     $        /dfloat((m2+2*r-3)*(m2+2*r+1)))*cc*cc/dfloat(m2+2*r-1)
 1002 r=r+2
      ib=0
      call sf02(n,cdiag,csubd,ib)
cm* this goes to setup for and execution of nag routine f02ajf and sorts
      if(ib.ne.0) goto 102
c sf02 returns eigenvalues of same parity in ascending order,
cm* with respect to the size of the real part (see onaka, asano's letter)
cm? is this completely general?
      j=0
      ls=s+1
      do 1030 i=ls,nl,2
         xlami(i)=cdiag(j+1)
 1030 j=j+1
      if(s.eq.1) return
c find eigenvalues of odd parity
      s=1
      goto 1005
cm*
cm* error returns
 100  if(nl.le.0) write(6,600)
      if (n.gt.maxdim) write(6,602) n
 600  format(20h-nl negative or zero)
 602  format(' -allowable array size exceeded. increase dimension 
     $   statements in sphwf to at least ',i4)
      return
 101  write(6,603) ia
 603  format(18h-tql1 error number,i4)
 102  write(6,604) ib
 604  format(20h-f02ajf error number,i4)
      return
      end
c
      subroutine tql1(n,d,e,ierr)
c
      integer i,j,l,m,n,ii,l1,mml,ierr
      real*8 d(n),e(n)
      real*8 b,c,f,g,h,p,r,s,machep
      real*8 dsqrt,dabs,dsign
c
c     this subroutine is a translation of the algol procedure tql1,
c     num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
c     wilkinson.
c     handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
c
c     this subroutine finds the eigenvalues of a symmetric
c     tridiagonal matrix by the ql method.
c
c     on input:
c
c        n is the order of the matrix;
c
c        d contains the diagonal elements of the input matrix;
c
c        e contains the subdiagonal elements of the input matrix
c          in its last n-1 positions.  e(1) is arbitrary.
c
c      on output:
c
c        d contains the eigenvalues in ascending order.  if an
c          error exit is made, the eigenvalues are correct and
c          ordered for indices 1,2,...ierr-1, but may not be
c          the smallest eigenvalues;
c
c        e has been destroyed;
c
c        ierr is set to
c          zero       for normal return,
c          j          if the j-th eigenvalue has not been
c                     determined after 30 iterations.
c
c     questions and comments should be directed to b. s. garbow,
c     applied mathematics division, argonne national laboratory
c
c     ------------------------------------------------------------------
c
c     :::::::::: machep is a machine dependent parameter specifying
c                the relative precision of floating point arithmetic.
c                machep = 16.0d0**(-13) for long form arithmetic
c                on s360 ::::::::::
c      data machep/1.d-14/
      real*8 epslon,epslo2
      common /accura/ epslon,epslo2
      machep=epslon
c
c
      ierr = 0
      if (n .eq. 1) go to 1001
c
      do 100 i = 2, n
  100 e(i-1) = e(i)
c
      f = 0.0d0
      b = 0.0d0
      e(n) = 0.0d0
c
      do 290 l = 1, n
         j = 0
         h = machep * (dabs(d(l)) + dabs(e(l)))
         if (b .lt. h) b = h
c     :::::::::: look for small sub-diagonal element ::::::::::
         do 110 m = l, n
            if (dabs(e(m)) .le. b) go to 120
c     :::::::::: e(n) is always zero, so there is no exit
c                through the bottom of the loop ::::::::::
  110    continue
c
  120    if (m .eq. l) go to 210
  130    if (j .eq. 30) go to 1000
         j = j + 1
c     :::::::::: form shift ::::::::::
         l1 = l + 1
         g = d(l)
         p = (d(l1) - g) / (2.0d0 * e(l))
         r = dsqrt(p*p+1.0d0)
         d(l) = e(l) / (p + dsign(r,p))
         h = g - d(l)
c
         do 140 i = l1, n
  140    d(i) = d(i) - h
c
         f = f + h
c     :::::::::: ql transformation ::::::::::
         p = d(m)
         c = 1.0d0
         s = 0.0d0
         mml = m - l
c     :::::::::: for i=m-1 step -1 until l do -- ::::::::::
         do 200 ii = 1, mml
            i = m - ii
            g = c * e(i)
            h = c * p
            if (dabs(p) .lt. dabs(e(i))) go to 150
            c = e(i) / p
            r = dsqrt(c*c+1.0d0)
            e(i+1) = s * p * r
            s = c / r
            c = 1.0d0 / r
            go to 160
  150       c = p / e(i)
            r = dsqrt(c*c+1.0d0)
            e(i+1) = s * e(i) * r
            s = 1.0d0 / r
            c = c * s
  160       p = c * d(i) - s * g
            d(i+1) = h + s * (c * g + s * d(i))
  200    continue
c
         e(l) = s * p
         d(l) = c * p
         if (dabs(e(l)) .gt. b) go to 130
  210    p = d(l) + f
c     :::::::::: order eigenvalues ::::::::::
         if (l .eq. 1) go to 250
c     :::::::::: for i=l step -1 until 2 do -- ::::::::::
         do 230 ii = 2, l
            i = l + 2 - ii
            if (p .ge. d(i-1)) go to 270
            d(i) = d(i-1)
  230    continue
c
  250    i = 1
  270    d(i) = p
  290 continue
c
      go to 1001
c     :::::::::: set error -- no convergence to an
c                eigenvalue after 30 iterations ::::::::::
 1000 ierr = l
 1001 return
c     :::::::::: last card of tql1 ::::::::::
      end
c
      subroutine bowkmp(n,m,egv,rn,begin)
      implicit real*8(a-h,o-z)
      common /accura/ epslon,epslo2
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
      dimension rn(100)
      integer begin
      beta(j,m)=dfloat(j*(j-1))*dfloat((2*m+j)*(2*m+j-1))/dfloat((2*m+2
     $     *j-3)*(2*m+2*j+1))/dfloat((2*m+2*j-1)**2)*c**4
      gamma(j,m)=dfloat((m+j)*(m+j+1))+(1.0d0-dfloat(4*m*m-1)/dfloat((2
     $     *m+2*j-1)*(2*m+2*j+3)))*.50d0*c**2*sgn
cm* this is the outside version
      c = cout
      icount=0
c kmax is read in as begin to set lower limit on d's,rn's in bowkmp
      kmax=begin
      k=3
      if(mod(n-m,2).eq.0) k=2
      idc=idint(c)/5
      if(mod(idc,2).ne.0) idc=idc+1
      nm2=n-m+2+idc
      nmk=nm2-k
      first1=0.0d0
      delta2=0.0d0
      nm5=16+idint((2*c+6)*dlog10(c))-idint(c)
      nm5=max0(nm5,kmax)
      if(mod(nm5,2).ne.0)nm5=nm5+1
      nm6=min0(102-k,nm5+nm2)
      do 7 i=nm6,100,2
         check=dsqrt(beta(i,m))/(gamma(i,m)-egv)
         if(dabs(check).le.5.0d-4*dabs(c)**2) goto 5
 7    continue
 5    begin=i
      call ndown(egv,nm2,begin,m,first1,rn)
      u2=-rn(nm2)
      start=gamma(k-2,m)-egv
      first=start
      rn(k) =first
      if(nm2.le.3) goto 37
      call nup(egv,nm2,k,m,first,rn)
 37   u1=rn(nm2)
      avg=(u1-u2)/2. 0d0
      rn(nm2)=avg
      denom1=1.0d0
      fe=1.0d0
      if(nmk.lt.2) goto 25
      do 22 j=2,nmk,2
         fe=fe*beta(nm2-j,m)/(rn(nm2-j  )*rn(nm2-j  ))
c     if(dabs(fe).le.1.0d-16) goto 25
         if(dabs(fe).le.epslon) goto 25
 22   denom1=denom1+fe
 25   denom2=avg*avg/beta(nm2,m)
      fe2=denom2
      nkm2=begin-nm2-2
      do 27 j=2,nkm2,2
         fe2=fe2*rn(nm2+j  )*rn(nm2+j  )/beta(nm2+j,m)
c     if(dabs(fe2).le.1.0d-16) goto 28
         if(dabs(fe2).le.epslon) goto 28
 27   denom2=fe2+denom2
 28   continue
      delta=(u1+u2)/(denom1+denom2)
      egv=egv+delta
      x=epslon*dmax1(dabs(egv),dabs(start))
      if(dabs(delta).le.x.and.dabs(delta2).le.x) goto 30
      icount=icount+1
      if(icount.lt.50) goto 501
 499  write(6,500) m,n
 500  format(1x, 51hbouwkamps method has not converged,50 iterations,m=
     $     ,i2,2x,2hn=,i2)
      goto 30
 501  continue
      i=min0(102-k,i+2)
      delta2=delta
      goto 5
 30   begin=begin-2
      return
      end
c
      subroutine sf02(n,cdiag,csubd,ifail)
      implicit real*8(a-h,o-z)
      logical wanted
      common/output/ wanted
      complex*16 cdiag(n),csubd(n)
      complex*16 dcmplx
      dimension ar(100,100),ai(100,100),rr(100),ri(100)
      dimension intger(100),irank(100)
      do 1 i=1,n
         do 1 j=1,n
            ar(i,j)=0.
            ai(i,j)=0.
 1    continue
cm* diagonal
      do 2 i=1,n
         ar(i,i)=dreal(cdiag(i))
         ai(i,i)=dimag(cdiag(i))
 2    continue
cm* subdiagonal, and superdiagonal since symmetric
      do 3 i=2,n
         arr=dreal(csubd(i))
         brr=dimag(csubd(i))
         j=i-1
         ar(i,j) = arr
         ar(j,i) = arr
         ai(i,j) = brr
         ai(j,i) = brr
 3    continue
      ifail = 1
      call f02ajf(ar,100,ai,100,n,rr,ri,intger,ifail)
      if(ifail.eq.0) go to 20
      write(6,*) ' f02ajf',ifail
      stop
cm* rank the real parts, then rearrange, using mark 12 nag routines
c
 20   continue
      ifail=0
      call m01daf(rr,1,n,'a',irank,ifail)
ccm* rank the real parts rr using a mark 1 nag routine
c20   ifail=1
c  call m01aaf(rr,1,n,irank,intger,ifail)
      if(ifail.eq.0) go to 21
      write(6,*) ' m01daf ',ifail
      stop
cm* sort into complex number
 21   do 9 i=1,n
         cdiag(irank(i)) =dcmplx(rr(i),ri(i))
 9    continue
c
      if(.not.wanted) return
      write(6,99994)
99994 format(1x,42hcomplex eigenvalues - even then odd parity)
      write(6,99995) (cdiag(i),i=1,n)
99995 format(1x,3(2e12.4,2x))
      return
      end
c
      subroutine nup(egv,lf,li,m,start,rn)
      implicit real*8(a-h,o-z)
      dimension rn(100)
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
      beta(j,m)=dfloat(j*(j-1))*dfloat((2*m+j)*(2*m+j-1))/dfloat((2*m+2
     $     *j-3)*(2*m+2*j+1))/dfloat((2*m+2*j-1)**2)*c**4
      gamma(j,m)=dfloat((m+j)*(m+j+1))+(1.0d0-dfloat(4*m*m-1)/dfloat((2
     $     *m+2*j-1)*(2*m+2*j+3)))*.50d0*c**2*sgn
      c = cout
      rn(li)=start
      if(lf.le.3)return
      lg=lf-2
      do 10 i=li,lg,2
 10   rn(i+2)=gamma(i,m)-egv-beta(i,m)/rn(i)
      return
      end
c
      subroutine ndown(egv,lf,li,m,start,rn)
      implicit real*8(a-h,o-z)
      dimension rn(100)
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
      beta(j,m)=dfloat(j*(j-1))*dfloat((2*m+j)*(2*m+j-1))/dfloat((2*m+2
     $     *j-3)*(2*m+2*j+1))/dfloat((2*m+2*j-1)**2)*c**4
      gamma(j,m)=dfloat((m+j)*(m+j+1))+(1.0d0-dfloat(4*m*m-1)/dfloat((2
     $     *m+2*j-1)*(2*m+2*j+3)))*.50d0*c**2*sgn
      c = cout
      rn(li)=start
      k=li-lf
      do 10 j=2,k,2
         jj=li-j
 10   rn(jj)=beta(jj,m)/(gamma(jj,m)-egv-rn(jj+2))
 30   return
      end
c
      subroutine cbowkm(n,m,egv,rn,begin)
      implicit real*8(a-h,o-z)
      common /accura/ epslon,epslo2
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
      integer begin
      complex*16 c,beta,gamma,egv,rn
      complex*16 egvsav,delsav
      dimension rn(100)
      complex*16 avg,czero,check,delta,delta2,denom1,denom2
      complex*16 fe,fe2,first,first1,start,u1,u2
      beta(j,m)=dfloat(j*(j-1))*dfloat((2*m+j)*(2*m+j-1))/dfloat((2*m+2
     $     *j-3)*(2*m+2*j+1))/dfloat((2*m+2*j-1)**2)*c**4
      gamma(j,m)=dfloat((m+j)*(m+j+1))+(1.0d0-dfloat(4*m*m-1)/dfloat((2
     $     *m+2*j-1)*(2*m+2*j+3)))*.50d0*c**2*sgn
      egvsav=egv
      c = cin
      icount=0
c kmax is read in as begin to set lower limit on d's,rn's in bowkmp
      kmax=begin
      k=3
      if(mod(n-m,2).eq.0) k=2
cm*
      ac=cdabs(c)
      idc=idint(ac)/5
      if(mod(idc,2).ne.0) idc=idc+1
      nm2=n-m+2+idc
      nmk=nm2-k
      czero= (0.,0.)
      first1=czero
      delta2=czero
      nm5=16+idint((2*ac+6)*dlog10(ac))-idint(ac)
      nm5=max0(nm5,kmax)
      if(mod(nm5,2).ne.0)nm5=nm5+1
      nm6=min0(102-k,nm5+nm2)
      do 7 i=nm6,100,2
         check=cdsqrt(beta(i,m))/(gamma(i,m)-egv)
cm? probably ok since just an integer estimate
         if(cdabs(check).le.5.0d-4*ac**2) goto 5
 7    continue
 5    begin=i
      call cndown(egv,nm2,begin,m,first1,rn)
      u2=-rn(nm2)
      start=gamma(k-2,m)-egv
      first=start
      rn(k) =first
      if(nm2.le.3) goto 37
      call cnup(egv,nm2,k,m,first,rn)
 37   u1=rn(nm2)
      avg=(u1-u2)/2. 0d0
      rn(nm2)=avg
      denom1=(1.0d0,0.)
      fe=(1.0d0,0.)
      if(nmk.lt.2) goto 25
      do 22 j=2,nmk,2
         fe=fe*beta(nm2-j,m)/(rn(nm2-j  )*rn(nm2-j  ))
cm* underflow on modulus
         if(cdabs(fe).le.epslon) goto 25
 22   denom1=denom1+fe
 25   denom2=avg*avg/beta(nm2,m)
      fe2=denom2
      nkm2=begin-nm2-2
      do 27 j=2,nkm2,2
         fe2=fe2*rn(nm2+j  )*rn(nm2+j  )/beta(nm2+j,m)
cm* underflow on modulus
         if(cdabs(fe2).le.epslon) goto 28
 27   denom2=fe2+denom2
 28   continue
      delta=(u1+u2)/(denom1+denom2)
      egv=egv+delta
      if(icount.eq.0) delsav=delta
cm*  x=1.e-14*dmax1(dabs(egv),dabs(start))
cm*  if(dabs(delta).le.x.and.dabs(delta2).le.x) goto 30
cm* check convergence on real and imaginary parts separately
      ereal=dreal(egv)
      eimag=dimag(egv)
      sreal=dreal(start)
      simag=dimag(start)
c slight relaxation
      xmagni=3.d0
      xreal=xmagni*epslon*dmax1(dabs(ereal),dabs(sreal))
      ximag=xmagni*epslon*dmax1(dabs(eimag),dabs(simag))
      dereal=dreal(delta)
      deimag=dimag(delta)
      derea2=dreal(delta2)
      deima2=dimag(delta2)
      der=dmax1(dabs(dereal),dabs(derea2))
      dei=dmax1(dabs(deimag),dabs(deima2))
      if(der.le.xreal .and. dei.le.ximag) goto 30
      icount=icount+1
      if(icount.lt.50) goto 501
 499  write(6,500) m,n
 500  format(1x,'cbowkm: 50 iterations, m=',i2,2x,'n=',i2)
      write(6,505) delta,delta2,xreal,ximag
 505  format(1x,3(1pd11.3,1x,1pd11.3),2x)
      goto 30
 501  continue
      i=min0(102-k,i+2)
      delta2=delta
      goto 5
 30   continue
c test if eigenvalue has jumped to a new root
      if(cdabs(egv-egvsav).gt. 2.d0*cdabs(delsav)) 
     $     write(6,506) egvsav,egv
 506  format(1x,'cbowkm warning: possible eigenvalue switch',
     $     2(1pd11.3,1x,1pd11.3),2x)
      begin=begin-2
      return
      end
c
      subroutine cnup(egv,lf,li,m,start,rn)
      implicit real*8(a-h,o-z)
cm* complex inside version
      complex*16 egv,start,rn(100)
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
      complex*16 c,beta,gamma
      beta(j,m)=dfloat(j*(j-1))*dfloat((2*m+j)*(2*m+j-1))/dfloat((2*m+2
     $     *j-3)*(2*m+2*j+1))/dfloat((2*m+2*j-1)**2)*c**4
      gamma(j,m)=dfloat((m+j)*(m+j+1))+(1.0d0-dfloat(4*m*m-1)/dfloat((2
     $     *m+2*j-1)*(2*m+2*j+3)))*.50d0*c**2*sgn
      c = cin
      rn(li)=start
      if(lf.le.3)return
      lg=lf-2
      do 10 i=li,lg,2
 10   rn(i+2)=gamma(i,m)-egv-beta(i,m)/rn(i)
      return
      end
c
      subroutine cndown(egv,lf,li,m,start,rn)
      implicit real*8(a-h,o-z)
cm* complex inside version
      complex*16 egv,start,rn(100)
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
      complex*16 c,beta,gamma
      beta(j,m)=dfloat(j*(j-1))*dfloat((2*m+j)*(2*m+j-1))/dfloat((2*m+2
     $     *j-3)*(2*m+2*j+1))/dfloat((2*m+2*j-1)**2)*c**4
      gamma(j,m)=dfloat((m+j)*(m+j+1))+(1.0d0-dfloat(4*m*m-1)/dfloat((2
     $     *m+2*j-1)*(2*m+2*j+3)))*.50d0*c**2*sgn
      c = cin
      rn(li)=start
      k=li-lf
      do 10 j=2,k,2
         jj=li-j
 10   rn(jj)=beta(jj,m)/(gamma(jj,m)-egv-rn(jj+2))
 30   return
      end
c
      function fact(n)
      implicit real*8(a-h,o-z)
      fact=1.0d0
      if(n.le.1)return
      nn=n-1
      fact=dfloat(n)
      do 10 i=1,nn
 10   fact=fact*dfloat(n-i)
      return
      end
c
      subroutine excoef(n,m,rn,jj,jdn,jup,d,ddlog,dsgn)
cm* calculates the expansion coefficients
c calculates d(r) and stores as d(r+2)
      implicit real*8(a-h,o-z)
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
      dimension rn(100),d(100)
      dimension fsum(100)
      logical danger,okay
      ff(k,r) =-dfloat((2*m+2*k-1)*(2*m+2*k+1))*r /
     $     (dfloat((2*m+k)*(2*m+k-1))*c**2)*sgn
cm* note that ssum is not used
cm*      ssum(j)=dfloat((-1)**((j-l)/2))*fact(j+2*m+l)/
cm*     $     (2.0d0**j*fact((j-l)/2)*fact((j+2*m+l)/2))
      c = cout
      jj=jj-2
      k=3
      if(mod(n-m,2).eq.0)k=2
      l=k-2
      if(c.ne.0)goto 5
      do 10 j=k,jj,2
 10   d(j)=0.0d0
      d(n-m+2)=1.0d0
      jup=jj
      jdn=k
      return
 5    f1=10.0d0**(m-n)
      f=f1
      d(k)=f
      do 20 j=k,jj,2
         if(dabs(f).le.1.0d-70)goto 21
         f=f*ff(j,rn(j))
 20   d(j+2)=f
 21   jj=j-2
      sum=0.0d0
      jj2=jj+2
      if(l.eq.0)fsum(k)=fact(2*m)/fact(m)
      if(l.eq.1)fsum(k)=fact(2*m+2)/2.0d0/fact(m+1)
      m2l=2*m+l
      do 25 jj=k,jj2,2
         sum=fsum(jj)*d(jj)+sum
         fsum(jj+2)=-fsum(jj)*dfloat((jj+m2l)*(jj-1+m2l))/
     $        dfloat((jj-l)*(jj+m2l))
 25   continue
      total=fsum(n-m+2)
      danger=.false.
      ddlog=dlog  (dabs(total))-dlog  (dabs(sum))+dlog(f1)
      ddlog1=ddlog-dlog(f1)
      dsgn=dsign(1.0d0,total)*dsign(1.0d0,sum)
      if(ddlog1.le.-177.0d0)danger=.true.
      if(.not.danger)change=total/sum
      jdn=k
      jup=min0(102-k,jj2)
      okay=.false.
      do 30 j=jdn,jup,2
         if(okay)goto 31
         if(ddlog1+dlog  (dabs(d(j))).ge.-170.0d0)goto 31
         d(j)=0.0d0
         goto 30
 31      if(danger)goto 32
         d(j)=change*d(j)
         if(dabs(d(j)).gt.1.0d-60)okay=.true.
         if(dabs(d(j)).le.1.0d-60)okay=.false.
         goto 30
 32      temp=d(j)*total
         d(j)=temp/sum
         if(dabs(d(j)).gt.1.0d-60)okay=.true.
         if(dabs(d(j)).le.1.0d-60)okay=.false.
 30   continue
      return
      end
c
      subroutine cexcoe(n,m,rn,jj,jdn,jup,d,ddlog,dsgn)
cm* calculates the expansion coefficients on inside, where ddlog and
cm* dsgn are used for normalization, and are not used except in this routine
c   calculates d(r) and stores as d(r+2)
      implicit real*8(a-h,o-z)
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
      complex*16 rn(100),d(100)
      dimension fsum(100)
      complex*16 c,change,czero,f,ff,r,sum,temp
      logical danger,okay
      ff(k,r) =-dfloat((2*m+2*k-1)*(2*m+2*k+1))*r /
     $     (dfloat((2*m+k)*(2*m+k-1))*c**2)*sgn
cm* note that ssum is not used
cm*      ssum(j)=dfloat((-1)**((j-l)/2))*fact(j+2*m+l)/
cm*     $     (2.0d0**j*fact((j-l)/2)*fact((j+2*m+l)/2))
      c = cin
      czero=(0.,0.)
      jj=jj-2
      k=3
      if(mod(n-m,2).eq.0)k=2
      l=k-2
      if(c.ne.czero)goto 5
      do 10 j=k,jj,2
 10   d(j)=czero
      d(n-m+2)=(1.0d0,0.)
      jup=jj
      jdn=k
      return
 5    f1=10.0d0**(m-n)
      f=f1
      d(k)=f
      do 20 j=k,jj,2
cm* underflow on modulus
         if(cdabs(f).le.1.0d-70)goto 21
         f=f*ff(j,rn(j))
 20   d(j+2)=f
 21   jj=j-2
      sum=czero
      jj2=jj+2
      if(l.eq.0)fsum(k)=fact(2*m)/fact(m)
      if(l.eq.1)fsum(k)=fact(2*m+2)/2.0d0/fact(m+1)
      m2l=2*m+l
      do 25 jj=k,jj2,2
         sum=fsum(jj)*d(jj)+sum
         fsum(jj+2)=-fsum(jj)*dfloat((jj+m2l)*(jj-1+m2l))/
     $        dfloat((jj-l)*(jj+m2l))
 25   continue
      total=fsum(n-m+2)
      danger=.false.
cm* renormalize
      ddlog1=dlog(dabs(total)) - dlog(cdabs(sum))
      ddlog=ddlog1+dlog(f1)
cm* dsgn=dsign(1.0d0,total)*dsign(1.0d0,sum)
cm* dsgn(i) not used here or elsewhere - set here to cause problems if used
      dsgn=1.e+99
      if(ddlog1.le.-177.0d0)danger=.true.
      if(.not.danger)change=total/sum
      jdn=k
      jup=min0(102-k,jj2)
      okay=.false.
      do 30 j=jdn,jup,2
         if(okay)goto 31
         if(ddlog1+dlog(cdabs(d(j))).ge.-170.0d0)goto 31
         d(j)=czero
         goto 30
 31      if(danger)goto 32
         d(j)=change*d(j)
         if(cdabs(d(j)).gt.1.0d-60)okay=.true.
         if(cdabs(d(j)).le.1.0d-60)okay=.false.
         goto 30
 32      temp=d(j)*total
         d(j)=temp/sum
         if(cdabs(d(j)).gt.1.0d-60)okay=.true.
         if(cdabs(d(j)).le.1.0d-60)okay=.false.
 30   continue
      return
      end
c
      subroutine dsubro(n,m,egv,d2m,d2m1,drho,ldrho)
c calculates d(rho!r) stores as drho(r)
c if m=0 d2m=d(2),d2m1=d(3) otherwise dneg(2m),dneg(2m-1)
      implicit real*8(a-h,o-z)
      common /accura/ epslon,epslo2
      dimension drho(100),aprx(100),q(250)
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
      a(j)=dfloat((2*m+j)*(2*m+j-1))/dfloat((2*m+2*j-1)*(2*m+2*j+1))*c*
     $     c*sgn
      b(j)=dfloat((m+j)*(m+j+1))-egv+dfloat(2*(m+j)*(m+j+1)-2*m*m-1)/df
     $     loat((2*m+2*j-1)*(2*m+2*j+3))*c*c*sgn
      cc(j)=dfloat((j+1)*(j+2))/dfloat((2*m+2*j+1)*(2*m+2*j+3))*c*c*sgn
cm* only used on outside
      c = cout
      lim=20+2*(idint((2*c+6)*dlog10(c))/2)
      lim=min0(lim,100)
      if(mod(n-m,2).eq.0) goto 5
      ink=2*m+1
      f=d2m1
      aa=-c*c/dfloat((1-2*m)*(3-2*m))*sgn
      lim=lim-1
      goto 10
 5    ink=2*m+2
      f=d2m
      aa=-c*c/dfloat((2*m+1)*(1-2*m))*sgn
 10   nuj=2*lim+2+ink
      do 15 j=ink,lim,2
 15   drho(j)=0.0d0
      if(f.eq.0.0d0) ldrho=lim
      if(f.eq.0.0d0) return
      nk=ink
 18   nuj2=nuj-2
      q(nuj)=b(-nuj)
      do 20 kk= nk,nuj2,2
         k=nuj2+ nk-kk
 20   q(k)=b(-k)-cc(-k-2)*a(-k)/q(k+2)
      aprx(ink)=aa/q(ink)
      lim2=lim-2
      do 45 j=ink,lim2,2
 45   aprx(j+2)=-a(-j)/q(j+2)
      do 50 j=ink,lim,2
         if(dabs(aprx(j)-drho(j)).le.epslon*dabs(aprx(j))) goto 50
         nk=j
         do 55 k=ink,lim,2
 55      drho(k)=aprx(k)
         nuj=nuj+10
         goto 18
 50   continue
      do 60 j=ink,lim,2
         if(dlog(dabs(drho(j)))+dlog(dabs(f)).le.-170) goto 65
         drho(j)=drho(j)*f
 60   f=drho(j)
 65   ldrho=j-2
      if(ldrho.ge.ink) return
      drho(ink)=0.0d0
      ldrho=ink
      return
      end
c
      subroutine xpaneg(n,m,egv,d2,d3,dneg,ddlog,dsgn,ddlog2,dsgn2)
c calculates d(-r) and stores as dneg(r)
      implicit real*8(a-h,o-z)
      dimension dneg(40),q(42)
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
      a(j)=dfloat((2*m+j)*(2*m+j-1))/dfloat((2*m+2*j-1)*(2*m+2*j+1))*c*
     $     c*sgn
      b(j)=dfloat((m+j)*(m+j+1))-egv+dfloat(2*(m+j)*(m+j+1)-2*m*m-1)/df
     $     loat((2*m+2*j-1)*(2*m+2*j+3))*c*c*sgn
      cc(j)=dfloat((j+1)*(j+2))/dfloat((2*m+2*j+1)*(2*m+2*j+3))*c*c*sgn
cm* only used on outside
      c = cout
      if(m.eq.0) return
      if(mod(n-m,2).eq.0) goto 5
      ink=1
      lim=2*m-1
      f=d3
      goto 10
 5    ink=2
      lim=2*m
      f=d2
 10   lim2=lim-2
      q(lim)=b(-lim)
      if(lim2.lt.ink) goto 25
      do 20 jj=ink,lim2,2
         j=lim2+ink-jj
 20   q(j)=b(-j)-cc(-j-2)*a(-j)/q(j+2)
 25   if(f.eq.0.0d0) goto 40
      do 30 j=ink,lim,2
         t=-a(-j+2)/q(j)
         if(t.eq.0d0) goto 40
         if(dlog(dabs(t))+dlog(dabs(f)).le.-175.0d0)goto 40
         dneg(j)=t*f
 30   f=dneg(j)
      return
 40   dsgn2=dsgn
      ddlog2=ddlog
      do 45 j=ink,lim,2
         term=-a(-j+2)/q(j)
         dsgn2=dsign(1.0d0,term)*dsgn2
         ddlog2=dlog(dabs(term))+ddlog2
         if(f.eq.0.0d0.or.ddlog2.le.-177.0d0) dneg(j)=0.0d0
 45   continue
      return
      end
c
      subroutine radial(n,m,jdn,jup,d,kbs,max,r,rp)
cm*this is for the outside only
c  calculates asymptotic form of radial spheroidal wavefunction
c  and its derivative
      implicit real*8(a-h,o-z)
      real*8 kbs(max)
      dimension d(jup),term(100)
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
      c = cout
      sum1=0.0d0
      sum2=0.0d0
      sum3=0.0d0
      m2=2*m
c scale factor for fjm,sum1,sum2,sum3,prevents overflows
      scale=10.0d0**(-m2)
      fjm=fact(m2+jdn-2)*scale
      do 10 jj=jdn,jup,2
         j=jj-2
         term(jj)=0.0d0
         if(d(jj).eq.0.0d0) goto 5
         if(dlog(dabs(d(jj)))+dlog(fjm).le.-170.0d0) goto 5
         term(jj)=fjm*d(jj)
 5       fjm=fjm*dfloat((m2+jj)*(m2+jj-1))/dfloat(jj*(jj-1))
         sum1=sum1+term(jj)
 10   continue
 12   do 14 jj=jdn,jup,2
         j=jj-2
         if(term(jj).eq.0.0d0.or.kbs(m+j+1).eq.0.0d0) goto 14
         i=(j+m-n)/2
         if(dlog(dabs(term(jj)))+dlog(dabs(kbs(m+j+1))).le.-177.0d0)
     $        goto14
         term2=term(jj)*dfloat((-1)**i)*kbs(m+j+1)
         sum2=sum2+term2
 14   continue
 16   do 18 jj=jdn,jup,2
         if(term(jj).eq.0.0d0) goto 18
         j=jj-2
         i=(j+m-n)/2
         if(j+m.ne.0)goto 17
         if(dlog(dabs(term(jj)))+dlog(dabs(kbs(2)*c)).le.-177.0d0)
     $        goto18
         sum3=-c*kbs(2)*term(jj)*dfloat((-1)**i)
         goto 18
 17      f1=c*dfloat(j+m)/dfloat(2*j+2*m+1)*kbs(j+m)
         f2=-c*dfloat(j+m+1)/dfloat(2*j+2*m+1)*kbs(j+m+2)
         f3=sgn*dfloat(m)/(zeta**3-sgn*zeta)*kbs(j+m+1)
         f=f1+f2+f3
         if(f.eq.0.0d0) goto 18
         if(dlog(dabs(f))+dlog(dabs(term(jj))).le.-177.0d0)goto 18
         term3=dfloat((-1)**i)*term(jj)*f
         sum3=sum3+term3
 18   continue
c loss of sig figs in prolate case if zeta near 1 could use zeta-1
 20   rnorm=((zeta*zeta-sgn)/zeta/zeta)**(dfloat(m)/2.0d0)/sum1
      r=rnorm*sum2
      rp=rnorm*sum3
      return
      end
c
      subroutine cradia(n,m,jdn,jup,d,kbs,max,r,rp)
cm* for the inside - complex form of radial
c calculates asymptotic form of radial spheroidal wavefunction
c and its derivative
      implicit real*8(a-h,o-z)
cm* new complex typing
      complex*16 d(jup),kbs(max),r,rp
      complex*16 term(100)
      complex*16 c,czero,f,f1,f2,f3,rnorm,sum1,sum2,sum3,term2,term3
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
      c = cin
      czero=(0.,0.)
      sum1=czero
      sum2=czero
      sum3=czero
      m2=2*m
c scale factor for fjm,sum1,sum2,sum3,prevents overflows
      scale=10.0d0**(-m2)
      fjm=fact(m2+jdn-2)*scale
      do 10 jj=jdn,jup,2
         j=jj-2
         term(jj)=czero
         if(d(jj).eq.czero) goto 5
cm* underflow on modulus
         if(dlog(cdabs(d(jj)))+dlog(fjm).le.-170.0d0) goto 5
         term(jj)=fjm*d(jj)
 5       fjm=fjm*dfloat((m2+jj)*(m2+jj-1))/dfloat(jj*(jj-1))
         sum1=sum1+term(jj)
 10   continue
 12   do 14 jj=jdn,jup,2
         j=jj-2
         if(term(jj).eq.czero.or.kbs(m+j+1).eq.czero) goto 14
         i=(j+m-n)/2
cm* underflow on modulus
         if(dlog(cdabs(term(jj)))+dlog(cdabs(kbs(m+j+1)))
     $        .le.-177.d0)goto14
         term2=term(jj)*dfloat((-1)**i)*kbs(m+j+1)
         sum2=sum2+term2
 14   continue
 16   do 18 jj=jdn,jup,2
         if(term(jj).eq.czero) goto 18
         j=jj-2
         i=(j+m-n)/2
         if(j+m.ne.0)goto 17
         if(dlog(cdabs(term(jj)))+dlog(cdabs(kbs(2)*c)).le.-177.d0)
     $        goto18
         sum3=-c*kbs(2)*term(jj)*dfloat((-1)**i)
         goto 18
 17      f1=c*dfloat(j+m)/dfloat(2*j+2*m+1)*kbs(j+m)
         f2=-c*dfloat(j+m+1)/dfloat(2*j+2*m+1)*kbs(j+m+2)
         f3=sgn*dfloat(m)/(zeta**3-sgn*zeta)*kbs(j+m+1)
         f=f1+f2+f3
         if(f.eq.czero) goto 18
         if(dlog(cdabs(f))+dlog(cdabs(term(jj))).le.-177.d0)
     $        goto 18
         term3=dfloat((-1)**i)*term(jj)*f
         sum3=sum3+term3
 18   continue
c loss of sig figs in prolate case if zeta near 1 could use zeta-1
 20   rnorm=((zeta*zeta-sgn)/zeta/zeta)**(dfloat(m)/2.0d0)/sum1
      r=rnorm*sum2
      rp=rnorm*sum3
      return
      end
c
      subroutine obzero(n,m,dneg,d,jup,r1,rp1,r2,rp2)
c      calculates oblate radial wavefunctions and derivatives
c in the special case zeta=0.0
      implicit real*8(a-h,o-z)
      dimension dneg(40),d(jup)
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
cm* only used on outside
      c = cout
      pi=3.141592653589793d0
      sum=0.0d0
      r1=0.0d0
      rp1=0.0d0
      nm=n-m
      m2=2*m
      if(mod(nm,2).eq.0) goto 15
      fj=fact(m2+1)
      if(m.eq.0) d2m1=d(3)
      if(m.ne.0) d2m1=dneg(2*m-1)
      rp1a=dfloat((-1)**((nm-1)/2))*2.0d0**m*fact(m)*c**(m+1)*d(3)
      do 10 jj=3,jup,2
         j=jj-2
         sum=sum+fj*d(jj)
 10   fj=fj*dfloat((m2+jj)*(m2+jj-1))/dfloat(jj*(jj-1))
      rp1=rp1a/sum/dfloat( 2*m+3)
      r2=-1.0d0/c/rp1
      rp2a=dfloat((-1)**((nm-1)/2)*(2*m-3)*(2*m-1))*fact(m)*c**(m-2)*pi
      rp2b=2.0d0**(2*n-m+1)*fact(2*m)*d2m1*sum
      rp2=rp2a/rp2b*(fact(n+m+1)/fact((nm-1)/2)/fact((n+m+1)/2))**2
      return
 15   r1a=dfloat((-1)**(nm/2))*2.0d0**m*fact(m)*c**m*d(2)
      fj=fact(m2)
      if(m.eq.0) d2m=d(2)
      if(m.ne.0) d2m=dneg(2*m)
      do 20 jj=2,jup,2
         j=jj-2
         sum=sum+fj*d(jj)
 20   fj=fj*dfloat((m2+jj)*(m2+jj-1))/dfloat(jj*(jj-1))
      r1=r1a/dfloat(2*m+1)/sum
      r2a=dfloat((-1)**(nm/2)*(2*m-1))*fact(m)*c**(m-1)*pi
      r2b=2.0d0**(2*n-m+1)*fact(2*m)*d2m*sum
      r2=r2a/r2b*(fact(n+m)/fact(nm/2)/fact((n+m)/2))**2
      rp2=1.0d0/c/r1
      return
      end
c
      subroutine baber(n,m,egv,an,ap,lima)
c calculates coefficients for the baber+hasse expansion of oblate
c radial wavefunctions of nonlarge c*zeta
c stores as a(-j)=an(j) a(j)=ap(j+1)
      implicit real*8 (a-h,o-z)
      dimension an(20),ap(100),q(20),p(100)
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
      aa(k)=2.0d0*c*dfloat((k+m)*(k+2*m))/dfloat(2*k+2*m+1)
      bb(k)=-(dfloat((k+m)*(k+m+1))-egv-c*c)
      cc(k)=-2.0d0*c*dfloat((k+1)*(k+m+1))/dfloat(2*k+2*m+1)
cm* only used on outside
      c = cout
      idc=idint(c)/5
      if(mod(idc,2).ne.0) idc=idc+1
      nm=n-m+idc
      nm1=nm+1
      lima=16+6*idint((.1*c+1)*dlog10(c))+nm
      lima=min0(98,lima)
      ap(lima+2)=-c/lima*10.0d0**(n-m)
      ap(lima+1)=10.0d0**(n-m)
      do 10 jj=nm1 ,lima
         j=lima+nm1-jj
 10   ap(j)=-((bb(j)+aa(j+1)*ap(j+2)/ap(j+1))/cc(j-1))*ap(j+1)
      if(nm.eq.0) goto 17
      p(1)=bb(0)
      do 12 j=1,nm
 12   p(j+1)=bb(j)-cc(j-1)*aa(j)/p(j)
      do 15 jj=1,nm
         j=nm+1-jj
 15   ap(j)=-aa(j)/p(j)*ap(j+1)
 17   continue
      lima=lima+2
      if(m.eq.0) return
      q(m)=bb(-m)
      if(m.eq.1) goto 25
      m1=m-1
      do 20 jj=1,m1
         j=m1-jj+1
 20   q(j)=bb(-j)-(cc(-j-1)*aa(-j))/q(j+1)
 25   an1=ap(1)
      do 30 j=1,m
         an(j)=-aa(-j+1)/q(j)*an1
 30   an1=an(j)
      return
      end
c
      subroutine radob(n,m,r1,rp1,qh,nmax,mmax,an,ap,lima,r2,rp2)
c calculates the baber+hasse expansions for the oblate radial
c wavefunction and its derivative for nonlarge c*zeta
cm* choice made vs result of radial
      implicit real*8(a-h,o-z)
      dimension qh(nmax,mmax),an(20),ap(lima)
      complex*16 i,r3,rp3,sum1,sum2,fr,term,term1
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
cm* only used on outside
      c = cout
      pi=3.141592653589793d0
      i=(0.0d0,1.0d0)
      fr=cdexp(i*(c*zeta-dfloat(n+1)*pi/2.0d0))*(-i)**m/fact(m)/c
c calculate the wavefunction
      sum1=(0.0d0,0.0d0)
      if(m.eq.0) anm=ap(1)
      if(m.eq.0) goto 15
      anm=an(m)
      do 10 j=1,m
 10   sum1=sum1+an(j)/anm*qh(m-j+1,m+1)/(-i)**j
 15   do 20 jj=1,lima
         j=jj-1
         term1=(-i)**j*ap(jj)*qh(m+j+1,m+1)
         term=term1/anm
 20   sum1=sum1+ term
      r3=fr*sum1
      r2=-i*r3
c calculate the derivative
      sum2=(0.0d0,0.0d0)
      zeta2=zeta*zeta+1.0d0
      if(m.eq.0) goto 35
      do 30 j=1,m
         sum2=sum2+an(j)/anm*(qh(m-j+1,m+1)*(i*c-dfloat(m-j+1)
     $        *zeta/zeta2)-dfloat(1-j)*qh(m-j+2,m+1)/zeta2)/(-i)**j
 30   continue
 35   do 40 jj=1,lima
         j=jj-1
         term1=(-i)**j*ap(jj)*(qh(m+j+1,m+1)*(i*c-dfloat(m+j+1)*
     $        zeta/zeta2)-dfloat(1+j)*qh(m+j+2,m+1)/zeta2)
         term=term1/anm
 40   sum2=sum2+term
      rp3=fr*sum2
      rp2=-i*rp3
      return
      end
c
      subroutine radpr(n,m,dneg,drho,ldrho,d,jup,qneg,q,p,maxdeg,
     $     maxord,lsig,y,r2,rp2,ddlog,dsgn,ddlog2,dsgn2)
c prolate radial wavefunctions for nonlarge c*zeta
cm* choice made vs result of radial
      implicit real*8(a-h,o-z)
      dimension dneg(40),drho(ldrho),d(jup),qneg(maxord)
      dimension q(maxdeg,maxord),p(maxdeg,maxord)
      real*8 kappa
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
cm* only used on outside
      c = cout
c for zeta near one
      if(lsig.eq.0) zeta1=zeta*zeta-1.0d0
      if(lsig.eq.1) zeta1=y*(y+2.0d0)
c joining factor kappa
      sum=0.0d0
      l=0
      ll=0
      if(mod(n-m,2).ne.0) l=1
      if(mod(m,2).ne.0) ll=1
      l2=l+2
      m2=2*m
c scale factor for fjm,prevents overflows
      scale=10.0d0**(-m2)
      fjm=fact(m2+l)*scale
      do 10 jj=l2,jup,2
         j=jj-2
         sum=sum+fjm*d(jj)
         fjm=fjm*dfloat((m2+jj)*(m2+jj-1))/dfloat(jj*(jj-1))
 10   continue
      sum=sum/scale
c dfact=fact((n-m-l)/2)*fact((n+m+l)/2)/fact(n+m+l)
      dfact=1.0d0
      lim1=n+m+l
      lim2=lim1/2
      lim3=(n-m-l)/2
      do 5 j=1,lim1
         dj1=dfloat(j)
         dj2=1.0d0
         dj3=1.0d0
         if(j.le.lim2) dj2=dj1
         if(j.le.lim3) dj3=dj1
 5    dfact=dj2*dj3/dj1*dfact
      top=2.0d0**(n-m)*fact(2*m)*dfact
      bot=dfloat(2*m-1)*fact(m)*c**(m-1-l)
      f1=d(l+2)
      if(m.ne.0) f1=dneg(2*m-l)
      f2=1.0d0
      if(l.ne.0) f2=dfloat(2*m-3)
      kappa=dfloat((-1)**l)*top/bot*f1/f2*sum
      if(f1.ne.0.0d0) goto 12
      term=top*sum/f2/bot
      if(m.eq.0) goto 11
      caplog=dlog(dabs(term))+ddlog2
      capsgn=dsign(1.0d0,dfloat((-1)**l)*term)*dsgn2
      goto 12
 11   caplog=dlog  (dabs(term))+ddlog
      capsgn=dsign(1.0d0,dfloat((-1)**l)*term)*dsgn
 12   continue
c calculate wavefunction r2
      sum1=0.0d0
      sum2=0.0d0
      sum3=0.0d0
      sum4=0.0d0
      if(m.eq.0) goto 40
      if(m.ne.1) goto 15
      qq=qneg(1)
      if(l.ne.0) qq=q(1,2)
      sum1=dneg(2*m-l)*qq
      goto 40
 15   continue
      if(l.eq.1.and.ll.eq.1) m2l=m+2
      if(l.ne.1.or.ll.ne.1) m2l=m+2-ll-l
      m2=2*m-l
      do 20 j=m2l,m2,2
 20   sum1=sum1+dneg(j)*qneg(j-m)
      ml=m2l-2
      ls=2-l
      do 30 j=ls,ml,2
 30   sum2=sum2+dneg(j)*q(m-j+1,m+1)
 40   do 50 jj=l2,jup,2
         if(d(jj).eq.0.0d0) goto 50
         j=jj-2
         if(dabs(q(m+j+1,m+1)).le.1.0d-76/dabs(d(jj)))goto 50
         term=d(jj)*q(m+j+1,m+1)
         sum3=sum3+term
 50   continue
 55   mp2=2*m+2-l
      do 60 j=mp2,ldrho,2
 60   sum4=sum4+drho(j)*p(j-m,m+1)
      tsum=sum1+sum2+sum3+sum4
      if(f1.ne.0.0d0)goto 62
      rlog=dlog(dabs(tsum))
      rsgn=dsign(1.0d0,tsum)
      r2=rsgn*capsgn*dexp(rlog-caplog)
      goto 65
 62   r2=tsum/kappa
 65   continue
c calculate the derivative of the wavefunction  rp2
      sum1=0.0d0
      sum2=0.0d0
      sum3=0.0d0
      if(m.ne.0) goto 80
      do 70 jj=l2,jup,2
         j=jj-2
 70   sum1=sum1+d(jj)*q(j+1,2)
      sum1=sum1/dsqrt(zeta1)
      goto 140
 80   if(m.ne.1) goto 85
      if(l.eq.0)sum1=(zeta*qneg(1)+dsqrt(zeta1)*(zeta*q(1,3)+q(2,3))/2.
     $     0d0)*dneg(2)/zeta1
      if(l.eq.1)sum1=-qneg(1)/zeta1*dneg(1)
      goto 120
 85   do 90 j=m2l,m2,2
         if(j-m.eq.1) qq=q(1,m+1)
         if(j-m.ne.1) qq=qneg(j-m-1)
         sum1=sum1+dneg(j)*(dfloat(1-j)*qq-dfloat(m-j+1)*
     $        zeta*qneg(j-m))/zeta1
 90   continue
      do 100 j=ls,ml,2
         sum1=sum1+dneg(j)*(-dfloat(m)*zeta*q(m-j+1,m+1)/
     $        zeta1+dfloat((2*m-j)*(1-j))*q(m-j+1,m)/dsqrt(zeta1))
 100  continue
 120  do 130 jj=l2,jup,2
         j=jj-2
         sum2=sum2+d(jj)*(-dfloat(m)*zeta*q(m+j+1,m+1)/
     $        zeta1+dfloat((2*m+j)*(j+1))*q(m+j+1,m)/dsqrt(zeta1))
 130  continue
 140  do 150 j=mp2,ldrho,2
         sum3=sum3+drho(j)*(p(j-m,m+2)/dsqrt(zeta1)+
     $        dfloat(m)*zeta*p(j-m,m+1)/zeta1)
 150  continue
      tsum=sum1+sum2+sum3
      if(f1.ne.0.0d0) goto 152
      rlog=dlog(dabs(tsum))
      rsgn=dsign(1.0d0,tsum)
      rp2=rsgn*capsgn*dexp(rlog-caplog)
      return
 152  rp2=tsum/kappa
      return
      end
c 
      subroutine capfun
      implicit real*8 (a-h,o-z)
      complex*16 uo3,yo3,uo1,yo1,ui1,yi1,vo3,xo3,vo1,xo1,vi1,xi1
      real*8 ii
      complex*16 czero,ar
      complex*16 dci,aai,bbi,cci,ddi,eei,ffi,ggi,hhi,iii
      dimension dc(106),dci(106)
      common/coefot/dn(40,100),rn(100),d(100),drho(100),ap(100),dneg
     $     (40),an(20),jupn(40),jdn,jup
      complex*16 di,rni
      common/coefin/rni(100),di(100),jdni,jupi
      common/limits/kmax,kmaxr,kmax1,kmaxmm,kmaxr4,m,n,nm
      complex*16 egvi,ri1,rip1,xlami
      common/rwaves/egv,r1,rp1,r2,rp2,egvi,ri1,rip1,xlam(40),xlami(40)
      complex*16 uvo3,xyo3,uvi1,xyi1
      common/uvxy/uvo3(40,40),xyo3(40,40),uvo1(40,40),xyo1(40,40),
     $     uvi1(40,40),xyi1(40,40)
      complex*16 cfindx
      common/param/cfindx,rfindx,anglin,y,lsig
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
      kmax11=kmax+11
c
      jup6=jup+6
      jdn2=jdn+2
      do 5 l=1,jdn2
 5    dc(l)=0.0d0
      do 10 l=jdn,jup,2
 10   dc(l+4)=d(l)
      do 15 l=jup6,kmax11,2
 15   dc(l)=0.0d0
c 
      czero = (0.,0.)
      jupi6=jupi+6
      jdni2=jdni+2
      do 6 l=1,jdni2
 6    dci(l)=czero
      do 11 l=jdni,jupi,2
 11   dci(l+4)=di(l)
      do 16 l=jupi6,kmax11,2
 16   dci(l)=czero
c 
      do 40 kk=1,kmax1
         k=kk-1
c for inside calculations
         call ccaple(n,m,k,jdni,jupi,dci,aai,bbi,cci,ddi,
     $        eei,ffi,ggi,hhi,iii)
c for outside calculations
         call caplet(n,m,k,jdn,jup,dc,aa,bb,cc,dd,ee,ff,gg,hh,ii)
         if(mod(nm+k,2).ne.0) goto 25
cm* nm+k even
cm* do the outside calculations
         call uy(m,r1,rp1,r2,rp2,aa,bb,ee,hh,uo3,yo3)
         uvo3(nm+1,k+1)=uo3
         xyo3(nm+1,k+1)=yo3
         call uy(m,r1,rp1,0.0d0,0.0d0,aa,bb,ee,hh,uo1,yo1)
         uvo1(nm+1,k+1)=uo1
         xyo1(nm+1,k+1)=yo1
cm* do the inside evaluation
         call cuy(m,ri1,rip1,aai,bbi,eei,hhi,ui1,yi1)
cm? 
         ar=ri1
         if(cdabs(rip1).gt.cdabs(ri1)) ar=rip1
         uvi1(nm+1,k+1)=ui1/ar
         xyi1(nm+1,k+1)=yi1/ar
         goto 40
cm* nm+k odd
 25      call vx(m,r1,rp1,r2,rp2,egv,cc,dd,ff,gg,ii,vo3,xo3)
         uvo3(nm+1,k+1)=vo3
         xyo3(nm+1,k+1)=xo3
         call vx(m,r1,rp1,0.0d0,0.0d0,egv,cc,dd,ff,gg,ii,vo1,xo1)
         uvo1(nm+1,k+1)=vo1
         xyo1(nm+1,k+1)=xo1
         call cvx(m,ri1,rip1,egvi,cci,ddi,ffi,ggi,iii,vi1,xi1)
cm?
         ar=ri1
         if(cdabs(rip1).gt.cdabs(ri1)) ar=rip1
         xyi1(nm+1,k+1)=xi1/ar
         uvi1(nm+1,k+1)=vi1/ar
 40   continue
      return
      end
c
      subroutine caplet(n,m,k,jdn,jup,dc,aa,bb,cc,dd,ee,ff,gg,hh,ii)
      implicit real*8(a-h,o-z)
      real*8 m2,i1,i2,i3,i4,i5,i6,ii
      integer t1
      dimension dc(106)
      t1=k+1
      jup2=jup-2
      m2=dfloat(2*m)
      it = k
      t=dfloat(k)
      t2=2.0d0*t
      if(mod(n-m+k,2).eq.0) goto 200
      aa=0.0d0
      bb=0.0d0
      ee=0.0d0
      hh=0.0d0
      if(m.eq.0) goto 100
      c1=(t+m2+1.0d0)*dc(it+1+6)/(t2+m2+3.0d0)
      c2=t*dc(it-1+6)/(t2+m2-1.0d0)
      c3=(t+m2-1.0d0)*(t+m2)*(c1+c2)/(t2+m2+1.0d0)
      c4=(t+m2-1.0d0)*dc(it-1+6)/(t2+m2-1.0d0)
      c5=(t-2.0d0)*dc(it-3+6)/(t2+m2-5.0d0)
      c6=t*(t-1.0d0)*(c4+c5)/(t2+m2-3.0d0)
      cc=c3-c6
      sum=0.0d0
      do 15 j=t1,jup2,2
 15   sum=sum+dc(j+6)
      dd=t*dc(it-1+6)+(t2+m2-1.0d0)*sum
      f1=dc(it+1+6)/(t2+m2+3.0d0)-dc(it+3+6)/(t2+m2+7.0d0)
      f2=(t+m2-1.0d0)*(t+m2)*(t+m2+1.0d0)*(t+m2+2.0d0)*(t+m2+3.0d0)
      f3=(t2+m2+1.0d0)*(t2+m2+3.0d0)*(t2+m2+5.0d0)
      f4=f1*f2/f3
      f5=dc(it-1+6)/(t2+m2-1.0d0)-dc(it+1+6)/(t2+m2+3.0d0)
      f6=t*(t-m2)*(t+m2+1.0d0)*(t+m2)*(t+m2-1.0d0)
      f7=(t2+m2-3.0d0)*(t2+m2+1.0d0)*(t2+m2+3.0d0)
      f8=f5*f6/f7
      f9=dc(it-3+6)/(t2+m2-5.0d0)-dc(it-1+6)/(t2+m2-1.0d0)
      f10=t*(t-1.0d0)*(t-2.0d0)*(t+m2-1.0d0)*(t+2*m2-1.0d0)
      f11=(t2+m2-5.0d0)*(t2+m2-3.0d0)*(t2+m2+1.0d0)
      f12=f9*f10/f11
      f13=dc(it-5+6)/(t2+m2-9.0d0)-dc(it-3+6)/(t2+m2-5.0d0)
      f14=t*(t-1.0d0)*(t-2.0d0)*(t-3.0d0)*(t-4.0d0)
      f15=(t2+m2-3.0d0)*(t2+m2-5.0d0)*(t2+m2-7.0d0)
      f16=f13*f14/f15
      ff=f4-f8-f12+f16
      gg=-t*(t+m-1.0d0)*dc(it-1+6)+m*(t2+m2-1.0d0)*sum
      i1=(t+m+2.0d0)*(t+m2+1.0d0)*dc(it+1+6)/(t2+m2+3.0d0)
      i2=t*(m+t-1.0d0)*dc(it-1+6)/(t2+m2-1.0d0)
      i3=(t+m2)*(t+m2-1.0d0)*(i1-i2)/(t2+m2+1.0d0)
      i4=(t+m)*(t+m2-1.0d0)*dc(it-1+6)/(t2+m2-1.0d0)
      i5=(t-2.0d0)*(t+m-3.0d0)*dc(it-3+6)/(t2+m2-5.0d0)
      i6=t*(t-1.0d0)*(i4-i5)/(t2+m2-3.0d0)
      ii=i3-i6
      return
 100  c1=(t+1.0d0)*dc(it+1+6)/(t2+3.0d0)+t*dc(it-1+6)/(t2-1.0d0)
      c2=(t+3.0d0)*dc(it+3+6)/(t2+7.0d0)+(t+2.0d0)*dc(it+1+6)
     1     /(t2+3.0d0)
      cc=c1/(t2+1.0d0)-c2/(t2+5.0d0)
      dd=0.0d0
      f1=dc(it+1+6)/(t2+3.0d0)/(t2+5.0d0)
      f2=2.0d0*dc(it+3+6)/(t2+5.0d0)/(t2+9.0d0)
      f3=dc(it+5+6)/(t2+9.0d0)/(t2+11.0d0)
      f4=(t+3.0d0)*(t+4.0d0)*(t+5.0d0)*(f1-f2+f3)/(t2+5.0d0)/(t2+7.0d0)
      f5=dc(it-1+6)/(t2-1.0d0)/(t2+1.0d0)
      f6=2.0d0*dc(it+1+6)/(t2+1.0d0)/(t2+5.0d0)
      f7=dc(it+3+6)/(t2+5.0d0)/(t2+7.0d0)
      f8=3.0d0*t*(t+3.0d0)*(f5-f6+f7)/(t2+1.0d0)/(t2+5.0d0)
      f9=dc(it-3+6)/(t2-5.0d0)/(t2-3.0d0)
      f10=2.0d0*dc(it-1+6)/(t2-3.0d0)/(t2+1.0d0)
      f11=dc(it+1+6)/(t2+1.0d0)/(t2+3.0d0)
      f12=t*(t-1.0d0)*(t-2.0d0)*(f9-f10+f11)/(t2-1.0d0)/(t2+1.0d0)
      ff=f4+f8-f12
      gg=dc(it+1+6)
      i1=(t+1.0d0)*(t+2.0d0)*dc(it+1+6)/(t2+3.0d0)
      i2=t*(t-1.0d0)*dc(it-1+6)/(t2-1.0d0)
      i3=(i1-i2)/(t2+1.0d0)
      i4=(t+3.0d0)*(t+4.0d0)*dc(it+3+6)/(t2+7.0d0)
      i5=(t+1.0d0)*(t+2.0d0)*dc(it+1+6)/(t2+3.0d0)
      i6=(i4-i5)/(t2+5.0d0)
      ii=i3-i6
      return
 200  cc=0.0d0
      dd=0.0d0
      ff=0.0d0
      gg=0.0d0
      ii=0.0d0
      if(m.ne.0) goto 250
      aa=0.0d0
      bb=0.0d0
      ee=0.0d0
      hh=0.0d0
      return
 250  a1=(t+m2-1.0d0)*(t+m2)*dc(it+6)/(t2+m2+1.0d0)
      a2=t*(t-1.0d0)*dc(it-2+6)/(t2+m2-3.0d0)
      aa=a1-a2
      sum=0.0d0
      jup2=jup-2
      do 255 j=t1,jup2,2
255   sum=sum+dc(j-1+6)
      bb=(t2+m2-1.0d0)*sum
      e1=dc(it+6)/(t2+m2+1.0d0)-dc(it+2+6)/(t2+m2+5.0d0)
      e2=(t+m2-1.0d0)*(t+m2)*(t+m2+2.0d0)*(t+m2+1.0d0)
      e3=(t2+m2+1.0d0)*(t2+m2+3.0d0)
      e4=e1*e2/e3
      e5=dc(it-2+6)/(t2+m2-3.0d0)-dc(it+6)/(t2+m2+1.0d0)
      e6=2.0d0*t*(t-1.0d0)*(t+m2)*(t+m2-1.0d0)
      e7=(t2+m2-3.0d0)*(t2+m2+1.0d0)
      e8=e5*e6/e7
      e9=dc(it-4+6)/(t2+m2-7.0d0)-dc(it-2+6)/(t2+m2-3.0d0)
      e10=t*(t-1.0d0)*(t-2.0d0)*(t-3.0d0)
      e11=(t2+m2-3.0d0)*(t2+m2-5.0d0)
      e12=e9*e10/e11
      ee=e4-e8+e12
      h1=-t*(t-1.0d0)*(t+m -2.0d0)*dc(it-2+6)/(t2+m2-3.0d0)
      h2=-(t*(t-1.0d0)*(t2+m2+1.0d0)+(t+m2)*(t+m2-1.0d0))
     $     *dc(it+6)/2.0d0/(t2+m2+1.0d0)
      h3=m*(t2+m2-1.0d0)*(sum-dc(it+6))
      hh=h1+h2+h3
      return
      end
c
      subroutine ccaple(n,m,k,jdn,jup,dc,aa,bb,cc,dd,ee,ff,gg,hh,ii)
      implicit real*8(a-h,o-z)
cm* complex version of caplet
      real*8 m2
      complex*16 dc,aa,bb,cc,dd,ee,ff,gg,hh,ii
      dimension dc(106)
      complex*16 czero,sum
      complex*16 a1,a2
      complex*16 c1,c2,c3,c4,c5,c6
      complex*16 e1,e4,e5,e8,e9,e12
      complex*16 f1,f4,f5,f8,f9,f12,f13,f16,i1,i2,i3,i4,i5,i6
      complex*16 ff2,ff3,ff6,ff7,ff10,ff11
      complex*16 h1,h2,h3
      integer t1
      czero = (0.,0.)
      t1=k+1
      jup2=jup-2
      m2=dfloat(2*m)
      it = k
      t=dfloat(k)
      t2=2.0d0*t
      if(mod(n-m+k,2).eq.0) goto 200
      aa=czero
      bb=czero
      ee=czero
      hh=czero
      if(m.eq.0) goto 100
      c1=(t+m2+1.0d0)*dc(it+1+6)/(t2+m2+3.0d0)
      c2=t*dc(it-1+6)/(t2+m2-1.0d0)
      c3=(t+m2-1.0d0)*(t+m2)*(c1+c2)/(t2+m2+1.0d0)
      c4=(t+m2-1.0d0)*dc(it-1+6)/(t2+m2-1.0d0)
      c5=(t-2.0d0)*dc(it-3+6)/(t2+m2-5.0d0)
      c6=t*(t-1.0d0)*(c4+c5)/(t2+m2-3.0d0)
      cc=c3-c6
      sum=czero
      do 15 j=t1,jup2,2
 15   sum=sum+dc(j+6)
      dd=t*dc(it-1+6)+(t2+m2-1.0d0)*sum
      f1=dc(it+1+6)/(t2+m2+3.0d0)-dc(it+3+6)/(t2+m2+7.0d0)
      f2=(t+m2-1.0d0)*(t+m2)*(t+m2+1.0d0)*(t+m2+2.0d0)*(t+m2+3.0d0)
      f3=(t2+m2+1.0d0)*(t2+m2+3.0d0)*(t2+m2+5.0d0)
      f4=f1*f2/f3
      f5=dc(it-1+6)/(t2+m2-1.0d0)-dc(it+1+6)/(t2+m2+3.0d0)
      f6=t*(t-m2)*(t+m2+1.0d0)*(t+m2)*(t+m2-1.0d0)
      f7=(t2+m2-3.0d0)*(t2+m2+1.0d0)*(t2+m2+3.0d0)
      f8=f5*f6/f7
      f9=dc(it-3+6)/(t2+m2-5.0d0)-dc(it-1+6)/(t2+m2-1.0d0)
      f10=t*(t-1.0d0)*(t-2.0d0)*(t+m2-1.0d0)*(t+2*m2-1.0d0)
      f11=(t2+m2-5.0d0)*(t2+m2-3.0d0)*(t2+m2+1.0d0)
      f12=f9*f10/f11
      f13=dc(it-5+6)/(t2+m2-9.0d0)-dc(it-3+6)/(t2+m2-5.0d0)
      f14=t*(t-1.0d0)*(t-2.0d0)*(t-3.0d0)*(t-4.0d0)
      f15=(t2+m2-3.0d0)*(t2+m2-5.0d0)*(t2+m2-7.0d0)
      f16=f13*f14/f15
      ff=f4-f8-f12+f16
      gg=-t*(t+m-1.0d0)*dc(it-1+6)+m*(t2+m2-1.0d0)*sum
      i1=(t+m+2.0d0)*(t+m2+1.0d0)*dc(it+1+6)/(t2+m2+3.0d0)
      i2=t*(m+t-1.0d0)*dc(it-1+6)/(t2+m2-1.0d0)
      i3=(t+m2)*(t+m2-1.0d0)*(i1-i2)/(t2+m2+1.0d0)
      i4=(t+m)*(t+m2-1.0d0)*dc(it-1+6)/(t2+m2-1.0d0)
      i5=(t-2.0d0)*(t+m-3.0d0)*dc(it-3+6)/(t2+m2-5.0d0)
      i6=t*(t-1.0d0)*(i4-i5)/(t2+m2-3.0d0)
      ii=i3-i6
      return
 100  c1=(t+1.0d0)*dc(it+1+6)/(t2+3.0d0)+t*dc(it-1+6)/(t2-1.0d0)
      c2=(t+3.0d0)*dc(it+3+6)/(t2+7.0d0)+(t+2.0d0)*dc(it+1+6)
     1     /(t2+3.0d0)
      cc=c1/(t2+1.0d0)-c2/(t2+5.0d0)
      dd=czero
      f1=dc(it+1+6)/(t2+3.0d0)/(t2+5.0d0)
      ff2=2.0d0*dc(it+3+6)/(t2+5.0d0)/(t2+9.0d0)
      ff3=dc(it+5+6)/(t2+9.0d0)/(t2+11.0d0)
      f4=(t+3.0d0)*(t+4.0d0)*(t+5.0d0)*(f1-ff2+ff3)
     $     /(t2+5.0d0)/(t2+7.0d0)
      f5=dc(it-1+6)/(t2-1.0d0)/(t2+1.0d0)
      ff6=2.0d0*dc(it+1+6)/(t2+1.0d0)/(t2+5.0d0)
      ff7=dc(it+3+6)/(t2+5.0d0)/(t2+7.0d0)
      f8=3.0d0*t*(t+3.0d0)*(f5-ff6+ff7)/(t2+1.0d0)/(t2+5.0d0)
      f9=dc(it-3+6)/(t2-5.0d0)/(t2-3.0d0)
      ff10=2.0d0*dc(it-1+6)/(t2-3.0d0)/(t2+1.0d0)
      ff11=dc(it+1+6)/(t2+1.0d0)/(t2+3.0d0)
      f12=t*(t-1.0d0)*(t-2.0d0)*(f9-ff10+ff11)/(t2-1.0d0)/(t2+1.0d0)
      ff=f4+f8-f12
      gg=dc(it+1+6)
      i1=(t+1.0d0)*(t+2.0d0)*dc(it+1+6)/(t2+3.0d0)
      i2=t*(t-1.0d0)*dc(it-1+6)/(t2-1.0d0)
      i3=(i1-i2)/(t2+1.0d0)
      i4=(t+3.0d0)*(t+4.0d0)*dc(it+3+6)/(t2+7.0d0)
      i5=(t+1.0d0)*(t+2.0d0)*dc(it+1+6)/(t2+3.0d0)
      i6=(i4-i5)/(t2+5.0d0)
      ii=i3-i6
      return
 200  cc=czero
      dd=czero
      ff=czero
      gg=czero
      ii=czero
      if(m.ne.0) goto 250
      aa=czero
      bb=czero
      ee=czero
      hh=czero
      return
 250  a1=(t+m2-1.0d0)*(t+m2)*dc(it+6)/(t2+m2+1.0d0)
      a2=t*(t-1.0d0)*dc(it-2+6)/(t2+m2-3.0d0)
      aa=a1-a2
      sum=czero
      jup2=jup-2
      do 255 j=t1,jup2,2
 255  sum=sum+dc(j-1+6)
      bb=(t2+m2-1.0d0)*sum
      e1=dc(it+6)/(t2+m2+1.0d0)-dc(it+2+6)/(t2+m2+5.0d0)
      e2=(t+m2-1.0d0)*(t+m2)*(t+m2+2.0d0)*(t+m2+1.0d0)
      e3=(t2+m2+1.0d0)*(t2+m2+3.0d0)
      e4=e1*e2/e3
      e5=dc(it-2+6)/(t2+m2-3.0d0)-dc(it+6)/(t2+m2+1.0d0)
      e6=2.0d0*t*(t-1.0d0)*(t+m2)*(t+m2-1.0d0)
      e7=(t2+m2-3.0d0)*(t2+m2+1.0d0)
      e8=e5*e6/e7
      e9=dc(it-4+6)/(t2+m2-7.0d0)-dc(it-2+6)/(t2+m2-3.0d0)
      e10=t*(t-1.0d0)*(t-2.0d0)*(t-3.0d0)
      e11=(t2+m2-3.0d0)*(t2+m2-5.0d0)
      e12=e9*e10/e11
      ee=e4-e8+e12
      h1=-t*(t-1.0d0)*(t+m -2.0d0)*dc(it-2+6)/(t2+m2-3.0d0)
      h2=-(t*(t-1.0d0)*(t2+m2+1.0d0)+(t+m2)*(t+m2-1.0d0))
     $     *dc(it+6)/2.0d0/(t2+m2+1.0d0)
      h3=m*(t2+m2-1.0d0)*(sum-dc(it+6))
      hh=h1+h2+h3
      return
      end
c
      subroutine uy(m,r1,rp1,r2,rp2,aa,bb,ee,hh,u,y)
      implicit real*8(a-h,o-z)
      complex*16 u,y,y1,y2,icm
      complex*16 cin
      complex*16 dcmplx
      common/cincou/ cin,cout,sgn,zeta
cm* outside only - 2 calls
cm* scattered wave uses r(3) = r(1) + i r(2)
cm* note that u is real on call with r(1), and y is really pure imaginary
cm* but is changed with icm to store in the real part
      c = cout
      if(r1.ne.0.0d0) dlr1=dlog(dabs(r1))
      if(rp1.ne.0.0d0) dlp1=dlog(dabs(rp1))
      dzeta=(zeta*zeta-sgn)
      zeta2=dzeta*dzeta
      term=(zeta2*bb+sgn*2.0d0*dzeta*aa+ee)*dfloat(m)*zeta
      ur=0.0d0
      if(term.eq.0.0d0.or.r1.eq.0.0d0) goto 10
      if(dlr1+dlog(dabs(term)).gt.-175.0d0) ur=term*r1
 10   ui=term*r2
      u=dcmplx(ur,ui)
cm* calculate y
      icm=dcmplx(0.0d0,dfloat(m)/c)
      if(r2.eq.0.0d0.and.rp2.eq.0.0d0) icm=dcmplx(dfloat(m)/c,0.0d0)
      term1=sgn/dzeta*(aa+hh)
      y1r=0.0d0
      if(term1.eq.0.0d0.or.r1.eq.0.0d0) goto 20
      if(dlr1+dlog(dabs(term1)).gt.-175.0d0) y1r=term1*r1
 20   y1i=term1*r2
      y1=dcmplx(y1r,y1i)
      y2r1=0.0d0
      if(bb.eq.0.0d0.or.r1.eq.0.0d0) goto 30
      if(dlr1+dlog(dabs(bb)).gt.-175.0d0) y2r1=bb*r1
 30   y2i1=bb*r2
      term2=zeta*bb
      y2r2=0.0d0
      if(term2.eq.0.0d0.or.rp1.eq.0.0d0) goto 40
      if(dlp1+dlog(dabs(term2)).gt.-175.0d0) y2r2=term2*rp1
 40   y2i2=term2*rp2
      y2=dcmplx(y2r1+y2r2,y2i1+y2i2)
      y=icm*(y1+y2)
      return
      end
c
      subroutine cuy(m,r1,rp1,aa,bb,ee,hh,u,y)
      implicit real*8(a-h,o-z)
      complex*16 u,y,y1,y2,y21,y22,icm
      complex*16 c,czero,r1,rp1,aa,bb,ee,hh,term,term1,term2
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
cm* inside only
cm* all reference to r2, rp2 deleted
      c = cin
      czero = (0.,0.)
      if(r1.ne.czero) dlr1=dlog(cdabs(r1))
      if(rp1.ne.czero) dlp1=dlog(cdabs(rp1))
cm* calculate u
      dzeta=(zeta*zeta-sgn)
      zeta2=dzeta*dzeta
      term=(zeta2*bb+sgn*2.0d0*dzeta*aa+ee)*dfloat(m)*zeta
      u=czero
      if(term.eq.czero.or.r1.eq.czero) goto 10
      if(dlr1+dlog(cdabs(term)).gt.-175.0d0) u=term*r1
 10   continue
cm* calculate y
cm* note extraction of i, since r2,rp2=0.
      icm= dfloat(m)/c
      term1=sgn/dzeta*(aa+hh)
      y1=czero
      if(cdabs(term1).eq.czero.or.r1.eq.czero) goto 20
      if(dlr1+dlog(cdabs(term1)).gt.-175.0d0) y1=term1*r1
 20   continue
      y21=czero
      if(bb.eq.czero.or.r1.eq.czero) goto 30
      if(dlr1+dlog(cdabs(bb)).gt.-175.0d0) y21=bb*r1
 30   continue
      term2=zeta*bb
      y22=czero
      if(cdabs(term2).eq.czero.or.rp1.eq.czero) goto 40
      if(dlp1+dlog(cdabs(term2)).gt.-175.0d0) y22=term2*rp1
 40   continue
      y2=y21+y22
      y=icm*(y1+y2)
      return
      end
c
      subroutine vx(m,r1,rp1,r2,rp2,egv,cc,dd,ff,gg,ii,v,x)
      implicit real*8(a-h,o-z)
      real*8 ii
      complex*16 x,v,ic
      complex*16 cin
      complex*16 dcmplx
      common/cincou/ cin,cout,sgn,zeta
cm* this is for outside only - 2 calls
      c = cout
      if(r1.ne.0.0d0) dlr1=dlog(dabs(r1))
      if(rp1.ne.0.0d0) dlp1=dlog(dabs(rp1))
      term=zeta*gg
      xr1=0.0d0
      if(term.eq.0.0d0.or.r1.eq.0.0d0) goto 10
      if(dlr1+dlog(dabs(term)).gt.-175.0d0) xr1=term*r1
 10   xi1=term*r2
      term1=-sgn*cc
      xr2=0.0d0
      if(term1.eq.0.0d0.or.rp1.eq.0.0d0) goto 20
      if(dlp1+dlog(dabs(term1)).gt.-175.0d0) xr2=term1*rp1
 20   xi2=term1*rp2
      x=dcmplx(xr1+xr2,xi1+xi2)
      ic=dcmplx(0.0d0,1.0d0/c)
      if(r2.eq.0.0d0.and.rp2.eq.0.0d0) ic=dcmplx(1.0d0/c,0.0d0)
      dzeta=zeta*zeta-sgn
      zeta2=dzeta*dzeta
      term1=(zeta2*dd+sgn*2.0d0*dzeta*cc+ff)*dfloat(m*m)/dzeta*sgn
      term2=-sgn*((egv-(c*zeta)**2+sgn*dfloat(m*m)/dzeta)*(dzeta*cc+sgn
     $     *ff))
      term3=zeta*dzeta*(sgn*2.0d0*cc+dzeta*gg+sgn*ii)
      term4=zeta2*gg+sgn*(3.0d0*zeta*zeta-sgn)*ii
      term=term1+term2+term4
      v1r=0.0d0
      if(term.eq.0.0d0.or.r1.eq.0.0d0) goto 30
      if(dlr1+dlog(dabs(term)).gt.-175.0d0) v1r=term*r1
 30   v1i=term*r2
      v2r=0.0d0
      if(term3.eq.0.0d0.or.rp1.eq.0.0d0) goto 40
      if(dlp1+dlog(dabs(term3)).gt.-175.0d0)v2r=term3*rp1
 40   v2i=term3*rp2
      v=ic*dcmplx(v1r+v2r,v1i+v2i)
      return
      end
c
      subroutine cvx(m,r1,rp1,egv,cc,dd,ff,gg,ii,v,x)
      implicit real*8(a-h,o-z)
      complex*16 x,v,ic
      complex*16 c,czero,r1,rp1,egv,cc,dd,ff,gg,ii
      complex*16 term,term1,term2,term3,term4,x1,x2,v1,v2
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
cm* this is for inside only
cm* delete all reference to r2, rp2
      c = cin
      czero = (0.,0.)
      if(r1.ne.czero) dlr1=dlog(cdabs(r1))
      if(rp1.ne.czero) dlp1=dlog(cdabs(rp1))
cm* calculate x
      term=zeta*gg
      x1=czero
      if(term.eq.czero.or.r1.eq.czero) goto 10
      if(dlr1+dlog(cdabs(term)).gt.-175.0d0) x1=term*r1
 10   continue
      term1=-sgn*cc
      x2=czero
      if(term1.eq.czero.or.rp1.eq.czero) goto 20
      if(dlp1+dlog(cdabs(term1)).gt.-175.0d0) x2=term1*rp1
 20   continue
      x= x1+x2
cm* calculate v
cm* note extraction of i, since r2,rp2=0.
      ic= 1./c
      dzeta=zeta*zeta-sgn
      zeta2=dzeta*dzeta
      term1=(zeta2*dd+sgn*2.0d0*dzeta*cc+ff)*dfloat(m*m)/dzeta*sgn
      term2=-sgn*((egv-(c*zeta)**2+sgn*dfloat(m*m)/dzeta)*(dzeta*cc+sgn
     $     *ff))
      term4=zeta2*gg+sgn*(3.0d0*zeta*zeta-sgn)*ii
      term=term1+term2+term4
      v1=czero
      if(term.eq.czero.or.r1.eq.czero) goto 30
      if(dlr1+dlog(cdabs(term)).gt.-175.0d0) v1=term*r1
 30   continue
      term3=zeta*dzeta*(sgn*2.0d0*cc+dzeta*gg+sgn*ii)
      v2=czero
      if(term3.eq.czero.or.rp1.eq.czero) goto 40
      if(dlp1+dlog(cdabs(term3)).gt.-175.0d0) v2=term3*rp1
 40   continue
      v=ic* (v1+v2)
      return
      end
c
      subroutine fandg(m,anglin,d,jdn,jup,pl,nmax,f,g)
c   calculates the f's and g's as functions of angle of incidence
cm* outside only
      implicit real*8 (a-h,o-z)
      common /accura/ epslon,epslo2
      dimension d(jup),pl(120)
      m2=2*m
c scale factor for sum prevents overflows
      scale=10.0d0**(-m2)
      l=0
      if(m.eq.0) l=1
      if(jdn.eq.2) term=fact(m2)*scale
      if(jdn.eq.3) term=fact(m2+1)*scale
      cosang=dcos(anglin)
      sinang=dsin(anglin)
      sqr=sinang*sinang
      sum=0.0d0
      do 10 jj=jdn,jup,2
         j=jj-2
         if(dabs(d(jj)).le.epslo2) goto 10
         sum=sum+term*d(jj)*d(jj)/dfloat(2*j+m2+1)
         term=term*dfloat((m2+jj)*(m2+jj-1))/dfloat(jj*(jj-1))
 10   continue
      sum=2.0d0*sum/scale
      if(anglin.eq.0.0d0) goto 40
      sum1=0.0d0
      sum2=0.0d0
      if(m.eq.0) goto 35
      do 20 jj=jdn,jup,2
         if(d(jj).eq.0.0d0) goto 20
         j=jj-2
         term=pl(j+m+1)/sinang/dfloat((j+m)*(j+m+1))
         if(term .eq.0.0d0) goto 20
         if(dlog(dabs(d(jj)))+dlog(dabs(term)).le.-177.0d0) goto 20
         sum1=sum1+term*d(jj)
 20   continue
 35   js=jdn
      do 30 jj=js,jup,2
         j=jj-2
         if(d(jj).eq.0.0d0) goto 30
         dpl=(dfloat(j+1)*pl(j+m+2)-dfloat(j+m+1)*
     $        cosang*pl(j+m+1))/sinang
         if(m.eq.0.and.j.eq.0) goto 30
         if(dabs(dpl).le.epslon) goto 30
         term=dpl/dfloat((j+m)*(j+m+1))
         if(dlog(dabs(d(jj)))+dlog(dabs(term)).le.-177.0d0) goto 30
         sum2=sum2+term*d(jj)
 30   continue
      f=dfloat(4*m)*sum1/sum
      g=dfloat(2*(2-l))*sum2/sum
      return
 40   if(m.eq.1) goto 50
      f=0.0d0
      g=0.0d0
      return
 50   sum1=0.0d0
      do 55 jj=jdn,jup,2
 55   sum1=sum1+d(jj)
      f=2.0d0*sum1/sum
      g=f
      return
      end
c
      subroutine angl(m,jdn,jup,d,theta,p,nmax,chi,sigma)
c calculates modified spheroidal angular wavefunctions
cm* for use in extinction cross-sections
cm* and in tangle --- outside only
cm* p comes from elegen
      implicit real*8(a-h,o-z)
      common /accura/ epslon,epslo2
      dimension d(jup),p(120)
      logical even
      even=.false.
      if(mod(jdn,2).eq.0) even=.true.
      pi=3.141592653589793d0
      cosang=dcos(theta)
      sinang=dsin(theta)
      if(dabs(theta-pi/2.0d0).le.epslon) cosang=0.0d0
      if(dabs(theta-pi/2.0d0).le.epslon) sinang=1.0d0
      if(theta.le.epslon) cosang=1.0d0
      if(theta.le.epslon) sinang=0.0d0
      sum1=0.0d0
      sum2=0.0d0
      if(m.ne.1.and.(theta.eq.0.0d0.or.theta.eq.pi)) goto 15
      if(m.eq.1.and.(theta.eq.0.0d0.or.theta.eq.pi)) goto 20
      do 10 jj=jdn,jup,2
         j=jj-2
         if(d(jj).eq.0.0d0) goto 10
         if(p(j+m+1).eq.0.0d0) goto 5
         if(dlog(dabs(d(jj)))+dlog(dabs(p(j+m+1))).le.-177.0d0) goto 5
         sum1=sum1+d(jj)*p(j+m+1)
 5       if(even.and.dabs(theta-pi/2.0d0).le.epslon) goto 10
         dpl=(dfloat(j+1)*p(j+m+2)-dfloat(j+m+1)*
     $        cosang*p(j+m+1))/sinang
         if(dabs(dpl).le.epslon) goto 10
         if(dlog(dabs(dpl))+dlog(dabs(d(jj))).le.-177.0d0) goto 10
         sum2=sum2+d(jj)*dpl
 10   continue
      sum1=sum1/sinang
 15   sigma=m*sum1
      chi=sum2
      return
 20   do 30 jj=jdn,jup,2
         j=jj-2
 30   sum1=sum1+dfloat((j+1)*(j+2))*d(jj)
      chi=.50d0*sum1
      sigma=chi
      return
      end
c
      subroutine scater
c making coefficients for scattering cross sections
      implicit real*8(a-h,o-z)
      common/misc/ar2(40),f(40),g(40),chi(40),sigma(40),pie(40,40)
      common/coefot/dn(40,100),rn(100),d(100),drho(100),ap(100),dneg
     $     (40),an(20),jupn(40),jdn,jup
      common/limits/kmax,kmaxr,kmax1,kmaxmm,kmaxr4,m,n,nm
      m2=2*m
c scale factor for term,prevents overflows
      scale=10.0d0**(-m2)
      do 505 j1=1,kmax1
         j=j1-1
         do 510 jp1=j1,kmax1
            jp=jp1-1
            sum=0.0d0
            if(mod(jp-j,2).ne.0) goto 512
            k=2
            if(mod(j,2).ne.0) k=3
            jmin= min0(jupn(j+1),jupn(jp+1))
            term=fact(k-2+m2)*scale
            do 515 jj=k,jmin,2
               l=jj-2
               f3=term*2.0d0*dfloat((l+m)*(l+m+1))/dfloat(2*l+m2+1)
               term=term*dfloat((jj+m2)*(jj-1+m2))/dfloat(jj*(jj-1))
               if(f3.eq.0.0d0) goto 515
               if(dn(j+1,jj).eq.0.0d0.or.dn(jp+1,jj).eq.0.0d0) goto 515
               if(dlog(f3)+dlog(dabs(dn(j+1,jj)))+
     $              dlog(dabs(dn(jp+1,jj))).le.-177.0d0) goto 515
               sum=sum+dn(j+1,jj)*f3*dn(jp+1,jj)
 515        continue
 512        pie(j+1,jp+1)=sum/scale
            pie(jp+1,j+1)=sum/scale
 510     continue
 505  continue
      return
      end
c
      subroutine tangle
      implicit real*8(a-h,o-z)
      common/coefot/dn(40,100),rn(100),d(100),drho(100),ap(100),dneg
     $     (40),an(20),jupn(40),jdn,jup
      common/limits/kmax,kmaxr,kmax1,kmaxmm,kmaxr4,m,n,nm
      complex*16 cfindx
      common/param/cfindx,rfindx,anglin,y,lsig
      complex*16 t11,t12,t22,t21
      common/ampfun/tsigma(40,10),tchi(40,10),t11(36,7),t12(36,7),t22
     $     (36,7),t21(36,7)
      dimension tp(120),da(100)
      inc=10
      maxdeg=0
      do 10 j=1,kmax1
 10   maxdeg=max0(maxdeg,jupn(j)+m)
      pi=3.141592653589793d0
      rad=pi/180.0d0
      tinc=dfloat(inc)*rad
      lstop=90/inc+1
      do 60 ltheta=1,lstop
         theta=tinc*dfloat(ltheta-1)
         call elegen(maxdeg,m,theta,tp)
         do 70 nsm=1,kmax1
            jupa=jupn(nsm)
            jdna=mod(nsm-1,2)+2
            do 80 jj=jdna,jupa,2
 80         da(jj)=dn(nsm,jj)
            call angl(m,jdna,jupa,da,theta,tp,maxdeg,chie,sigmae)
            tsigma(nsm,ltheta)=sigmae
            tchi(nsm,ltheta)=chie
 70      continue
 60   continue
      return
      end
c
      subroutine temod1
c te mode matrix 1
c note that even alpha are combined with odd beta, and vice versa
      implicit real*8(a-h,o-z)
      complex*16 i
      complex*16 uvo3,xyo3,uvi1,xyi1
      common/uvxy/uvo3(40,40),xyo3(40,40),uvo1(40,40),xyo1(40,40),
     $     uvi1(40,40),xyi1(40,40)
      complex*16 a,b,alpha,beta
c
      complex*16 sa
      dimension sa(40)
c
      common/matrix/a(80,80),b(80),alpha(40),beta(40),matmax
      common/misc/ar2(40),f(40),g(40),chi(40),sigma(40),pie(40,40)
      complex*16 cexte1,cexte2,cextm3,cextm4,cscte2,csctm4
      common/xsectn/ cexte1,cexte2,cextm3,cextm4,cscte2,csctm4
      common/limits/kmax,kmaxr,kmax1,kmaxmm,kmaxr4,m,n,nm
      complex*16 cfindx
      complex*16 dcmplx
      common/param/cfindx,rfindx,anglin,y,lsig
      dimension rb(160)
      equivalence (b(1),rb(1))
      matmax=80
      i=(0.0d0,1.0d0)
      do 520 kks=1,kmaxr
         ks=kks-1
         k=2*ks
         l=4*ks
         b(l+1)=(0.0d0,0.0d0)
         b(l+2)=(0.0d0,0.0d0)
         b(l+3)=(0.0d0,0.0d0)
         b(l+4)=(0.0d0,0.0d0)
         do 530 jj=1,kmaxr
            j=jj-1
            j2=2*j
            j3=j2+1
            a(l+1,j+1)=i**(m+j2)*(uvo3(j2+1,k+1)+
     $           xyo3(j2+1,k+1))/ar2(j2+1)
            a(l+1,j+1+kmaxr)=i**(m+j3)*(uvo3(j3+1,k+1)+
     $           xyo3(j3+1,k+1))/ar2(j3+1)
            a(l+1,j+1+2*kmaxr)=-i**(m+j2)*(cfindx*uvi1(j2+1,k+1)+
     $           i*xyi1(j2+1,k+1))
            a(l+1,j+1+3*kmaxr)=-i**(m+j3)*(i*cfindx*uvi1(j3+1,k+1)+
     $           xyi1(j3+1,k+1))
            a(l+2,j+1)=i**(m+j2)*(uvo3(j2+1,k+1)-
     $           xyo3(j2+1,k+1))/ar2(j2+1)
            a(l+2,j+1+kmaxr)=i**(m+j3)*(uvo3(j3+1,k+1)-
     $           xyo3(j3+1,k+1))/ar2(j3+1)
            a(l+2,j+1+2*kmaxr)=-i**(m+j2)*(cfindx*uvi1(j2+1,k+1)-
     $           i*xyi1(j2+1,k+1))
            a(l+2,j+1+3*kmaxr)=-i**(m+j3)*(i*cfindx*uvi1(j3+1,k+1)-
     $           xyi1(j3+1,k+1))
            a(l+3,j+1)=i**(m+j2)*(uvo3(j2+1,k+2)+
     $           xyo3(j2+1,k+2))/ar2(j2+1)
            a(l+3,j+1+kmaxr)=i**(m+j3)*(uvo3(j3+1,k+2)+
     $           xyo3(j3+1,k+2))/ar2(j3+1)
            a(l+3,j+1+2*kmaxr)=-i**(m+j2)*(i*uvi1(j2+1,k+2)+
     $           cfindx*xyi1(j2+1,k+2))
            a(l+3,j+1+3*kmaxr)=-i**(m+j3)*(uvi1(j3+1,k+2)+
     $           i*cfindx*xyi1(j3+1,k+2))
            a(l+4,j+1)=i**(m+j2)*(uvo3(j2+1,k+2)-
     $           xyo3(j2+1,k+2))/ar2(j2+1)
            a(l+4,j+1+kmaxr)=i**(m+j3)*(uvo3(j3+1,k+2)-
     $           xyo3(j3+1,k+2))/ar2(j3+1)
            a(l+4,j+1+2*kmaxr)=-i**(m+j2)*(i*uvi1(j2+1,k+2)-
     $           cfindx*xyi1(j2+1,k+2))
            a(l+4,j+1+3*kmaxr)=-i**(m+j3)*(uvi1(j3+1,k+2)-
     $           i*cfindx*xyi1(j3+1,k+2))
            b(l+1)=(-i**(j2+m)*f(j2+1)*(uvo1(j2+1,k+1)+
     $           i*xyo1(j2+1,k+1))-i**(j3+m)*g(j3+1)*(i*uvo1(j3+1,k+1)+
     $           xyo1(j3+1,k+1)))+b(l+1)
            b(l+2)=(-i**(j2+m)*f(j2+1)*(uvo1(j2+1,k+1)-
     $           i*xyo1(j2+1,k+1))-i**(j3+m)*g(j3+1)*(i*uvo1(j3+1,k+1)-
     $           xyo1(j3+1,k+1)))+b(l+2)
            b(l+3)=(-i**(j2+m)*f(j2+1)*(i*uvo1(j2+1,k+2)+
     $           xyo1(j2+1,k+2))-i**(j3+m)*g(j3+1)*(uvo1(j3+1,k+2)+
     $           i*xyo1(j3+1,k+2)))+b(l+3)
            b(l+4)=(-i**(j2+m)*f(j2+1)*(i*uvo1(j2+1,k+2)-
     $           xyo1(j2+1,k+2))-i**(j3+m)*g(j3+1)*(uvo1(j3+1,k+2)-
     $           i*xyo1(j3+1,k+2)))+b(l+4)
 530     continue
 520  continue
      call cmateq(a,b,kmaxr4,1,matmax)
      do 526 j=1,kmax,2
         jj=(j+1)/2
         dar2=dabs(ar2(j))
         if(dar2.lt.1.0d65) goto 524
         jj2=2*jj
         ralpha=0.0d0
         if(dabs(rb(jj2-1)).ge.1.0d-76*dar2) ralpha=rb(jj2-1)/ar2(j)
         alphai=0.0d0
         if(dabs(rb(jj2)).ge.1.0d-76*dar2) alphai=rb(jj2)/ar2(j)
         alpha(j)=dcmplx(ralpha,alphai)
         goto 525
 524     alpha(j)=b(jj)/ar2(j)
 525     continue
         sa(j)=alpha(j)*sigma(j)
         cexte1=cexte1+sa(j)
 526  continue
      do 536 jj=1,kmaxr
         jk=jj+kmaxr
         j=2*jj
         dar2=dabs(ar2(j))
         if(dar2.lt.1.0d65) goto 534
         jk2=2*jk
         rbeta=0.0d0
         if(dabs(rb(jk2-1)).ge.1.0d-76*dar2) rbeta=rb(jk2-1)/ar2(j)
         betai=0.0d0
         if(dabs(rb(jk2)).ge.1.0d-76*dar2) betai=rb(jk2)/ar2(j)
         beta(j)=dcmplx(rbeta,betai)
         goto 535
 534     beta(j)=b(jk)/ar2(j)
 535     continue
         sa(j)=beta(j)*chi(j)
         cexte1=cexte1+sa(j)
 536  continue
c
      mode=1
      call conn(sa,cexte1,mode)
c
      return
      end
c
      subroutine conn(sa,cext,mode)
c convergence test within n series from given m, for extinction and phase
c this is used in 4 parts, 2 for te mode, 2 for tmmod
c comparison is of the last alpha+beta contribution to the current total cext
c cext is the current total over all m for the partial complex extinction
      implicit real*8(a-h,o-z)
      common/limits/kmax,kmaxr,kmax1,kmaxmm,kmaxr4,m,n,nm
      logical wanted
      common/output/ wanted
      common/conver/ cnvrgn, cnvrgs
      complex*16 sa,ccon,cext
      dimension sa(40)
c index for alpha(mn)
      moden=mod(mode+1,2)-1+m
      if(wanted) then
         do 555 jj=1,kmaxr
            j=2*jj-1
            write(6,556) j+moden,sa(j),sa(j+1)
 556        format(1x,i3,2(1pd11.3,1x,1pd11.3),2x)
 555     continue
      endif
      j=2*kmaxr-1
      ccon=sa(j)+sa(j+1)
cdec modification revealed by dec; probably not exact zero unless all sa are
      if(dreal(cext).ne.0.d0) then
         cconr=dreal(ccon)/dreal(cext)
      else
         cconr=0.d0
      endif
      if(dimag(cext).ne.0.d0) then
         cconi=dimag(ccon)/dimag(cext)
      else
         cconi=0.d0
      endif
      cri=dmax1(dabs(cconr),dabs(cconi))
      if(cri.gt.cnvrgn) write(6,557) mode,m,cext,cconr,cconi
 557  format(1x,'warning: n series cext convergence: mode, m',
     $        2i3,2(1pd11.3,1x,1pd11.3),2x)
      return
      end
c
      subroutine temod2
c te mode matrix 2 + scattering
c note that even alpha are combined with odd beta, and vice versa
      implicit real*8(a-h,o-z)
c
      complex*16 sa
      dimension sa(40)
c
      complex*16 i
      complex*16 uvo3,xyo3,uvi1,xyi1
      common/uvxy/uvo3(40,40),xyo3(40,40),uvo1(40,40),xyo1(40,40),
     $     uvi1(40,40),xyi1(40,40)
      complex*16 a,b,alpha,beta
      common/matrix/a(80,80),b(80),alpha(40),beta(40),matmax
      common/misc/ar2(40),f(40),g(40),chi(40),sigma(40),pie(40,40)
      complex*16 cexte1,cexte2,cextm3,cextm4,cscte2,csctm4
      common/xsectn/ cexte1,cexte2,cextm3,cextm4,cscte2,csctm4
      common/limits/kmax,kmaxr,kmax1,kmaxmm,kmaxr4,m,n,nm
      complex*16 cfindx
      complex*16 dconjg,dcmplx
      common/param/cfindx,rfindx,anglin,y,lsig
      logical wanted
      common/output/ wanted
      dimension rb(160)
      equivalence (b(1),rb(1))
      matmax=80
      i=(0.0d0,1.0d0)
      do 540 kks=1,kmaxr
         ks=kks-1
         k=2*ks
         l=4*ks
         b(l+1)=(0.0d0,0.0d0)
         b(l+2)=(0.0d0,0.0d0)
         b(l+3)=(0.0d0,0.0d0)
         b(l+4)=(0.0d0,0.0d0)
         do 550 jj=1,kmaxr
            j=jj-1
            j2=2*j
            j3=j2+1
            a(l+1,j+1)=i**(m+j3)*(uvo3(j3+1,k+2)+xyo3(j3+1,k+2))/
     $           ar2(j3+1)
            a(l+1,j+1+kmaxr)=i**(m+j2)*(uvo3(j2+1,k+2)+
     $           xyo3(j2+1,k+2))/ar2(j2+1)
            a(l+1,j+1+2*kmaxr)=-i**(m+j3)*(cfindx*uvi1(j3+1,k+2)+
     $           i*xyi1(j3+1,k+2))
            a(l+1,j+1+3*kmaxr)=-i**(m+j2)*(i*cfindx*uvi1(j2+1,k+2)+
     $           xyi1(j2+1,k+2))
            a(l+2,j+1)=i**(m+j3)*(uvo3(j3+1,k+2)-
     $           xyo3(j3+1,k+2))/ar2(j3+1)
            a(l+2,j+1+kmaxr)=i**(m+j2)*(uvo3(j2+1,k+2)-
     $           xyo3(j2+1,k+2))/ar2(j2+1)
            a(l+2,j+1+2*kmaxr)=-i**(m+j3)*(cfindx*uvi1(j3+1,k+2)-
     $           i*xyi1(j3+1,k+2))
            a(l+2,j+1+3*kmaxr)=-i**(m+j2)*(i*cfindx*uvi1(j2+1,k+2)-
     $           xyi1(j2+1,k+2))
            a(l+3,j+1)=i**(m+j3)*(uvo3(j3+1,k+1)+xyo3(j3+1,k+1))/
     $           ar2(j3+1)
            a(l+3,j+1+kmaxr)=i**(m+j2)*(uvo3(j2+1,k+1)+
     $           xyo3(j2+1,k+1))/ar2(j2+1)
            a(l+3,j+1+2*kmaxr)=-i**(m+j3)*(i*uvi1(j3+1,k+1)+
     $           cfindx*xyi1(j3+1,k+1))
            a(l+3,j+1+3*kmaxr)=-i**(m+j2)*(uvi1(j2+1,k+1)+
     $           i*cfindx*xyi1(j2+1,k+1))
            a(l+4,j+1)=i**(m+j3)*(uvo3(j3+1,k+1)-xyo3(j3+1,k+1))/
     $           ar2(j3+1)
            a(l+4,j+1+kmaxr)=i**(m+j2)*(uvo3(j2+1,k+1)-xyo3(j2+1,k+1))/
     $           ar2(j2+1)
            a(l+4,j+1+2*kmaxr)=-i**(m+j3)*(i*uvi1(j3+1,k+1)-
     $           cfindx*xyi1(j3+1,k+1))
            a(l+4,j+1+3*kmaxr)=-i**(m+j2)*(uvi1(j2+1,k+1)-
     $           i*cfindx*xyi1(j2+1,k+1))
            b(l+1)=(-i**(j3+m)*f(j3+1)*(uvo1(j3+1,k+2)+
     $           i*xyo1(j3+1,k+2))-i**(j2+m)*g(j2+1)*(i*uvo1(j2+1,k+2)+
     $           xyo1(j2+1,k+2)))+b(l+1)
            b(l+2)=(-i**(j3+m)*f(j3+1)*(uvo1(j3+1,k+2)-
     $           i*xyo1(j3+1,k+2))-i**(j2+m)*g(j2+1)*(i*uvo1(j2+1,k+2)-
     $           xyo1(j2+1,k+2)))+b(l+2)
            b(l+3)=(-i**(j3+m)*f(j3+1)*(i*uvo1(j3+1,k+1)+
     $           xyo1(j3+1,k+1))-i**(j2+m)*g(j2+1)*(uvo1(j2+1,k+1)+
     $           i*xyo1(j2+1,k+1)))+b(l+3)
            b(l+4)=(-i**(j3+m)*f(j3+1)*(i*uvo1(j3+1,k+1)-
     $           xyo1(j3+1,k+1))-i**(j2+m)*g(j2+1)*(uvo1(j2+1,k+1)-
     $           i*xyo1(j2+1,k+1)))+b(l+4)
 550     continue
 540  continue
      call cmateq(a,b,kmaxr4,1,matmax)
      do 5466 j=1,kmax,2
         jj=(j+1)/2+kmaxr
         dar2=dabs(ar2(j))
         if(dar2.lt.1.0d65) goto 544
         jj2=2*jj
         rbeta=0.0d0
         if(dabs(rb(jj2-1)).ge.1.0d-76*dar2) rbeta=rb(jj2-1)/ar2(j)
         betai=0.0d0
         if(dabs(rb(jj2)).ge.1.0d-76*dar2) betai=rb(jj2)/ar2(j)
         beta(j)=dcmplx(rbeta,betai)
         goto 545
 544     beta(j)=b(jj)/ar2(j)
 545     continue
         sa(j)=beta(j)*chi(j)
         cexte2=cexte2+sa(j)
 5466 continue
      do 556 jj=1,kmaxr
         j=2*jj
         dar2=dabs(ar2(j))
         if(dar2.lt.1.0d65) goto 554
         jj2=2*jj
         ralpha=0.0d0
         if(dabs(rb(jj2-1)).ge.1.0d-76*dar2) ralpha=rb(jj2-1)/ar2(j)
         alphai=0.0d0
         if(dabs(rb(jj2)).ge.1.0d-76*dar2) alphai=rb(jj2)/ar2(j)
         alpha(j)=dcmplx(ralpha,alphai)
         goto 555
 554     alpha(j)=b(jj)/ar2(j)
 555     continue
         sa(j)=alpha(j)*sigma(j)
         cexte2=cexte2+sa(j)
 556  continue
c
      mode=2
      call conn(sa,cexte2,mode)
c
c scattering
      do 547 j1=1,kmax1
         j=j1-1
         if(cdabs(alpha(j+1)).eq.0.0d0) goto 446
         do 557 jp1=1,kmax1
            jp=jp1-1
            if(cdabs(alpha(jp+1)).eq.0.0d0) goto 557
            if(pie(j+1,jp+1).eq.0.0d0) goto 557
            plog=dlog(dabs(pie(j+1,jp+1)))
            if(dlog(cdabs(alpha(j+1)))+dlog(cdabs(alpha(jp+1)))+
     $           plog.le.-100.0d0) goto 557
            cscte2=cscte2+pie(j+1,jp+1)*alpha(j+1)*dconjg(alpha(jp+1))
 557     continue
 446     if(cdabs(beta(j+1)).eq.0.0d0) goto  547
         do 457 jp1=1,kmax1
            jp=jp1-1
            if(cdabs(beta(jp+1)).eq.0.0d0) goto 457
            if(pie(j+1,jp+1).eq.0.0d0)goto 457
            plog=dlog(dabs(pie(j+1,jp+1)))
            if(dlog(cdabs(beta(j+1)))+dlog(cdabs(beta(jp+1)))+
     $           plog.le.-100.0d0) goto 457
            cscte2=cscte2+pie(j+1,jp+1)*beta(j+1)*dconjg(beta(jp+1))
 457     continue
 547  continue
c
      call conns(cscte2,mode)
c
      if(.not.wanted) goto 610
      write(6,601)
 601  format(1x,7hte mode,t27,5halpha,t74,4hbeta)
      do 546 j=1,kmax1
         write(6,548) alpha(j),beta(j)
 548     format(1x,4d24.16)
 546  continue
 610  continue
      return
      end
c
      subroutine conns(cscat,mode)
c term by term convergence check using last terms added in double series
c when increments of kmaxp are by 2
c this is very nonlinear, and changes in the alpha(nm) for small nm
c can be just as important, i.e.
c it appears as that as kmaxp (max n) increases, so does the solution for
c all of the alpha(nm), including the lower order ones already included.
c these changes in alpha (small n) are the same size as the whole new 
c alpha (high n) being added. thus the use of cnvrgn linearly for cext
c is a good strategy. however, the alpha's enter cscat in squared terms
c so that the effect of increasing n is to first order through the low n
c terms, but second order at the highest n. thus if we want to judge what
c the magnitude of the correction is on the fly, we must take the square 
c root of highest contribution (as normalized; dimensionless)
      implicit real*8(a-h,o-z)
      logical wanted
      common/output/ wanted
      common/limits/kmax,kmaxr,kmax1,kmaxmm,kmaxr4,m,n,nm
      complex*16 a,b,alpha,beta
      common/matrix/a(80,80),b(80),alpha(40),beta(40),matmax
      common/misc/ar2(40),f(40),g(40),chi(40),sigma(40),pie(40,40)
      complex*16 dconjg
      complex*16 cscat
      common/conver/ cnvrgn,cnvrgs
      cons=0.d0
      do 1547 j1=kmax1-1,kmax1
         j=j1-1
         do 1557 jp1=1,kmax1
            jp=jp1-1
            if(mod(iabs(jp-j),2).eq.1) go to 1557
            conss=dabs(dreal(pie(j+1,jp+1)*
     $           (alpha(j+1)*dconjg(alpha(jp+1))+
     $           beta(j+1)*dconjg(beta(jp+1)))))
            if(wanted) write(6,1559) j1,jp1,conss
 1559       format(1x,'cscat convergence diagnostic j, jp',
     $        2i3,2(1pd11.3,1x,1pd11.3),2x)
            if(conss.gt.cons) cons=conss
 1557     continue
 1547 continue
      cconr=cons/dreal(cscat)
c the kmax1 is a guess to get cnvrgs in same range as cnvrgn when each 1.d-4
c in the test case, cscat converged more slowly than cext
      cfact=dfloat(kmax1)*dsqrt(cconr)
      if(cfact.gt.cnvrgs) write(6,1558) mode,m,cscat,cconr,cfact
 1558 format(1x,
     $     ' warning: n series cscat conver, last elements: mode, m',
     $        2i3,2(1pd11.3,1x,1pd11.3),2x)
c
      return
      end

      subroutine teamp
c calculation of amplitude functions t11,t12
      implicit real*8(a-h,o-z)
      common /accura/ epslon,epslo2
      complex*16 a,b,alpha,beta
      complex*16 t11,t12,t22,t21
      complex*16 t1term,t2term
      common/matrix/a(80,80),b(80),alpha(40),beta(40),matmax
      common/limits/kmax,kmaxr,kmax1,kmaxmm,kmaxr4,m,n,nm
      complex*16 cfindx
      common/param/cfindx,rfindx,anglin,y,lsig
      common/ampfun/tsigma(40,10),tchi(40,10),t11(36,7),t12(36,7),t22
     $     (36,7),t21(36,7)
      logical flag
      pi=3.141592653589793d0
      rad=pi/180.0d0
      inceta=10
      incphi=15
      uphi=rad*dfloat(incphi)
      lstop=90/inceta+1
      limphi=90/incphi+1
      flag=.false.
      if(m.eq.0) flag=.true.
      if(m.eq.1.and.anglin.eq.0.0d0) flag=.true.
      if(.not.flag) goto3
      jset=360/inceta
      do 2 j=1,jset
         do 1 k=1,limphi
            t11(j,k)=(0.0d0,0.0d0)
 1       t12(j,k)=(0.0d0,0.0d0)
 2    continue
 3    continue
      if(anglin.eq.0.0d0) limphi=1
      do 10 iphi=1,limphi
         phi=uphi*dfloat(iphi-1)
         phim=dfloat(m)*phi
         cosphm=dcos(phim)
         sinphm=dsin(phim)
         if(dabs(cosphm).le.epslon) cosphm=0.0d0
         if(dabs(sinphm).le.epslon) sinphm=0.0d0
         jstop=4
         kstart=1
         kstop=lstop
         lunit=lstop-1
       do 20 jquad=1,jstop
          if(jquad.eq.4) kstop=lstop-1
          do 30 ktheta=kstart,kstop
             ltheta=(jquad-1)*lunit+ktheta
             mtheta=ktheta
             if(jquad.eq.2) mtheta=lstop+1-ktheta
             if(jquad.eq.4) mtheta=lstop+1-ktheta
             fsig=1.0d0
             fchi=1.0d0
             if(jquad.eq.4) fsig=dfloat((-1)**m)
             if(jquad.eq.4) fchi=dfloat((-1)**m)
             do 40 nsm=1,kmax1
                nm=nsm-1
                if(jquad.eq.2) fsig=dfloat((-1)**nm)
                if(jquad.eq.2) fchi=dfloat((-1)**nsm)
                if(jquad.eq.3) fsig=dfloat((-1)**(nm+m))
                if(jquad.eq.3) fchi=dfloat((-1)**(nsm+m))
                sigmae=fsig*tsigma(nsm,mtheta)
                chie=fchi*tchi(nsm,mtheta)
                t1term=alpha(nsm)*sigmae+beta(nsm)*chie
                t11(ltheta,iphi)=t11(ltheta,iphi)+t1term*cosphm
                t2term=alpha(nsm)*chie+beta(nsm)*sigmae
c  minus sign appears for t12,mistake in japanese paper
                t12(ltheta,iphi)=t12(ltheta,iphi)-t2term*sinphm
 40          continue
 30       continue
          kstart=2
 20    continue
 10   continue
      return
      end
c
      subroutine tmmod3
c tm mode matrix 3
      implicit real*8(a-h,o-z)
c
      complex*16 sa
      dimension sa(40)
c
      complex*16 i
      complex*16 uvo3,xyo3,uvi1,xyi1
      common/uvxy/uvo3(40,40),xyo3(40,40),uvo1(40,40),xyo1(40,40),
     $     uvi1(40,40),xyi1(40,40)
      complex*16 a,b,alpha,beta
      common/matrix/a(80,80),b(80),alpha(40),beta(40),matmax
      common/misc/ar2(40),f(40),g(40),chi(40),sigma(40),pie(40,40)
      common/limits/kmax,kmaxr,kmax1,kmaxmm,kmaxr4,m,n,nm
      complex*16 cexte1,cexte2,cextm3,cextm4,cscte2,csctm4
      common/xsectn/ cexte1,cexte2,cextm3,cextm4,cscte2,csctm4
      complex*16 cfindx
      complex*16 dcmplx
      common/param/cfindx,rfindx,anglin,y,lsig
      dimension rb(160)
      equivalence (b(1),rb(1))
      matmax=80
      i=(0.0d0,1.0d0)
      do 560 kks=1,kmaxr
         ks=kks-1
         k=2*ks
         l=4*ks
         b(l+1)=(0.0d0,0.0d0)
         b(l+2)=(0.0d0,0.0d0)
         b(l+3)=(0.0d0,0.0d0)
         b(l+4)=(0.0d0,0.0d0)
         do 570 jj=1,kmaxr
            j=jj-1
            j2=2*j
            j3=j2+1
            a(l+1,j+1)=i**(m+j2)*(uvo3(j2+1,k+1)+xyo3(j2+1,k+1))/
     $           ar2(j2+1)
            a(l+1,j+1+kmaxr)=i**(m+j3)*(uvo3(j3+1,k+1)+
     $           xyo3(j3+1,k+1))/ar2(j3+1)
            a(l+1,j+1+2*kmaxr)=-i**(m+j2)*(uvi1(j2+1,k+1)+
     $           i*cfindx*xyi1(j2+1,k+1))
            a(l+1,j+1+3*kmaxr)=-i**(m+j3)*(i*uvi1(j3+1,k+1)+
     $           cfindx*xyi1(j3+1,k+1))
            a(l+2,j+1)=i**(m+j2)*(uvo3(j2+1,k+1)-
     $           xyo3(j2+1,k+1))/ar2(j2+1)
            a(l+2,j+1+kmaxr)=i**(m+j3)*(uvo3(j3+1,k+1)-
     $           xyo3(j3+1,k+1))/ar2(j3+1)
            a(l+2,j+1+2*kmaxr)=-i**(m+j2)*(uvi1(j2+1,k+1)-
     $           i*cfindx*xyi1(j2+1,k+1))
            a(l+2,j+1+3*kmaxr)=-i**(m+j3)*(i*uvi1(j3+1,k+1)-
     $           cfindx*xyi1(j3+1,k+1))
            a(l+3,j+1)=i**(m+j2)*(uvo3(j2+1,k+2)+xyo3(j2+1,k+2))/
     $           ar2(j2+1)
            a(l+3,j+1+kmaxr)=i**(m+j3)*(uvo3(j3+1,k+2)+
     $           xyo3(j3+1,k+2))/ar2(j3+1)
            a(l+3,j+1+2*kmaxr)=-i**(m+j2)*(i*cfindx*uvi1(j2+1,k+2)+
     $           xyi1(j2+1,k+2))
            a(l+3,j+1+3*kmaxr)=-i**(m+j3)*(cfindx*uvi1(j3+1,k+2)+
     $           i*xyi1(j3+1,k+2))
            a(l+4,j+1)=i**(m+j2)*(uvo3(j2+1,k+2)-xyo3(j2+1,k+2))/
     $           ar2(j2+1)
            a(l+4,j+1+kmaxr)=i**(m+j3)*(uvo3(j3+1,k+2)-
     $           xyo3(j3+1,k+2))/ar2(j3+1)
            a(l+4,j+1+2*kmaxr)=-i**(m+j2)*(i*cfindx*uvi1(j2+1,k+2)-
     $           xyi1(j2+1,k+2))
            a(l+4,j+1+3*kmaxr)=-i**(m+j3)*(cfindx*uvi1(j3+1,k+2)-
     $           i*xyi1(j3+1,k+2))
            b(l+1)=(-i**(j2+m)*f(j2+1)*(uvo1(j2+1,k+1)+
     $           i*xyo1(j2+1,k+1))-i**(j3+m)*g(j3+1)*(i*uvo1(j3+1,k+1)+
     $           xyo1(j3+1,k+1)))+b(l+1)
            b(l+2)=(-i**(j2+m)*f(j2+1)*(uvo1(j2+1,k+1)-
     $           i*xyo1(j2+1,k+1))-i**(j3+m)*g(j3+1)*(i*uvo1(j3+1,k+1)-
     $           xyo1(j3+1,k+1)))+b(l+2)
            b(l+3)=(-i**(j2+m)*f(j2+1)*(i*uvo1(j2+1,k+2)+
     $           xyo1(j2+1,k+2))-i**(j3+m)*g(j3+1)*(uvo1(j3+1,k+2)+
     $           i*xyo1(j3+1,k+2)))+b(l+3)
            b(l+4)=(-i**(j2+m)*f(j2+1)*(i*uvo1(j2+1,k+2)-
     $           xyo1(j2+1,k+2))-i**(j3+m)*g(j3+1)*(uvo1(j3+1,k+2)-
     $           i*xyo1(j3+1,k+2)))+b(l+4)
 570     continue
 560  continue
      call cmateq(a,b,kmaxr4,1,matmax)
      do 566 j=1,kmax,2
         jj=(j+1)/2
         dar2=dabs(ar2(j))
         if(dar2.lt.1.0d65) goto 564
         jj2=2*jj
         ralpha=0.0d0
         if(dabs(rb(jj2-1)).ge.1.0d-76*dar2) ralpha=rb(jj2-1)/ar2(j)
         alphai=0.0d0
         if(dabs(rb(jj2)).ge.1.0d-76*dar2) alphai=rb(jj2)/ar2(j)
         alpha(j)=dcmplx(ralpha,alphai)
         goto 565
 564     alpha(j)=b(jj)/ar2(j)
 565     continue
         sa(j)=alpha(j)*sigma(j)
         cextm3=cextm3+sa(j)
 566  continue
      do 576 jj=1,kmaxr
         jk=jj+kmaxr
         j=2*jj
         dar2=dabs(ar2(j))
         if(dar2.lt.1.0d65) goto 574
         jk2=2*jk
         rbeta=0.0d0
         if(dabs(rb(jk2-1)).ge.1.0d-76*dar2) rbeta=rb(jk2-1)/ar2(j)
         betai=0.0d0
         if(dabs(rb(jk2)).ge.1.0d-76*dar2) betai=rb(jk2)/ar2(j)
         beta(j)=dcmplx(rbeta,betai)
         goto 575
 574     beta(j)=b(jk)/ar2(j)
 575     continue
         sa(j)=beta(j)*chi(j)
         cextm3=cextm3+sa(j)
 576  continue
c
      mode=3
      call conn(sa,cextm3,mode)
c
      return
      end
c
      subroutine tmmod4
c tm mode matrix 4
      implicit real*8(a-h,o-z)
c
      complex*16 sa
      dimension sa(40)
c
      complex*16 i
      complex*16 uvo3,xyo3,uvi1,xyi1
      common/uvxy/uvo3(40,40),xyo3(40,40),uvo1(40,40),xyo1(40,40),
     $     uvi1(40,40),xyi1(40,40)
      complex*16 a,b,alpha,beta
      common/matrix/a(80,80),b(80),alpha(40),beta(40),matmax
      common/misc/ar2(40),f(40),g(40),chi(40),sigma(40),pie(40,40)
      common/limits/kmax,kmaxr,kmax1,kmaxmm,kmaxr4,m,n,nm
      complex*16 cexte1,cexte2,cextm3,cextm4,cscte2,csctm4
      common/xsectn/ cexte1,cexte2,cextm3,cextm4,cscte2,csctm4
      logical wanted
      common/output/ wanted
      complex*16 cfindx
      complex*16 dconjg,dcmplx
      common/param/cfindx,rfindx,anglin,y,lsig
      dimension rb(160)
      equivalence (b(1),rb(1))
      matmax=80
      i=(0.0d0,1.0d0)
      do 580 kks=1,kmaxr
         ks=kks-1
         k=2*ks
         l=4*ks
         b(l+1)=(0.0d0,0.0d0)
         b(l+2)=(0.0d0,0.0d0)
         b(l+3)=(0.0d0,0.0d0)
         b(l+4)=(0.0d0,0.0d0)
         do 590 jj=1,kmaxr
            j=jj-1
            j2=2*j
            j3=j2+1
            a(l+1,j+1)=i**(m+j3)*(uvo3(j3+1,k+2)+xyo3(j3+1,k+2))/
     $           ar2(j3+1)
            a(l+1,j+1+kmaxr)=i**(m+j2)*(uvo3(j2+1,k+2)+
     $           xyo3(j2+1,k+2))/ar2(j2+1)
            a(l+1,j+1+2*kmaxr)=-i**(m+j3)*(uvi1(j3+1,k+2)+
     $           i*cfindx*xyi1(j3+1,k+2))
            a(l+1,j+1+3*kmaxr)=-i**(m+j2)*(i*uvi1(j2+1,k+2)+
     $           cfindx*xyi1(j2+1,k+2))
            a(l+2,j+1)=i**(m+j3)*(uvo3(j3+1,k+2)-xyo3(j3+1,k+2))/
     $           ar2(j3+1)
            a(l+2,j+1+kmaxr)=i**(m+j2)*(uvo3(j2+1,k+2)-
     $           xyo3(j2+1,k+2))/ar2(j2+1)
            a(l+2,j+1+2*kmaxr)=-i**(m+j3)*(uvi1(j3+1,k+2)-
     $           i*cfindx*xyi1(j3+1,k+2))
            a(l+2,j+1+3*kmaxr)=-i**(m+j2)*(i*uvi1(j2+1,k+2)-
     $           cfindx*xyi1(j2+1,k+2))
            a(l+3,j+1)=i**(m+j3)*(uvo3(j3+1,k+1)+xyo3(j3+1,k+1))/
     $           ar2(j3+1)
            a(l+3,j+1+kmaxr)=i**(m+j2)*(uvo3(j2+1,k+1)+
     $           xyo3(j2+1,k+1))/ar2(j2+1)
            a(l+3,j+1+2*kmaxr)=-i**(m+j3)*(i*cfindx*uvi1(j3+1,k+1)+
     $           xyi1(j3+1,k+1))
            a(l+3,j+1+3*kmaxr)=-i**(m+j2)*(cfindx*uvi1(j2+1,k+1)+
     $           i*xyi1(j2+1,k+1))
            a(l+4,j+1)=i**(m+j3)*(uvo3(j3+1,k+1)-xyo3(j3+1,k+1))/
     $           ar2(j3+1)
            a(l+4,j+1+kmaxr)=i**(m+j2)*(uvo3(j2+1,k+1)-
     $           xyo3(j2+1,k+1))/ar2(j2+1)
            a(l+4,j+1+2*kmaxr)=-i**(m+j3)*(i*cfindx*uvi1(j3+1,k+1)-
     $           xyi1(j3+1,k+1))
            a(l+4,j+1+3*kmaxr)=-i**(m+j2)*(cfindx*uvi1(j2+1,k+1)-
     $           i*xyi1(j2+1,k+1))
            b(l+1)=(-i**(j3+m)*f(j3+1)*(uvo1(j3+1,k+2)+
     $           i*xyo1(j3+1,k+2))-i**(j2+m)*g(j2+1)*(i*uvo1(j2+1,k+2)+
     $           xyo1(j2+1,k+2)))+b(l+1)
            b(l+2)=(-i**(j3+m)*f(j3+1)*(uvo1(j3+1,k+2)-
     $           i*xyo1(j3+1,k+2))-i**(j2+m)*g(j2+1)*(i*uvo1(j2+1,k+2)-
     $           xyo1(j2+1,k+2)))+b(l+2)
            b(l+3)=(-i**(j3+m)*f(j3+1)*(i*uvo1(j3+1,k+1)+
     $           xyo1(j3+1,k+1))-i**(j2+m)*g(j2+1)*(uvo1(j2+1,k+1)+
     $           i*xyo1(j2+1,k+1)))+b(l+3)
            b(l+4)=(-i**(j3+m)*f(j3+1)*(i*uvo1(j3+1,k+1)-
     $           xyo1(j3+1,k+1))-i**(j2+m)*g(j2+1)*(uvo1(j2+1,k+1)-
     $           i*xyo1(j2+1,k+1)))+b(l+4)
 590     continue
 580  continue
      call cmateq(a,b,kmaxr4,1,matmax)
      do 586 j=1,kmax,2
         jj=(j+1)/2 +kmaxr
         dar2=dabs(ar2(j))
         if(dar2.lt.1.0d65) goto 584
         jj2=2*jj
         rbeta=0.0d0
         if(dabs(rb(jj2-1)).ge.1.0d-76*dar2) rbeta=rb(jj2-1)/ar2(j)
         betai=0.0d0
         if(dabs(rb(jj2)).ge.1.0d-76*dar2) betai=rb(jj2)/ar2(j)
         beta(j)=dcmplx(rbeta,betai)
         goto 585
 584     beta(j)=b(jj)/ar2(j)
 585     continue
         sa(j)=beta(j)*chi(j)
         cextm4=cextm4+sa(j)
 586  continue
      do 596 jj=1,kmaxr
         j=2*jj
         dar2=dabs(ar2(j))
         if(dar2.lt.1.0d65) goto 594
         jj2=2*jj
         ralpha=0.0d0
         if(dabs(rb(jj2-1)).ge.1.0d-76*dar2) ralpha=rb(jj2-1)/ar2(j)
         alphai=0.0d0
         if(dabs(rb(jj2)).ge.1.0d-76*dar2) alphai=rb(jj2)/ar2(j)
         alpha(j)=dcmplx(ralpha,alphai)
         goto 595
 594     alpha(j)=b(jj)/ar2(j)
 595     continue
         sa(j)=alpha(j)*sigma(j)
         cextm4=cextm4+sa(j)
 596  continue
c
      mode=4
      call conn(sa,cextm4,mode)
c
      do 587 j1=1,kmax1
         j=j1-1
         if(cdabs(alpha(j+1)).eq.0.0d0) goto 486
         do 597 jp1=1,kmax1
            jp=jp1-1
            if(cdabs(alpha(jp+1)).eq.0.0d0) goto 597
            if(pie(j+1,jp+1).eq.0.0d0) goto 597
            plog=dlog(dabs(pie(j+1,jp+1)))
            if(dlog(cdabs(alpha(j+1)))+dlog(cdabs(alpha(jp+1)))+
     $           plog.le.-100.0d0) goto 597
            csctm4=csctm4+pie(j+1,jp+1)* alpha(j+1)*dconjg(alpha(jp+1))
 597     continue
 486     if(cdabs(beta(j+1)).eq.0.0d0) goto 587
         do 497 jp1=1,kmax1
            jp=jp1-1
            if(cdabs(beta(jp+1)).eq.0.0d0) goto 497
            if(pie(j+1,jp+1).eq.0.0d0) goto 497
            plog=dlog(dabs(pie(j+1,jp+1)))
            if(dlog(cdabs(beta(j+1)))+dlog(cdabs(beta(jp+1)))+
     $           plog.le.-100.0d0) goto 497
            csctm4=csctm4+pie(j+1,jp+1)*beta(j+1)*dconjg(beta(jp+1))
 497     continue
 587  continue
c
      call conns(csctm4,mode)
c
      if(.not.wanted) goto 610
      write(6,601)
 601  format(1x,7htm mode,t27,5halpha,t74,4hbeta)
      do 599 j=1,kmax1
         write(6,598) alpha(j),beta(j)
 598     format(1x,4d24.16)
 599  continue
 610  continue
      return
      end
c
      subroutine cmateq(a,b,iii,jjj,id)
      implicit real*8(a-h,o-z)
      complex*16 a,b,s
      dimension a(id,id),b(id,1)
      kk=iii
      nv=jjj
      kkm=kk-1
      do 9 i=1,kkm
         s=(0.0d0,0.0d0)
         do 1 j=i,kk
            r=cdabs(a(j,i))
            if(r.lt.cdabs(s)) goto 1
            s=r
            l=j
 1       continue
         if(l.eq.i) goto 5
         do 2 j=i,kk
            s=a(i,j)
            a(i,j)=a(l,j)
 2       a(l,j)=s
         do 3 j=1,nv
            s=b(i,j)
            b(i,j)=b(l,j)
 3       b(l,j)=s
 5       if(cdabs(a(i,i)).eq.0.0d0) goto 9
         ipo=i+1
         do 8 j=ipo,kk
            if(cdabs(a(j,i)).eq.0.0d0) goto 8
            s=a(j,i)/a(i,i)
            a(j,i)=(0.0d0,0.0d0)
            do 6 k=ipo,kk
 6          a(j,k)=a(j,k)-a(i,k)*s
            if(nv.le.0) goto 8
            do 7 k=1,nv
 7          b(j,k)=b(j,k)-b(i,k)*s
 8       continue
 9    continue
      kmo=kk-1
      do 12  k=1,nv
         b(kk,k)=b(kk,k)/a(kk,kk)
         do 12 i=1,kmo
            n=kk-i
            do 11 j=n,kmo
 11         b(n,k)=b(n,k)-a(n,j+1)*b(j+1,k)
 12   b(n,k)=b(n,k)/a(n,n)
      return
      end
c
      subroutine tmamp
c calculation of amplitude functions t22,t21
      implicit real*8(a-h,o-z)
      common /accura/ epslon,epslo2
      complex*16 a,b,alpha,beta
      complex*16 t11,t12,t22,t21
      complex*16 t1term,t2term
      common/matrix/a(80,80),b(80),alpha(40),beta(40),matmax
      common/limits/kmax,kmaxr,kmax1,kmaxmm,kmaxr4,m,n,nm
      complex*16 cfindx
      common/param/cfindx,rfindx,anglin,y,lsig
      common/ampfun/tsigma(40,10),tchi(40,10),t11(36,7),t12(36,7),t22
     $     (36,7),t21(36,7)
      logical flag
      pi=3.141592653589793d0
      rad=pi/180.0d0
      inceta=10
      incphi=15
      uphi=rad*dfloat(incphi)
      lstop=90/inceta+1
      limphi=90/incphi+1
      flag=.false.
      if(m.eq.0) flag=.true.
      if(m.eq.1.and.anglin.eq.0.0d0) flag=.true.
      if(.not.flag) goto3
      jset=360/inceta
      do 2 j=1,jset
         do 1 k=1,limphi
            t22(j,k)=(0.0d0,0.0d0)
 1       t21(j,k)=(0.0d0,0.0d0)
 2    continue
 3    continue
      if(anglin.eq.0.0d0) limphi=1
      do 10 iphi=1,limphi
         phi=uphi*dfloat(iphi-1)
         phim=dfloat(m)*phi
         cosphm=dcos(phim)
         sinphm=dsin(phim)
         if(dabs(cosphm).le.epslon) cosphm=0.0d0
         if(dabs(sinphm).le.epslon) sinphm=0.0d0
         jstop=4
         kstart=1
         kstop=lstop
         lunit=lstop-1
         do 20 jquad=1,jstop
            if(jquad.eq.4) kstop=lstop-1
            do 30 ktheta=kstart,kstop
               ltheta=(jquad-1)*lunit+ktheta
               mtheta=ktheta
               if(jquad.eq.2) mtheta=lstop+1-ktheta
               if(jquad.eq.4) mtheta=lstop+1-ktheta
               fsig=1.0d0
               fchi=1.0d0
               if(jquad.eq.4) fsig=dfloat((-1)**m)
               if(jquad.eq.4) fchi=dfloat((-1)**m)
               do 40 nsm=1,kmax1
                  nm=nsm-1
                  if(jquad.eq.2) fsig=dfloat((-1)**nm)
                  if(jquad.eq.2) fchi=dfloat((-1)**nsm)
                  if(jquad.eq.3) fsig=dfloat((-1)**(nm+m))
                  if(jquad.eq.3) fchi=dfloat((-1)**(nsm+m))
                  sigmae=fsig*tsigma(nsm,mtheta)
                  chie=fchi*tchi(nsm,mtheta)
                  t2term=alpha(nsm)*sigmae+beta(nsm)*chie
                  t22(ltheta,iphi)=t22(ltheta,iphi)+t2term*cosphm
                  t1term=alpha(nsm)*chie+beta(nsm)*sigmae
                  t21(ltheta,iphi)=t21(ltheta,iphi)+t1term*sinphm
 40            continue
 30         continue
            kstart=2
 20      continue
 10   continue
      return
      end
c
      subroutine trans
c calculates elements of the transformation matrix for the stokes
c parameters as defined in van de hulst
      implicit real*8(a-h,o-z)
      complex*16 cfindx
      complex*16 dconjg
      common/param/cfindx,rfindx,anglin,y,lsig
      complex*16 t11,t12,t22,t21
      common/ampfun/tsigma(40,10),tchi(40,10),t11(36,7),t12(36,7),t22
     $     (36,7),t21(36,7)
      complex*16 i,a
      real*8 m(4)
      dimension a(4),s(4,4),d(4,4),f(4,4)
      i=(0.0d0,1.0d0)
      pi=3.141592653589793d0
      rad=pi/180.0d0
      incphi=15
      inceta=10
      limphi=90/incphi+1
      iplane=180/incphi
      lstop=360/inceta
      lhalf=lstop/2+1
      if(anglin.eq.0.0d0) limphi=1
      lcount=0
      do 10 iphi=1,limphi
         do 20 ltheta=1,lstop
            if(lcount.eq.0) write(6,99)
 99         format('1',1x,3hphi,2x,5htheta,2x,10hscattering,
     $           t40,8h!m!(1,4),t70,5hphase,
     $           t105,3h!f!/16x,5hangle,t40,8h!a!(3,2),
     $           t69,8hfunction/)
            a(1)=t11(ltheta,iphi)
            a(2)=t22(ltheta,iphi)
            a(3)=t12(ltheta,iphi)
            a(4)=t21(ltheta,iphi)
            do 30 k=1,4
               m(k)=a(k)*dconjg(a(k))
               do 40 j=k,4
                  s(k,j)=.50d0*(a(j)*dconjg(a(k))+a(k)*dconjg(a(j)))
                  s(j,k)=s(k,j)
                  d(k,j)=-.50d0*i*(a(j)*dconjg(a(k))-
     $                 a(k)*dconjg(a(j)))
                  d(j,k)=-d(k,j)
 40            continue
 30         continue
            f11=.50d0*(m(2)+m(3)+m(4)+m(1))
            f(1,1)=1.0d0
            f(1,2)=.50d0*(m(2)-m(3)+m(4)-m(1))/f11
            f(1,3)=(s(2,3)+s(4,1))/f11
            f(1,4)=(-d(2,3)-d(4,1))/f11
            f(2,1)=.50d0*(m(2)+m(3)-m(4)-m(1))/f11
            f(2,2)=.50d0*(m(2)-m(3)-m(4)+m(1))/f11
            f(2,3)=(s(2,3)-s(4,1))/f11
            f(2,4)=(-d(2,3)+d(4,1))/f11
            f(3,1)=(s(2,4)+s(3,1))/f11
            f(3,2)=(s(2,4)-s(3,1))/f11
            f(3,3)=(s(2,1)+s(3,4))/f11
            f(3,4)=(-d(2,1)+d(3,4))/f11
            f(4,1)=(d(2,4)+d(3,1))/f11
            f(4,2)=(d(2,4)-d(3,1))/f11
            f(4,3)=(d(2,1)+d(3,4))/f11
            f(4,4)=(s(2,1)-s(3,4))/f11
            jphi=(iphi-1)*incphi
            jtheta=(ltheta-1)*inceta
            if(ltheta.le.lhalf) goto 50
            jtheta=(lstop-ltheta+1)*inceta
            jphi=(iplane+iphi-1)*incphi
 50         theta=rad*dfloat(jtheta)
            phi=rad*dfloat(jphi)
            scat=dacos(dcos(anglin)*dcos(theta)+dsin(anglin)*
     $           dsin(theta)*dcos(phi))/rad
            write(6,100) jphi,jtheta,scat,m(1),m(4),
     $           f11,(f(1,j),j=1,4)
            write(6,101) m(3),m(2),(f(2,j),j=1,4)
            write(6,102) a(1),a(4),(f(3,j),j=1,4)
            write(6,103) a(3),a(2),(f(4,j),j=1,4)
 100        format(2x,i3,3x,i3,5x,f6.2,t29,2h! ,d10.4,
     $           2x,d10.4,2h !,t68,d10.4,t80,1h!,4d12.4,3h !i)
 101        format(t29,2h! ,d10.4,2x,d10.4,2h !,
     $           t80,1h!,4d12.4,3h !q)
 102        format(t18,1h!,4d12.4,2h !,t80,1h!,4d12.4,3h !u)
 103        format(t18,1h!,4d12.4,2h !,t80,1h!,4d12.4,3h !v/)
            lcount=lcount+1
            if(lcount.lt.12) goto 20
            lcount=0
 20      continue
 10   continue
      return
      end
c 
      subroutine lwl
      implicit real*8(a-h,o-z)
c extinction and polarization from homogeneous spheroids in the long
c wavelength limit
c
c derived from conf.f which was for
c extinction and polarization from confocal spheroids
c see gilra 1972 in results from oao2 ed. code 
c draine and lee 1985 apj 285,89
c martin 1974 apj 187, 461
c martin 1975 apj 202, 390(p?)
c cosmic dust, table 2.1 and sections on long-wave limit p.20 and p.65
c
      logical verbose
c *** input****
      complex*16 cfindx
      common/param/cfindx,dum3,anglin,dum4,idum
      complex*16 cdumb
      common/cincou/ cdumb,dum1,sgn,dum2
      common/lwlim/ a,b,rl
c local to lwl routines
      complex*16 epc,cpar,cper
      common /cls/ rlparc,rlperc
      common /ev/ epc
      common /volume/ vc
      common /wave/ rkay
      common /basic/ eppar,epper,sigpar,sigper
      common /spinf/ c2psi
      common /derivd/ delta,dsig,dep,vop,dsde,pot
      common /derivq/ qepar,qeper,qapar,qaper,qspar,qsper,
     $     qeppar,qepper,deltaq,dsigq,depq,dsdeq,albedo
c
c ac is the symmetry axis, bc is the equal axis
c always enter from spheroid program with a.gt.b
c
      verbose=.false.
c
      if(sgn.eq.1.d0) then
c prolate
         ac=a
         bc=b
      else
         ac=b
         bc=a
      endif
c ****** for thin disk, put ac .eq. 0. 
c ****** for long cylinder, put ac .lt. 0. 
c
c     rl=wavelength in microns
      pi=3.141592653589793d0
      rkay=1.d+4*2.d0*pi/rl 
      ak=ac*2.d0*pi/rl
      bk=bc*2.d0*pi/rl
c in cm-1
c input complex refractive index
      rnc=dreal(cfindx)
      rkc=dimag(cfindx)
c the complex refractive index is used as rnc - i*rkc (opposite convention)
c convert to dielectric function
      call emtoep(epc,rnc,rkc)
      if(verbose) write(6,5) rl,rnc,rkc,epc
 5    format(5(1x,1pe13.6))
c
c oblique orientation: 90 is symmetry axis in plane of sky
      c2psi = dcos(anglin)
      c2psi = c2psi*c2psi
c
c obliquity etc.
      if(ac.ge.0.d0 .and. ac.le.bc) then
c thin disk, oblate or sphere
c note: spinning oblate = static one, so field inclination = axis inclination
         if(verbose) write(6,777) anglin*180.d0/pi
 777     format(1x,' oblate with field at (or inclined by)'
     $        ,f7.3,' degrees')
      else
c prolate options
         if(verbose) write(6,779) anglin*180.d0/pi
 779     format(1x,' prolate static with inclination',f6.2,' degrees')
      endif
c
      vc = 4.d0*pi*ac*bc*bc/3.d0
      if(ac.le.0.d0) vc = 2.d0*pi*bc*bc/3.d0 
c volume/2a = vol/length
c volume per unit length for thin disk or for long cyl
c this is proportional to the area perp to symmetry axis
      if(verbose) write(6,*) ac,bc,vc
c
c calculate depolarization factors
      call valol(ac,bc,rlparc,rlperc)
      if(verbose) write(6,*) ' depolarization factors'
      if(verbose) write(6,5) ac,bc,rlparc,rlperc
c
c calculate complex cross section cext + i*cphase
      call ralyc(ac,cpar,cper,spar,sper)
      if(verbose) write(6,*) ' cpar,cper ', cpar,cper
      if(ac.le.0.) write(6,*) 
     $' ***scat cross-section is per unit vol **per unit length 2a**'
      if(verbose) write(6,*) ' spar,sper ', spar,sper
 667  continue
c
c calculate cross sections relating to
c phase shift and extinction per unit volume
      call ralyes(ac,cpar,cper,spar,sper,apar,aper)
c
c allow for spinning and calculate polarization
      call spinc(ac,bc)
      if(ac.le.0.d0) 
     1     write(6,*) ' extinction does not include scattering'
      if(verbose) write(6,300)
 300  format(' sigpar,sigper,eppar,epper, ',
     1     'delta,dsig,dep,dsde, ',
     2     'vop,pot')
      if(verbose) write(6,3) sigpar,sigper,eppar,epper
      if(verbose) write(6,3) delta,dsig,dep,dsde
      if(verbose) write(6,3) vop,pot
 3    format(5(1x,1pd13.6))
c
      call spinq(spar,sper,apar,aper,ac,bc)
      if(ac.lt.0.d0) 
     1  write(6,*) ' ***to get lind-greenberg normalization of q(a,ep),
     2     multiply by (pi/2)/(4/3)***'
      if(ac.le.0.d0) 
     1     write(6,*) ' ***scattering q is **per unit length 2a**'
      if(verbose) write(6,301)
 301  format(' qapar,qaper,qspar,qsper, ',
     1     'qepar,qeper,qeppar,qepper, ',
     2     'deltaq,dsigq,depq,dsdeq, ',
     3     'alb')
      if(verbose) write(6,3) qapar,qaper,qspar,qsper
      if(verbose) write(6,3) qepar,qeper,qeppar,qepper
      if(verbose) write(6,3) deltaq,dsigq,depq,dsdeq
      if(verbose) write(6,3) albedo
      if(verbose) write(6,*)
     $     ' tm mode; par: e parallel to symmetry axis'
      if(verbose) write(6,*) 
     $     ' te mode; per: e perpendicular to symmetry axis'
c
      akbk=dmax1(ak,bk)
      if(akbk.gt.1.d-2) 
     $     write(6,507) ak,bk
 507  format(1x,'lw limit warning: size parameters large',
     $     2(1x,1pd10.3))
c effective qepar for tilted orientation
      qepare=qepar*(1.d0-c2psi)+qeper*c2psi
      qppare=qeppar*(1.d0-c2psi)+qepper*c2psi
      qspare=qspar*(1.d0-c2psi)+qsper*c2psi
      write(6,501) qeper,qepare,deltaq,dsigq,
     $     qepper,qppare,depq,qsper,qspare
 501  format(' lw limit side-on normal',t25,9(1x,1pd10.3))
      write(6,506) pot,vop,albedo
 506  format(1x,3(1x,d11.4))
c     
      return
      end
c
      subroutine valol(sa,sb,srlpar,srlper)
c gives values of l the depolarization factors for two modes for spheroids
c a/b gt 1 means prolate
      implicit real*8(a-h,o-z)
      a = sa
      b = sb
      if(dabs((a-b)/b).lt.1.d-06) go to 5 
c sphere
      if(a.le.0.d0) go to 6
      c=b*b/(a*a)
      e2=dabs(1.d0-c)
      e=dsqrt(e2)
c a is symmetry axis; parallel
      if(a.gt.b) 
     1     rlpar=(1.d0-e2)*(-1.d0+dlog((1.d0+e)/(1.d0-e))/(2.d0*e))/e2 
c prolate
      if(a.lt.b) 
     1     rlpar=(1.d0+e2)*(1.d0-datan(e)/e)/e2 
c oblate
      rlper=(1.d0-rlpar)/2.d0
      go to 999
c sphere
 5    rlpar=1.d0/3.d0
      rlper=rlpar
      go to 999
 6    continue 
c thin disk or long cylinder
      rlpar = 0.d0
      rlper = 0.d0
      if(a.eq.0.d0) rlpar = 1.d0 
c thin disk
      if(a.lt.0.d0) rlper = 0.5d0 
c long cylinder
 999  continue
      srlpar = rlpar
      srlper = rlper
      return
      end
c
      subroutine emtoep(ep,rn,rk)
      implicit real*8(a-h,o-z)
c complex refractive index rn - i*rk 
c as in gilra 1972
c change to complex dielectric constant
      complex*16 ci,cm,ep
      ci=(0.0d0,1.d0)
      cm=rn-rk*ci
      ep = cm*cm
      return
      end
c
      subroutine ralyc(ac,cpar,cper,spar,sper)
      implicit real*8(a-h,o-z)
c complex cross-sections
      complex*16 cpar,cper
      common /cls/ rlparc,rlperc
c
c e parallel to symmetry axis
      call cov(ac,cpar,spar,rlparc)
c e pependicular to symmetry axis
      call cov(ac,cper,sper,rlperc)
      return
      end
c
      subroutine cov(ac,cv,sv,rlc)
      implicit real*8(a-h,o-z)
c compute complex cross section (from i*alpha over v(mantle))
      complex*16 epc,ci,cv,av,dconjg 
c diel consts, cross section, 4pi alpha/v
      complex*16 b,c
      common /ev/ epc
      common /volume/ vc
      common /wave/ rkay
      pi=3.141592653589793d0
      b = epc-1.d0
      c = rlc*(epc-1.d0)+1.d0
      av = b/c
c as in gilra 1972, without 4pi
      ci = (0.0d0,1.d0)
      cv = ci*rkay*av 
c as in gilra 1972, 4pi cancels
      sv = dreal( av*dconjg(av) )
      rr = rkay*rkay 
c already in cm-1
      spi = 6.d0*pi
      if(ac.gt.0.d0) sv = sv*rr*rr*(vc*1.d-12)/spi 
c in cm-1
      if(ac.lt.0.d0) sv = sv*rr*rr*(vc*1.d-8)/spi 
c in cm-2
      if(ac.eq.0.d0) sv = sv*rr*rr*(vc*1.d-8)/spi 
c in cm-2
c note that for long cylinders and thin disks, vc is the area 2/3 pi*b*b,
c (rather than pi*b*b), since this*length (2a) gives the volume
c the volume is area*l where l is the total length, so this gives
c the scattering cross-section per unit volume **** per unit length ***
c so do not add it directly to absorption which is per unit volume
c
c here the basic difference of long wrt infinite cylinders becomes important: 
c for long cylinders, still assume that wavelength .gt. length
c scattering cross-section varies as vol*2/wavelength**4,
c or (lb**2)**2/wavelength**4
c or as (l/b)*(l*2b)*(b/wavelength)**4
c scattering cross-section per vol per unit length varies as b*b/wavelength**4
c this is calculated originally, and is ok even when l is undefined
c
c scattering cross-section per unit length varies as (l/b)*b*(b/wavelength)**4
c scattering efficiency (q=cu/2b) varies as (l/b)*(b/wavelength)**4
c i.e., like rayleigh scattering dependences when l is proportional to b
c or alternatively q varies as (l/wavelength)*(b/wavelength)**3
c when l is undefined, can only calculate scattering efficiency per unit length
c q/l varies as b**3/wavelength**4
c
c
c for infinite cylinders, wavelength .lt. l and l is essentially undefined
c scattering efficiency (q=cu/2b) varies as (b/wavelength)**3
c this is calculated originally
c
c scattering cross-section varies as l*b*(b/wavelength)**3
c scattering cross-section per unit length varies as b*(b/wavelength)**3
c scattering cross-section per vol per unit length varies as 
c b*(b/wavelength)**3/lb*b = b**2/wavelength**4 * (wavelength/l)
      return
      end
c
      subroutine ralyes(ac,cpar,cper,spar,sper,apar,aper)
      implicit real*8(a-h,o-z)
      complex*16 cpar,cper
      common /basic/ eppar,epper,sigpar,sigper
c use c = cabs + i*cphase + cscatt
c and sig = cext/2 and eps = cphase/2
c par: e parallel to symmetry axis
c per: e perpendicular to symmetry axis
c
c absorption per unit volume
      apar = real(cpar) 
      sigpar=0.5d0*apar 
      aper = dreal(cper)
      sigper=0.5d0*aper
      eppar=0.5d0*dimag(cpar)
      epper=0.5d0*dimag(cper)
      if(ac.gt.0.d0) then
c note that for long cylinders and thin disks, the scattering
c cross-section (spar or sper) is per unit volume **** per unit length ***
c so do not add it directly to absorption which is per unit volume, i.e. 
c only add scattering contribution if not a long cylinder or a thin disk
         sigpar = sigpar+0.5d0*spar
         sigper = sigper+0.5d0*sper 
      endif
      return
      end
c
      subroutine spinc(ac,bc)
      implicit real*8(a-h,o-z)
      common /basic/ eppar,epper,sigpar,sigper
      common /spinf/ c2psi
      common /derivd/ delta,dsig,dep,vop,dsde,pot
c
c now allow for oblique orientation in forming cross-sections per unit volume
c
c phi is angle between line of sight and symmetry axis
c for spheroids for either absorption or scattering 
c (not for infinite cylinders)
c (see cosmic dust, equations 2.34 and 4.17)
c this is for arbitrary oblique cylinder
c       sigpar-tipped = sigpar*s2psi+sigper*c2psi
c	sigper-tipped = sigper
c
      delta=sigpar*(1.d0-c2psi)+sigper*(1.d0+c2psi) 
      dsig = (sigpar-sigper)*(1.d0-c2psi) 
      dep = (eppar-epper)*(1.d0-c2psi) 
      if(ac.ge.0.d0 .and. ac.le.bc) then
c oblate (allows for inf cyl; ac<0)
c         dsig = -dsig
c         dep = -dep
c don't change sign, so this now is as in main spheroid program
         dsig = dsig
         dep = dep
      endif
      dsde = dsig*dep
      vop = -1.d0
      if(dsig.ne.0.d0)	vop = dep/dsig
      pot = -1.d0
      if(delta.ne.0.d0)	pot =dsig/delta
      return
      end
c
      subroutine spinq(spar,sper,apar,aper,ac,bc)
      implicit real*8(a-h,o-z)
      common /basic/ eppar,epper,sigpar,sigper
      common /spinf/ c2psi
      common /derivq/ qepar,qeper,qapar,qaper,qspar,qsper,
     1     qeppar,qepper,deltaq,dsigq,depq,dsdeq,albedo
c
c no modification to area for c2psi.ne.0 : side-on area used
      fact = (4.0d0/3.d0)*bc*1.d-4
c
c pi=3.141592653589793d0
c if(ac.le.0.) fact = (3.1415927/2.)*bc*1.e-4 
c
c use this if q (abs and phase) for long cylinders (and thin disks) is to 
c agree with lind and greenberg normalization for infinite cylinders
      qapar = fact*apar
      qaper = fact*aper
      qeppar = 2.d0*fact*eppar
      qepper = 2.d0*fact*epper
c
c scattering contribution
c	fact = (4./3.)*bc*1.e-4 
c restore lind and greenberg corr for continuity with finite sized spheroids
c of length 2*ac
      qspar = fact*spar
      qsper = fact*sper
c
      qepar = qapar 
      qeper = qaper
      if(ac.gt.0.d0) then
c this is only approximate, since terms of same order as 
c qs are being ignored in qa
         qepar = qepar + qspar
         qeper = qeper + qsper
      endif
c note that for long cylinders and thin disks, the scattering
c efficiency factor (qspar or qsper) is **** per unit length ***
c so do not add it directly to absorption which is dimensionless, i.e. 
c only add scattering contribution if not a long cylinder or a thin disk
c
      deltaq = (qepar*(1.d0-c2psi)+qeper*(1.d0+c2psi))/2.d0
      deltas = (qspar*(1.d0-c2psi)+qsper*(1.d0+c2psi))/2.d0 
      dsigq = (qepar-qeper)*(1.d0-c2psi)/2.d0 
      depq = (qeppar-qepper)*(1.d0-c2psi)/2.d0
      albedo = -1.d0
c just a signal that irrelevant for thin disks and inf cyl
      if(ac.gt.0.d0) albedo = deltas/deltaq
      dsdeq = dsigq*depq
      return
      end
c
