!
! STATUS: Oct 28, 2010
!
!------------------------------------------------------------
! swl3d.f: computes potential btw 2 unlike swl particles
! The code is in generic fortran and should be readable 
! for any compiler
!
! The code closely follows the explanations in 
! http://homepage.univie.ac.at/franz.vesely/notes/swl/swlmix.html 
! where we generally replace uppercase letters by double lower case
! characters, thus: I_0 --> ii0, J_34 --> jj34 etc. 
!
!======  Main program ======
!
      program swl3d
!----------------------------------------------------
! Main program, only for reading in the particle and
! configuration parameters, and for calling the 
! subroutines section and swlpot
!----------------------------------------------------
      implicit real*8(a-h,o-z)
      real*8 ll1,ll2   ! meaning: L1, L2
      real*8 lam1,lam2,lam3,lam4,laml,lamu,mul,muu
      character*80 txt,txt1,txt2
      logical lovr,lint
      pi=4.d0*datan(1.d0)
!
      open(unit=5,file='swl3d.par',form='formatted') ! input parameters
!
      write(6,100) 
100   format(/1x,60('-')
     1   /1x,"SWL Mixture in 3d, arbitrary orientations:"/)  
!
! Read in particle and config parameters:
      read(5,501) txt1,ll1,ll2,s1,s2,txt2,z2,ct1,ct2,phi
      write(6,501) txt1,ll1,ll2,s1,s2,txt2,z2,ct1,ct2,phi
501   format(a/4f10.5/a/4f10.5)
! ll1,ll2 ... axis lengths
! s1,s2 ..... inner/outer SW limit)
! ct1,ct2 ... cos(theta1..2)
! phi ....... phi12 in units of pi 
!
! A typical input file 5 would look like this:
!       ll1       ll2        s1        s2
!   2.00000   3.00000   1.00000   1.50000
!        z2       ct1       ct2       phi
!   1.20000   0.00000   0.30000   0.30000
!
! Now prepare the parameters needed whenn calling the subroutines 
! section, swlpot, and sint:
! h1,h2 (half lengths); s1:
      h1=0.5d0*ll1  
      h2=0.5d0*ll2  ! half lengths
      s1sq=s1*s1
      s2sq=s2*s2    
      rsq=z2*z2
      phi=pi*phi
      ex1=sqrt(1.d0-ct1*ct1)
      ey1=0.0d0
      ez1=ct1
      st2=sqrt(1.d0-ct2*ct2)
      ex2=st2*cos(phi)
      ey2=st2*sin(phi)
      ez2=ct2
      d=ex1*ex2+ey1*ey2+ez1*ez2
      if (d.lt.0.0d0) then ! turn e2 around
        ez2=-ez2
        d=-d
      endif
      b=z2*ez2
      c=z2*ez1   
!
! Check for overlap (-> lovr=true if overlap):
      call section(rsq,s1sq,h1,h2,b,c,d,aa,bb,cc
     1 ,lam1,lam2,lam3,lam4,laml,lamu,mul,muu,lovr) 
! (Note that just now we are interested only in lovr; the values aa-cc 
! and lam1-4 etc. will be overwritten in the next call to SR section)
!
! Check for interaction (-> lint=false if no interaction),
! and determine lam_1-4, mu_l/u:
      lint=.false.
      if(.not.lovr) then
        call section(rsq,s2sq,h1,h2,b,c,d,aa,bb,cc
     1   ,lam1,lam2,lam3,lam4,laml,lamu,mul,muu,lint) 
      endif
!
! Now compute the potential:
      call swlpot(rsq,lovr,lint,b,c,d,aa,bb,cc
     1 ,lam1,lam2,lam3,lam4,laml,lamu,mul,muu,aai,aan)  
      ua=-0.25d0*aai/h1/h2  
      un=-0.25d0*aan/h1/h2  
!
      write(6,143) ua,un
143   format(/1x,"Result: potential normalized by 1/L1/L2"
     1 //10x,'ua= ',f11.6,' analytical'
     1  /10x,'un= ',f11.6,' numerical')  
!
      write(6,104) lovr,lint
104   format(10x,"overlap=     ",2x,l3
     1      /10x,"interaction= ",2x,l3,//1x,60('-'))  
!
123   continue
!
      stop
      end
!--------------------------------------------------------------
!
!======  SR section ======
!
      subroutine section(rsq,ssq,h1,h2,b,c,d,aa,bb,cc
     1 ,lam1,lam2,lam3,lam4,laml,lamu,mul,muu,lsec)
!
! Determine the intersection btw a rectangle and an ellipse
!
! Input: rsq,ssq,b,c,d
! Output: aa,bb,cc,lam1-lam4,laml,lamu,mul,muu,lsec
!
! Let R(2*h1,2*h2) be a rectangle in the y-x plane, centered at (0,0);
! let E(b,c,d,aa,bb,cc,gg,hh) be an ellipse in that plane, such that
! E: y(x) = x*d-b +- sqrt(aa*x**2 + 2*bb*x + cc)
! or
! E: x(y) = y*d+c +- sqrt(aa*y**2 + 2*gg*y + hh).
! *** Note:
! The case of parallel particles (d=1) is treated correctly. However, this 
! case is much simpler, since many tediously calculated parameters take on 
! trivial values. The strictly parallel case should therefore be treated by a 
! separate code, e.g. swlparmix, see
! http://homepage.univie.ac.at/franz.vesely/notes/swl/swlpmfortran.html 
! 
! SR "section" checks if E is real and E \sect R is non-empty: lsec=true. If so, 
! the relevant integration limits x1-x4,xl,xu and the lower/upper limits yl,yu 
! are calculated.
!
      implicit real*8(a-h,o-z)
      real*8 lam1,lam2,lam3,lam4,laml,lamu,lamns,lamps
      real*8 mul,muu,muns,muos,muqs,mups
      character*80 txt
      logical lsec
!
      aa=d*d-1.d0
      bb=c-b*d
      cc=ssq-rsq+b*b
      gg=c*d-b
      hh=ssq-rsq+c*c
! Initialize lambda, mu:
      lam1=-h1      
      lam2=-h1
      lam3=-h1      
      lam4=-h1
      laml=-h1      
      lamu= h1
      mul=-h2
      muu= h2
!
      emd2=-aa ! always >= 0.
      if(emd2.ne.0.d0) then ! nonparallel case:
!-------------------------------------------------------------------
! Nonparallel case:
! First, find out if E is real: are there any x for which y is real?
! aa*x**2+2*bb*x+cc = 0 leads to x=(1/aa)*(-bb +- sqrt(bb*bb-aa*cc))
! x is real if bb*bb>=aa*cc:
!-------------------------------------------------------------------
        if(bb*bb.lt.aa*cc)  then ! non real ellipse E
          lsec=.false.    
        else  ! real ellipse:
!--------------------------------------------------------------------------
! If E is real, check if its container rectangle Rc intersects R; 
! if no: no intersection of E and R; 
! if yes, check if xl<xu and yl<yu (==> intersection):
!
! Container rectangle Rc: points N,P (left/rightmost) and O,Q (bottom/top):
!--------------------------------------------------------------------------
          lamns=bb-sqrt(bb*bb-aa*cc)  ! =lam(N)*emd2
          lamps=bb+sqrt(bb*bb-aa*cc)  ! =lam(P)*emd2
          muos=gg-sqrt(gg*gg-aa*hh)   ! =mu(O)*emd2
          muqs=gg+sqrt(gg*gg-aa*hh)   ! =mu(Q)*emd2 
          if (lamns.gt.(-h1*emd2)) laml=min( h1,lamns/emd2)    ! clip rect at N
          if (lamps.lt.(h1*emd2))  lamu=max(-h1,lamps/emd2)   ! clip rectangle at P
          if (muos.gt.(-h2*emd2))  mul=min( h2,muos/emd2)     ! clip rect at O
          if (muqs.lt.(h2*emd2))   muu=max(-h2,muqs/emd2)     ! clip rect at Q
          if((laml.ge.lamu).or.(mul.ge.muu)) then ! empty intersection of Rc and R
            lsec=.false.
          else  ! non-empty rect intersection; there may still be no E/R inters.:
! Intersect mu=mul (lower rectangle side) with E --> points 1,2:
            arg1=hh-mul*mul*emd2+2.d0*mul*gg   
            term1=0.d0
            if (arg1.gt.0.d0) term1=sqrt(arg1)        
            lam1=min(lamu,max(laml,mul*d+c-term1))   
            lam2=min(lamu,max(laml,mul*d+c+term1))  
! Intersect mu=muu (upper rectangle side) with E --> points 3,4:
            arg1=hh-muu*muu*emd2+2.d0*muu*gg   
            term1=0.d0
            if (arg1.gt.0.d0) term1=sqrt(arg1)         
            lam3=min(lamu,max(laml,muu*d+c-term1))   
            lam4=min(lamu,max(laml,muu*d+c+term1))  
! Clip rectangle further: 
! if pt. 2 or 4 is on upper branch of Ei, lamu=lam2 or lam4:
            mups=gg+d*sqrt(bb*bb-aa*cc)   ! =mu(P)*emd2 
            if (mul*emd2.gt.mups) lamu=lam2    
            if (muu*emd2.le.mups) lamu=lam4  
! if pt. 1/3 on lower branch of Ei, xl=x1/x3 etc.
! Same with pts. 1,3:
            muns=gg-d*sqrt(bb*bb-aa*cc)   ! =mu(N)*emd2  
            if (mul*emd2.gt.muns) laml=lam1  
            if (muu*emd2.le.muns) laml=lam3  
!
            lsec=.true. 
            if(laml.ge.lamu) lsec=.false.
          endif ! end non-empty interval
        endif  ! end real ellipse
      else  ! end nonparallel case
!----------------------------------------------------------------------
! Parallel case: ellipse is now a pair of lines,
! y(x)=x-b +- sqrt(C); therefore N,P and O,Q are somewhere in infinity,
! and xl, xu as well as yl, yu are -+h.
! First, find out if E is real: are there any x for which y is real?
! aa*x**2+2*bb*x+cc > 0? x=(1/aa)*(-bb +- sqrt(bb*bb-aa*cc)) real
! if cc >=0:
!-----------------------------------------------------------------------
        if(cc.lt.0.0d0)  then ! no real pair of lines
          lsec=.false.    ! E is not real
        else  ! real pair of lines
! Intersect mu=mul (lower rectangle side) with E --> points 1,2:
          arg1=hh
          term1=0.0
          if (arg1.gt.0.0d0) term1=sqrt(arg1)      
          lam1=min(lamu,max(mul*d+c-term1,laml)) 
          lam2=min(lamu,max(mul*d+c+term1,laml))
! Intersect mu=muu (upper rectangle side) with E --> points 3,4:
          lam3=min(lamu,max(muu*d+c-term1,laml)) 
          lam4=min(lamu,max(muu*d+c+term1,laml))
! Clip rectangle further: we certainly have laml=lam1, and lamu=lam4:
          laml=lam1
          lamu=lam4 
!
          lsec=.true.
          if(laml.ge.lamu) lsec=.false.
        endif  ! end real p.o.lines
      endif ! end parallel case
!
      return 
      end 
!--------------------------------------------------------------
!
!======  SR swlpot ======
!
      subroutine swlpot(rsq,lovr,lint,b,c,d,aa,bb,cc
     1 ,lam1,lam2,lam3,lam4,laml,lamu,mul,muu,aai,aan)  
!------------------------------------------
! Compute SWL potential analytically (ua)
! and compare to a quick numerical sum (un)
!
! Input: rsq,lovr,lint,b,c,d,lam1-lam4,laml,lamu,mul,muu
! Output: ua,un
!------------------------------------------
      implicit real*8(a-h,o-z)
      real*8 lam1,lam2,lam3,lam4,laml,lamu,mul,muu
      real*8 lami,muil,muiu
      real*8 ii0,ii12,ii34,iilup,iilum ! meaning the integrals I_0, I_12, I_34
! meaning: lambda_1-4, lambda_lower/upper, mu_lower/upper                 
      character*8 text
      logical lovr,lint
!
      aan=0.d0   ! area of intersection (numerical)
      aai=0.d0   ! area of intersection (analytical)
      if(.not.lovr.and.lint) then ! start calculating the potential
!
! Compute the integrals ii0, ii12, ii34 that make up 
! the intersection area of rectangle and ellipse:
!
        call sint(aa,bb,cc,b,d,laml,lamu,1,mul,iilup)  ! I_lower_upper^(+)
        call sint(aa,bb,cc,b,d,laml,lamu,-1,mul,iilum) ! I_lower_upper^(-)
        ii0=iilup-iilum  
        call sint(aa,bb,cc,b,d,lam1,lam2,-1,mul,ii12)  ! I_12
        call sint(aa,bb,cc,b,d,lam3,lam4,1,muu,ii34)   ! I_34
! Now combine these integrals:
        aai=ii0+ii12-ii34 ! aai = A is the intersection area 
!
! check: calculate a simple numerical estimate un for the potential:
        nlam=10000
        dlam=(lamu-laml)/(float(nlam-1))
        do 1 i=1,nlam
          lami=laml+float(i-1)*dlam 
          arg2=aa*lami*lami+2.0d0*lami*bb+cc   
          xincr=0.d0
          term1=0.d0
          if(arg2.gt.0.0d0) then
            term1=sqrt(arg2)
! E: y(x) = x*d-b +- sqrt(aa*x**2 + 2*bb*x + cc)
            muiu=max(mul,min(muu,lami*d-b + term1))
            muil=max(mul,min(muu,lami*d-b - term1))
            xincr=muiu-muil  ! meaning mu_i_upper - mu_i_lower at that lam_i
          endif
          aan=aan+xincr  
1       continue
        aan=aan*dlam
      endif
!
      return 
      end 
!-----------------------------------------------------------
!
!======  SR sint ======
!
      subroutine sint(aa,bb,cc,b,d,lama,lamb,moda,base,iiab)
!---------------------------------------------------------------------
! Do integrations to compute the overlap area of rectangle and ellipse
!
! Input: aa,bb,cc,dd,lama,lamb,moda,base
! Output: iiab (may be I_lower_upper, I_12, or I_34) 
!---------------------------------------------------------------------
      implicit real*8(a-h,o-z)
      real*8 lama,lamb,iiab,jja,jjb,jjab 
!
      del=aa*cc-bb*bb
      delsq=del*del 
      xkk=1.d0
      emd2=-aa
      emd=sqrt(emd2)
      if(delsq.gt.0.0d0) xkk=aa/del  
! value at limit lamb:
      arg1=aa*lamb*lamb+2.d0*bb*lamb+cc   
      term1=0.0d0
      if(arg1.ge.0.d0) term1=sqrt(arg1)  
      arg2=1.d0
      if(del.lt.0.d0.and.delsq.gt.0.0d0) then
        arg2=(aa*lamb+bb)/sqrt(-del)
        if(arg2.gt.1.0d0) arg2=1.0d0
        if(arg2.lt.-1.0d0) arg2=-1.0d0
      endif  
      term2=0.5d0*pi
      if(arg2.le.1.0d0.and.arg2.ge.-1.0d0) term2=asin(arg2)  
      jjb=0.d0
      if(cc.gt.0.d0) jjb=sqrt(cc)*lamb  ! parallel case  
      if(aa.lt.0.0d0) then
        bra1=(aa*lamb+bb)*term1/2.0d0/aa 
        bra2=-0.5d0*term2/xkk/emd   
        jjb=bra1+bra2   
      endif
! value at limit lama:
      arg1=aa*lama*lama+2.0d0*bb*lama+cc   
      term1=0.d0
      if(arg1.gt.0.d0) term1=sqrt(arg1)  
      arg2=0.d0
      if(del.lt.0.d0.and.delsq.gt.0.0d0) arg2=(aa*lama+bb)/sqrt(-del)  
      term2=0.5d0*pi   
      if(arg2.lt.1.0d0) term2=asin(arg2)  
      jja=sqrt(cc)*lama  
      if(aa.lt.0.0d0) then
        bra3=(aa*lama+bb)*term1/2.0d0/aa
        bra4=-0.5d0*term2/xkk/emd   
        jja=bra3+bra4   
      endif
      jjab=jjb-jja     ! this is the pure sqrt, integrated
      iiab=0.5d0*d*(lamb*lamb-lama*lama)
     1     -(b+base)*(lamb-lama)+moda*jjab    
!
      return
      end
!----------------------------------------------------------

