!-------------------------------------------------
      program swlparmix
!     
! Compute energy of two unlike but parallel square well lines
! dec-09
!
! Simple main routine, calls sr potential(...):
!
      implicit real*8(a-h,o-z)
      common pi,pih,rho,z,zm2h,a,a2,rlo
     1,rup,r,dr,rhmin,rhmax,drh,gam,gam1,gam2,gam3,gam4,del,delmax
     2,delmin,dgd,dz,sx12,sx23,sx34,sx32,sx24,sx,nrh,nr
      logical lovr,lint
!
      pi=4.d0*datan(1.d0) 
      pih=0.5d0*pi
!
      open(unit=7,file='swlparmix.out',form='formatted')
      open(unit=8,file='swlparmix.dat',form='formatted')
!
!   Parallel square well lines:
      write(7,100) 
100   format(1x,"SW Lines, parallel, mixed:")
!
!   Particle parameters:
      xl1=2.d0 
      xl2=2.d0 
      s2=1.5d0 ! (s1=1. by assumption)
!
!   Configuration: Particle 1 at origin
      x1=0.d0
      z1=0.d0
      x2=1.0001d0
      z2=0.d0 

      x12=x2-x1
      z12=z2-z1 
      call potential(xl1,xl2,s2,x12,z12,lovr,lint,ua)
!
      stop
      end
!-----------------------------------------------------------
!
      subroutine potential(xl1,xl2,s2,x12,z12,lovr,lint,ua)
!
! Does the actual calculation according to the recipe in ./swlparmix.html
!
      implicit real*8(a-h,o-z)
      common pi,pih,xl,h,rho,z,zm2h,a,a2,rlo
     1,rup,r,dr,rhmin,rhmax,drh,gam,gam1,gam2,gam3,gam4,del,delmax
     2,delmin,dgd,dz,sx12,sx23,sx34,sx32,sx24,sx,nrh,nr
      logical lovr,lint
!
      s2sq=s2*s2  
      h1=0.5d0*xl1 
      h2=0.5d0*xl2 
!
! If necessary, make z12 positive:
      if (z12.lt.0.d0) z12=-z12
! Calculate useful parameters:
      b=z12
      bsq=b*b
      asq=s2sq-x12*x12-z12*z12
      rtba=dsqrt(bsq+asq)
!
! Initialize potential:
      ua=0.d0  
! Exclude hard overlap:
      lovr=.true.
      xsq=x12*x12
      bra=z12-(h1+h2)
      if (xsq > 1.d0) then
        lovr=.false.
      else
        if ((z12.gt.(h1+h2)).and.((xsq+bra*bra).gt.1.d0)) lovr=.false.
      endif
! Check for interaction:
      lint=.true.
      if (xsq > s2sq) then
        lint=.false.
      else
        if ((z12.gt.(h1+h2)).and.((xsq+bra*bra).gt.s2sq)) then 
          lint=.false.
        endif
      endif
!
      if (.not.lovr.and.lint) then
!
! Compute points 1-4:
        xlam1=min(h1,max(-h1,b-h2-rtba))
        xlam2=min(h1,max(-h1,b-h2+rtba))
        xlam3=min(h1,max(-h1,b+h2-rtba))
        xlam4=min(h1,max(-h1,b+h2+rtba))  
!
! Compute shaded (=shared) area:
        sx0=2.d0*rtba*(xlam4-xlam1)
        sx12=0.5d0*(xlam2*xlam2-xlam1*xlam1)
        sx12=sx12-(b-h2+rtba)*(xlam2-xlam1)
        sx34=0.5d0*(xlam4*xlam4-xlam3*xlam3)
        sx34=sx34-(b+h2-rtba)*(xlam4-xlam3)
        sx=sx0+sx12-sx34
        ua=-sx/xl1/xl2
!
      write(6,104) sx0,rtba,xlam4,xlam1
104   format(1x,"sx0,rtba,xlam4,xlam1:",4f10.5)
      endif 
      write(6,103) x12,z12,ua,lovr,lint
      write(7,103) x12,z12,ua,lovr,lint
103   format(1x,"x12,z12,ua=:",3f10.5,2l5)
!
      return 
      end