c     corrstarter.for   version 1.21 Jan. 4, 1994.
      program corrstarter
      implicit none
      integer itype,iq,npart,l1,l2,nkmax,kk,ii,jj,ll
      double precision mom,delmom,mass1,mass2,akk
      double precision wsym1,wsym2,wanti1,wanti2,root2
      double precision rmass,phase1,phase2,eta,qz,pi
      double complex phicon(1000)
      double complex cd1(1000), cd2(1000)
      common eta,rmass,npart,itype
c     ------------------------------------
c     When you vary the particle, e.g. pp, pi-pi, nn ,KK to do using
c     interferometry, you should change these parameters.
c     rmass is the reduced mass.
c     qz is the charge. = 1 for like charge, -1 for opposite..., 0
c     for Z1*Z2 = 0
c     Note: also go into vr and adjust the potential for appropriate
c     pair.
      root2=dsqrt(2.d0)
      pi=dacos(-1.d0)
      write(6,*) 'enter an integer...'
      write(6,*) ' 1 : for (pi-,pi-) or (pi+,pi+)'
      write(6,*) ' 2 : for (pi0,pi0)'
      write(6,*) ' 3 : for (K-,K-) or (K+,K+)'
      write(6,*) ' 4 : for (p,p)'
      write(6,*) ' 5 : for (n,n)'
      write(6,*) ' 6 : for (p,n)'
      write(6,*) ' 7 : for (d,d)'
      write(6,*) ' 8 : for (K+,K-)'
      write(6,*) ' 9 : for (K-,pi+)'
      write(6,*) '10 : for (K+,p)'
      write(6,*) '11 : for (d,alpha)'
      write(6,*) '12 : for (p,pi+)'
      write(6,*) '13 : for (p,alpha)'
      read(5,*) itype
c     First set masses
c     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      if(itype.eq.1) then
         mass1=139.58d0
         mass2=139.58d0
      endif
      if(itype.eq.2) then
         mass1=134.98d0
         mass2=134.98d0
      endif
      if(itype.eq.3.or.itype.eq.8) then
         mass1=493.646d0
         mass2=493.646d0
      endif
      if(itype.eq.4) then
         mass1=938.28d0
         mass2=938.28d0
      endif
      if(itype.eq.5) then
         mass1=939.56d0
         mass2=939.56d0
      endif
      if(itype.eq.6) then
         mass1=938.28d0
         mass2=939.56d0
      endif
      if(itype.eq.7) then
         mass1=1875.65d0
         mass2=1875.65d0
      endif	
      if(itype.eq.9) then
         mass1=493.646d0
         mass2=139.58d0
      endif
      if(itype.eq.10) then
         mass1=493.646d0
         mass2=938.28d0
      endif
      if(itype.eq.11) then
         mass1=1875.585d0
         mass2=3727.323d0
      endif
      if(itype.eq.12) then
         mass1=938.28d0
         mass2=139.58d0
      endif
      if(itype.eq.13) then
         mass1=938.28d0
         mass2=3726.35d0
      endif
      rmass=(mass1*mass2)/(mass1+mass2)
c     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c     Now set symmetrization weights and no.s of partial waves to calc.
c     Strong int. corrections for
c     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      wsym2=0.d0
      wanti2=0.d0
      l1=0
      l2=0
      if(itype.eq.1) then
         wsym1=1.d0
c     A non-relativistic potential is used here, it is not a
c     very good treatment, but the effect is negligible for HI sources.
c     With some work, one could do better, which might be important
c     for very small sources.
         npart=1
c         npart=0
         wanti1=0.d0
         l1=0
      endif
      if(itype.eq.2) then
         wsym1=1.d0
         npart=0
         wanti1=0.d0
         l1=0
      endif

      if(itype.eq.3) then
         wsym1=1.d0
         wanti1=0.d0
         npart=0
      endif
      if(itype.eq.4.or.itype.eq.5) then
         wsym1=.25d0
         wanti2=.75d0
         npart=1
         l1=0
         l2=1
      endif
      if(itype.eq.6) then
         wsym2=.125d0
         wsym1=.375d0
         wanti1=1.d0-wsym1-wsym2
         l1=0
         l2=0
         npart=2
      endif
      if(itype.eq.7) then
         wsym1=.55555555555555555d0
         wsym2=.11111111111111111d0
         wanti1=1.-wsym1-wsym2
         wanti2=0.d0
         npart=2
         l1=0
         l2=0
      endif
      if(itype.eq.8) then
         wsym1=.5d0
         wsym2=0.d0
         wanti1=0.d0
         wanti2=.5d0
         iq=-1
         npart=0
      endif
      if(itype.eq.9.or.itype.eq.10) then
         wsym1=.5d0
         wsym2=0.d0
         wanti1=0.d0
         wanti2=.5d0
         npart=0
      endif
      if(itype.eq.11) then
         wsym1=7.d0/30.d0
         wanti1=0.d0
         wsym2=8.d0/30.d0
         wanti2=.5d0
         npart=0
c     npart=1
         l1=2
         iq=2
      endif
      if(itype.eq.12) then
         wsym1=.5d0
         wsym2=0.d0
         wanti1=0.d0
         wanti2=.5d0
         npart=0
      endif
      if(itype.eq.13) then
         wsym1=.5d0
         wsym2=0.d0
         wanti1=0.d0
         wanti2=.5d0
         npart=0
         iq=2
      endif	
c     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c     Now set charges
c     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      iq=0
      if(itype.eq.1.or.itype.eq.3.or.itype.eq.4.or.itype.eq.7.or.
     +     itype.eq.9.or.itype.eq.10.or.itype.eq.12) iq=1
      if(itype.eq.8) iq=-1
      if(itype.eq.11.or.itype.eq.13) iq=2
      qz=dfloat(iq)
c     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c     -------------------------------------
c     mom is the momentum, mom is double precision mom, eta=usual eta
c     akk is the upper boundary of the kth bin, while mom is 
c     the actual momentum used for the calculations of the
c     wavefunctions.
c     cds and cdp are the change in the relative wave functions for
c     s and p waves due to the relative potential.
c     cds and cdp are stored as a function of x=kr for delx=.d01 
c     in an array
c     of 1000, (x>0,x<10) The program diffy calculates these values as
c     a function of the relative momentum, mom(MeV/c).
c     ------------------
c     nkmax is the number of momentum values

      write(6,*) 'what is nkmax?'
      read(5,*) nkmax
      mom=0.d0
      write(6,*) 'What is the size of the mometum bin? (in MeV/c)'
      write(6,*) 'This program works with reduced rel. momentum, so put'
      write(6,*) 'in 1/2 of (p1-p2) if (p1-p2) is what you want:'
      read(5,*) delmom
c     This loop calculates the correction to the rel. wave func.
c     for various
c     momenta.  jj refers to kr and kk refers to the momentum.
      write(25) mass1,mass2
      write(25) wsym1,wsym2,wanti1,wanti2
      write(25) iq,itype
      write(25) npart,l1,l2
      write(25) nkmax
      do 990 kk=1,nkmax
c     if(kk.gt.10) delmom=5.d0
c     if(kk.gt.15) delmom=10.d0
         mom=mom+delmom
         akk=mom+delmom/2.d0
         eta=qz*(rmass/137.036d0)/mom
         if(npart.gt.0) call diffy(mom,cd1,cd2,phase1,phase2,l1,l2)
         if(phase1.lt.0.d0) phase1=phase1+2.d0*3.14159265358979323844d0
         if(phase2.lt.0.d0) phase2=phase2+2.d0*3.14159265358979323844d0
         write(6,165) mom,phase1,phase2
 165     format(' mom=',f8.2,3x,'phase shifts :',f9.4,5x,f9.4)
         write(25) mom,akk
         if(npart.gt.0) write(25) (cd1(ii),ii=1,1000)
         if(npart.gt.1) write(25) (cd2(jj),jj=1,1000)
c     #########################################
c     if(itype.eq.8) call phistarter(mom,phicon)
c     if(itype.eq.8) write(25) (phicon(jj),jj=1,1000)
c     if(itype.eq.8) write(56,*) '---',kk,(phicon(jj),jj=1,40)
c     #########################################
         if(npart.gt.0) write(25) phase1
         if(npart.gt.1) write(25) phase2
 990  continue
      stop
      end
**********************************
      double complex function cw2(r,ll)
      implicit none
      double complex z,a,b,f1,top1,top2,bot,delf,i
      double precision r,eta,l,rmass,arg,pi
      integer npart,itype,n,j,ll
      common eta,rmass,npart,itype
      pi=dacos(-1.d0)
      i=(0.d0,1.d0)
      l=dfloat(ll)
      z=-2*r*i
      a=l+1+i*eta
      b=2.*(l+1.d0)
      f1=(1.d0,0.d0)
      top1=(1.d0,0.d0)
      top2=(1.d0,0.d0)
      bot=(1.d0,0.d0)
      n=5
      do 75 j=1,n
         top1=top1*(dfloat(j)-a)
         top2=top2*(dfloat(j)+b-a-1.d0)
         bot=bot*dfloat(j)
         delf=(top1*top2)/(bot*((z)**j))
         f1=f1+delf
 75   continue
      arg=eta*dlog(2.d0*r)
      arg=arg-2.d0*pi*int(arg/(2.d0*pi))
      cw2=cdexp(i*arg)*f1
      return
      end
*************************************
      double precision function vr(X,jj,lpass,p)
      implicit none
c     This is the Reid potential stolen out of Koonin's code
c     (Capital letters have not been edited).  I did check 
c     the phase shifts
c     to compare to exp. to see that it looks o.k.
      double precision eta,rmass,p,pmux,con,f1,f2,f3,f4,f6,f7,vrd
      double precision x,fudge,fudge4,fudge6,r
      integer lpass,jj,npart,itype
      common eta,rmass,npart,itype
      CON=2.*rmass/(p*p)
      if(jj.eq.1.or.(X*197.323d0/p).gt.14.d0) then
         vr=0.d0
         return
      endif
c     --------------------
      if(itype.eq.4.or.itype.eq.5) then
         if(lpass.eq.1) then
            PMUX=X*197.323d0*.7d0/p
            f1=dexp(-PMUX)
            f4=f1**4
            f7=f4*(f1**3)
            VRD=-10.463d0*f1/PMUX-1650.6d0*f4/PMUX
            VRD=VRD+6484.2*f7/PMUX
         endif
         if(lpass.eq.2) then
            PMUX=X*197.323d0*.7d0/p
            f1=dexp(-PMUX)
            f2=f1**2
            f4=f2**2
            f3=f2*f1
            VRD=(1.d0+2.d0/PMUX+2.d0/(PMUX**2))*f1
            VRD=VRD-(8.d0/PMUX+2./(PMUX**2))*f4
            VRD=10.463d0*VRD/PMUX
            VRD=VRD-135.25d0*f2/PMUX+472.81d0*f3/PMUX
         endif
      endif
      if(itype.eq.1.or.itype.eq.2) then
         PMUX=X*760.d0/p
         VRD=2500.d0*dexp(-PMUX)/PMUX 
      endif
      if(itype.eq.6) then
         f1=dexp(-PMUX)
         if(lpass.eq.2) then
            PMUX=X*197.323*.7d0/p
            f1=dexp(-PMUX)
            f4=f1**4
            f7=f4*(f1**3)
            VRD=-10.463d0*f1/PMUX-1650.6d0*f4/PMUX
            VRD=VRD+6484.2d0*f7/PMUX
         endif
         if(lpass.eq.1) then
            fudge4=1.4579d0
            fudge6=1.15d0
            PMUX=X*197.323d0*.7d0/p
            f1=dexp(-PMUX)
            f2=f1**2
            f4=f2**2
            f6=f2*f4
            VRD=-10.463d0*f1+105.468d0*f2-
     +           fudge4*3187.8d0*f4+fudge6*9924.3d0*f6
            VRD=VRD/PMUX
         endif                  
      endif
      if(itype.eq.7) then
         r=X*197.323d0/p
         if(lpass.eq.1) VRD=29.8d0/(1.+dexp((r-4.21d0)/.134d0))
         if(lpass.eq.2) VRD=26.d0/(1.+dexp((r-1.08d0)/1.25d0))
      endif
      if(itype.eq.11) then
         r=X*197.323d0/p
         VRD=0.d0
         if(lpass.eq.1) VRD=-42.045d0/(1.+dexp((r-2.7648d0)/.7d0))
      endif
      vr=VRD*CON
      return
      end
*********************************************
      subroutine diffy(mom,delf0,delf1,phase11,phase22,l1,l2)
      implicit none
      double precision vr,delr,r00,rmax,mom,l,phase11,phase22
      double precision vreid,r,r1,r2,r0,vcs,rvcs
      double precision b,eta,rmass
      double precision phase0,phase1,phase2,oldphase
      double complex  cf(1000), cfin(2001),cf1(1000),cf2(1000)
      double complex cy0,cy1,cy2,cw2,ceip,cww2
      double complex delf0(1000),delf1(1000)
      double complex i,phasec
      integer mm,lll,l2,l1,jj,nmax,jjj,lpass,j,itype,ll,npart
      common eta,rmass,npart,itype
      i=(0.d0,1.d0)
c     Wave functions are found as a function of kr, referred to as
c     r here.  Values are calculated for delr=.d0005*krmax,
c     krmax=40fm*mom
c     other value is passed on to the main prog.  
c     So the main prog. uses
c     a mesh of delr=.d001*krmax
      rmax=40.*mom/197.323
      delr=-.0005d0*rmax
      r00=rmax
c     -----------------
c     if(r00.lt.6.and.eta.gt.1E-5) then
c     r00=-delr*int(-6./delr)
c     endif
c     -----------------
      nmax=nint(r00/(-delr))
      if(nmax.gt.30000) write(6,*) 'out of bounds, nmax too big=',nmax
********************ANGULAR MOMENTUM LOOP**********************
c     We do calc. for l=0 and l=1,  ll=l=angular momentum
      do 619 lll=1,npart
         lpass=lll
         if(lpass.eq.1) ll=l1
         if(lpass.eq.2) ll=l2
         l=dfloat(ll)
******************W & WO REID POTENTIAL LOOP*******************
c     We perform calc with and with out the Reid pot.
c     We are only interested in the correction to the wave function
c     due to interactions. The first
c     time through we turn off the Reid pot.(jj=1,off) (jj=2,on)
         do 577 jjj=1,2
            jj=jjj
            r0=r00
            r1=r0+delr
c     This initiates an incoming coul. wave function at
c     two points at r=10.  See Messiah's appendix on incoming partial
c     Coul. waves for explanation.  We then calc. inward until r=0.
            if(r0.gt.10.d0.or.eta.lt.0.0001d0) then
	       cy0=cw2(r0,ll)*cdexp(-i*r0)
	       cy1=cw2(r1,ll)*cdexp(-i*r1)
               if(ll.eq.1) cy0=cy0*i
               if(ll.eq.1) cy1=cy1*i
               if(ll.eq.2) cy0=-cy0
               if(ll.eq.2) cy1=-cy1
            else
               cy0=cww2(ll,r0,eta)
               cy1=cww2(ll,r1,eta)
               if(ll.eq.1) cy0=cy0*i
               if(ll.eq.1) cy1=cy1*i
               if(ll.eq.2) cy0=-cy0
               if(ll.eq.2) cy1=-cy1
            endif
            cfin(nmax)=cy0
            cfin(nmax-1)=cy1
c     Given the wave func. at two pts the wave func. can be found at
c     third pt. from diff. eq.
            do 111 j=1,nmax-2
               r2=r1+delr
               vreid=vr(r1,jj,lpass,mom)
               rvcs=(r1/mom)*197.323d0
               vcs=1.d0
               if(itype.eq.11.and.rvcs.le.4.d0) vcs=(rvcs/4.d0)**3
               b=1.d0-((l*(l+1))/(r1**2))-(2*eta*vcs/r1)-vreid
               cy2=(-b*cy1*(delr**2)+2*cy1-cy0)
               cfin(nmax-1-j) = cy2
               r1=r2
               cy0=cy1
               cy1=cy2
 111        continue
c     phase0 is the phase of the wavefunction at r=0.  
c     Subtracting the complex
c     conjugate of the wavefunction *exp(-2i*phase0) 
c     from the wavefunction
c     yields a sol. which satisfies the bound. conditions.Multiply by i
c     to correspond to the correct incoming partia-wave phase(
c     See Messiah).
c     finally, when coming through after doing it for both with and 
c     withoutinteractions, take the difference.
            phase1=-i*cdlog(cfin(1)/cdabs(cfin(1)))
            phase2=-i*cdlog(cfin(2)/cdabs(cfin(2)))
            phase0=2*phase1-phase2
            ceip=cdexp(i*phase0)
            j=0
c     Note that we send back to the main prog. the answer with a less
c     dense mesh(delr=.d01)
            do 112 mm=1,1999,2
               j=j+1
               cf(j)=i*(cfin(mm)-dconjg(cfin(mm))*(ceip**2))/2
               if(jj.eq.1) cf1(j)=cf(j)
               if(jj.eq.2) cf2(j)=cf(j)
               if(jj.eq.2.and.lpass.eq.1) delf0(j)=cf2(j)-cf1(j)
               if(jj.eq.2.and.lpass.eq.2) delf1(j)=cf2(j)-cf1(j)
 112        continue
            if(jj.eq.2.and.lpass.eq.1) phase11=phase0-oldphase
            if(jj.eq.2.and.lpass.eq.2) phase22=phase0-oldphase
            oldphase=phase0
 577     continue
**************end of big loop****************************************
 619  continue
      return
      end

c     ***************************************************
      double complex function cww2(l,rho,deta)
c     The notation is like Gr. + R, page 1063.
c     The coulomb wave function is the same as W(i*deta,l+1/2,2*i*rho)
      implicit none
      double precision pi,euler
      parameter (pi=3.14159265358979323844,
     +     euler=0.5772156649015328606)
      double complex cdgamma,i,factor,lterm,fact1,fact2,
     +     psi1,psi2,psi3,cx,sum1,sum2,delsum1,delp,cdcon1,cdcon2
      parameter(i=(0.d0,1.d0))
      double precision deta,rho,dgamma
      integer l,ii,k,twol
      cdcon1=cdgamma(-dfloat(l)-i*deta)
      cdcon2=cdgamma(dfloat(l+1)-i*deta)
      factor=((-1)**(2*l+1))*((2*i*rho)**(l+1))*cdexp(-i*rho)/
     +     (cdcon1*cdcon2)
      psi1=-euler
      psi2=-euler
      do 456 ii=1,2*l+1
         psi2=psi2+1.d0/dfloat(ii)
 456  continue
      cx=-i*deta
      psi3=-euler
      do 457 ii=1,10000000
         delp=-1.d0/(cx+dfloat(ii))+1.d0/dfloat(ii)
         psi3=psi3+delp
         if(cdabs(delp).lt.1E-9) go to 987
 457  continue
      write(6,*) 'psi3 did not converge', cdabs(delp)
 987  continue
      do 458 ii=1,l
         psi3=psi3+1.d0/(cx+dfloat(ii))
 458  continue
c     write(6,*) 'the 3 psis', psi1,psi2,psi3
      lterm=cdlog(2*i*rho)
      fact1=cdcon2/dgamma(2*l+2)
      sum1=fact1*(psi1+psi2-psi3-lterm)
      do 123 k=1,1000000
         fact1=fact1*(2*i*rho)*(dfloat(l+k)-i*deta)/
     +        (dfloat(k)*dfloat(2*l+1+k))
         psi1=psi1+1.d0/dfloat(k)
         psi2=psi2+1.d0/dfloat(2*l+1+k)
         psi3=psi3+1.d0/(dfloat(l+k)-i*deta)
         delsum1=fact1*(psi1+psi2-psi3-lterm)
         sum1=sum1+delsum1
         if(cdabs(delsum1).lt.1E-15) go to 888
 123  continue
      write(6,*) 'never escaped the loop!'
 888  fact2=dgamma(2*l+1)*cdcon1/((-2*i*rho)**(2*l+1))
      sum2=fact2
      twol=2*l
      do 124 k=1,twol
         fact2=fact2*(dfloat(k-l-1)-i*deta)*
     +        (-2*i*rho)/(dfloat(k)*dfloat(2*l-k+1))
         sum2=sum2+fact2
 124  continue
      sum1=factor*sum1
      sum2=factor*sum2
      cww2=(sum1+sum2)*dexp(pi*deta/2.d0)
c     write(6,*) 'cw2,sum1,sum2',cw2,sum1,sum2
      return
      end
*****************************************
      double complex function cdgamma(c)
c     This calc.s gamma functions which are in the form gamma(n+i*y)
c     where n is an integer and y is real.
      implicit none
      double complex c,ci,cg,cphase
      integer mm,j
      parameter (ci=(0.d0,1.d0))
      double precision x,y,phase,euler,delp,cgmag,pi
      parameter (pi=3.14159265358979323844,euler=0.5772156649011673)
      x=c
      y=-ci*c
      phase=-euler*y
      do 123 j=1,50000
         delp=(y/dfloat(j))-datan(y/dfloat(j))
         phase=phase+delp
         if(dabs(delp).lt.1E-9) go to 223
 123  continue
      write(6,*) 'oops not accurate enough, increase jmax'
 223  phase=phase-2.d0*pi*int(phase/(2.d0*pi))
      cphase=cdexp(ci*phase)
      cgmag=dsqrt(pi*y/dsinh(pi*y))
      mm=nint(x)
      cg=cgmag*cphase
      if(mm.lt.1) then
         do 444 j=1,-mm+1
            cg=cg/(1.d0+dfloat(-j)+ci*y)
 444     continue
      endif
      if(mm.gt.1) then
         do 333 j=1,mm-1
            cg=cg*(dfloat(j)+ci*y)
 333     continue
      endif
      cdgamma=cg
      return
      end
******************************************

      double precision function dgamma(mm)
c     This calc.s gamma functions which are in the form gamma(n)
c     where n is an integer > 0.
      implicit none
      double precision cg
      integer mm,j
      cg=1.d0
      if(mm.lt.1) write(6,*) 'whoa boy, mm <1'
      if(mm.gt.1) then
         do 333 j=1,mm-1
            cg=cg*(dfloat(j))
 333     continue
      endif
      dgamma=cg
      return
      end
c     
      subroutine phistarter(p,phicon)
      implicit none
      complex ci,u
      complex phicon(1000)
      double precision p,g,pi,mphi,gamma,pc,mk,kzero,delx,x,qc
      double precision int,eq,q,delq,arg
      integer i,ii
      mk=493.646
      mphi=1019.41
      pi=3.141592654
      pc=dsqrt((mphi/2.)**2-mk**2)
      g=6.5
      ci=(0.d0,1.d0)
      kzero=2.d0*dsqrt(mk**2+p**2)
      gamma=(g**2)/(6*pi*kzero)
      gamma=gamma*((kzero/2)**2-mk**2)**1.5
      delx=40.d0/(197.323*1000.d0)
      x=-delx/2.d0
      do 727 i=1,1000
         x=x+delx
         phicon(i)=0.d0
         qc=dsqrt((kzero/2.)**2-mk**2)
         delq=qc/20.d0
         q=-delq/2.d0
         int=0.d0
         do 123 ii=1,250
            q=q+delq
            eq=dsqrt(mk**2+q**2)
            arg=cos(2*q*x)-sin(2*q*x)/(2*q*x)
            arg=arg/((kzero/2)**2-eq**2)
            arg=delq*arg*(q**2)/eq
            int=int+arg
 123     continue
         u=pi*pc*(cos(2*pc*x)-sin(2*pc*x)/(2*pc*x))
         u=u+ci*kzero*int
         u=u*(g**2)*p/((2*pi)**2)
         u=u/(kzero**2-mphi**2+ci*gamma)
         u=u/(kzero*x)
         phicon(i)=u
 727  continue
      return
      end