!
! STATUS: Oct 28, 2010
!
!------------------------------------------------------------
! swlpmint.f ("SWL, parallel, mixture, integration over rho")
!
! () For various z displacements btw two square well lines,
!    do the area integral J(z) = 2*pi*int[rho*u(rho,z)*drho]
! () compute the Discrete Fourier Transform F(n) of J(z)
!
! The code follows the explanations given in
! http://homepage.univie.ac.at/franz.vesely/notes/swl/swlparmix.html#areaint
! where we generally replace uppercase letters by double lower case
! characters, thus: J_0 --> jj0, F(n) --> ff(n) etc. 
!-------------------------------------------------
!
!======  Main program ======
!
      program swlpmint
!----------------------------------------------------------------------------
! Main program, only for reading in the particle parameters, then
! () calling SR "jayofzee" which produces J(z) for a number of z values;
! () calling SR "dfourier" to compute the discrete Fourier transform of J(z) 
!----------------------------------------------------------------------------
!     
      implicit real*8(a-h,o-z)
      real*8 jj,jj0,jj12,jj34
      real*8 ll1,ll2 ! particle lengths
      dimension fk(0:1200),fkt(0:1200),ffr(0:1200),ffi(0:1200)
     1 ,fnr(-600:600),fni(-600:600)
      character*80 txt1
!
      pi=4.d0*datan(1.d0)         
      pi2=2.d0*pi
!
      open(unit=5,file='swlpmint.par',form='formatted')
      open(unit=8,file='swlpmint.dat',form='formatted')
      open(unit=9,file='swlpmint.rhi',form='formatted')
      open(unit=11,file='swlpmint.joz',form='formatted')
      open(unit=12,file='swlpmint.dft',form='formatted')
      open(unit=99,file='swlpmint.tst',form='formatted')
!
      write(6,100) 
100   format(1x,"SW Lines, parallel, mixed:")
!
!--------------------------------------------------
! Read in particle parameters:
      read(5,501) txt1,ll1,ll2,s1,s2,nz
      write(6,501) txt1,ll1,ll2,s1,s2,nz
501   format(a/4f10.5,i10)
! ll1,ll2 ... axis lengths
! s1,s2 ..... inner/outer SW limit)
! nz ........ (even) no. of z intervals btw 0 and zmax
!
! A typical input file 5 would look like this:
!       ll1       ll2        s1        s2        nz
!   2.00000   3.00000   1.00000   1.50000       100
!
!--------------------------------------------------
!
      if(ll2.lt.ll1) stop 66 ! see website for explanation
! 
      h1=0.5d0*ll1      
      h2=0.5d0*ll2
      s1sq=s1*s1      
      s2sq=s2*s2       
!
!   Configuration: Particle 1 at origin
      x1=0.d0      
      z1=0.d0
!
! z12 loop:
      zmax=h1+h2+s2      
      dz=0.d0      
      if (nz.gt.0) dz=zmax/float(nz)
      do 98 k=0,nz
        z=k*dz
        call jayofzee(z,h1,h2,s1sq,s2sq,jj)
        fk(k)=jj      
        write(6,876) k,z,fk(k)
        write(11,876) k,z,fk(k) ! table entry [k, z(k), J(z_k)] for FT
876     format(1x,' k,z,-2*pi*sj/L1/L2=',i5,7f10.5)
!
98    continue
!
! Now do the FT of f_k:
      call dfourier(nz,dz,fk,fnr,fni) ! Discrete Fourier Transform
!
      stop
      end
!---------------------------------------------------------------------
!
!======  SR jayofzee ======
!
      subroutine jayofzee(z,h1,h2,s1sq,s2sq,jj)
!-----------------------------------------------------------
!  Compute area integrals J0(z),J12(z),J34(z) and total J(z)
!-----------------------------------------------------------
      implicit real*8(a-h,o-z)
      real*8 jj,jj0,jj12,jj34
      real*8 jj0l,jj0u,jj12l,jj12u,jj34l,jj34u
      pi=4.d0*datan(1.d0)         
!
! Determine boundaries of rho intervals:
        f=z-h2+h1
        f0=min(0.d0,f)
        f1=max(0.d0,f)       
        e=z-h2-h1
        e0=min(0.d0,e)
        e1=max(0.d0,e)      
        v=z+h2+h1      ! always >0 
        u=z+h2-h1      ! always >0 
        rhmnsq=max(0.d0,s1sq-e1*e1)
        rhmn=dsqrt(rhmnsq)
        rhmxsq=max(0.d0,s2sq-e1*e1)
        rhmx=dsqrt(rhmxsq)
        rh1sq=min(rhmxsq,max(rhmnsq,s2sq-f1*f1))
        rh1=dsqrt(rh1sq)
        rh2sq=min(rhmxsq,max(rhmnsq,s2sq-e1*e1))
        rh2=dsqrt(rh2sq)
        rh3sq=min(rhmxsq,max(rhmnsq,s2sq-f0*f0))
        rh3=dsqrt(rh3sq)
        rh4sq=min(rhmxsq,max(rhmnsq,s2sq-e0*e0))
        rh4=dsqrt(rh4sq)
        rh5sq=min(rhmxsq,max(rhmnsq,s2sq-v*v))
        rh5=dsqrt(rh5sq)
        rh6sq=min(rhmxsq,max(rhmnsq,s2sq-u*u))
        rh6=dsqrt(rh6sq)
!
!        write(9,704) z,rhmn,rh1,rh2,rh3,rh4,rh5,rh6,rhmx
704     format(1x," z,rhmn,rh1-6,rhmx=",9f10.5)
!        write(6,602) ll1,ll2,z
602     format(//1x,"Area integral jj(z) of u(z;rho) for L1,L2=",2f8.2
     1  /12x,'   z=',f10.5)
!
      jj0=0.d0
      jj12=0.d0
      jj34=0.d0      
!
! Compute J0:
!
! 1A-4C: J0=-4/3 h1 t^3      
!        I0=4.d0*t*h1
      rhl=rhmn
      tl=dsqrt(s2sq-rhl*rhl)
      jj0l=-(4.d0/3.d0)*h1*tl*tl*tl
      rhu=rh1
      tu=dsqrt(s2sq-rhu*rhu)
      jj0u=-(4.d0/3.d0)*h1*tu*tu*tu
      jj0=jj0+(jj0u-jj0l)
!
! 1B-4C: J0=-4/3*h1*t^3      
!        I0=-2*t*e+2*(s2sq-rhsq)
      rhl=rh1
      rhlsq=rhl*rhl
      tl=dsqrt(s2sq-rhlsq)
      jj0l=s2sq*rhlsq+(2.d0/3.d0)*e*tl*tl*tl-0.5d0*rhlsq*rhlsq
      rhu=rhmx
      rhusq=rhu*rhu
      tu=dsqrt(s2sq-rhusq)
      jj0u=s2sq*rhusq+(2.d0/3.d0)*e*tu*tu*tu-0.5d0*rhusq*rhusq
      jj0=jj0+(jj0u-jj0l)
!
!------------
! Compute s12:
!
! 1A-2B: J12=-(1/4)(s2^2+f^2)rh^2+(1/3)*f*t^3+(1/8)rh^4
!        I12=-0.5*(s2sq+f^2)-f*t+0.5*rhsq
      rhl=rh4
      rhlsq=rhl*rhl
      tl=dsqrt(s2sq-rhlsq)
      jj12l=-0.25d0*(s2sq+f*f)*rhlsq+(1.d0/3.d0)*f*tl*tl*tl
        jj12l=jj12l+0.125d0*rhlsq*rhlsq
      rhu=max(rhl,min(rh1,rh3))       
      rhusq=rhu*rhu
      tu=dsqrt(s2sq-rhusq)
      jj12u=-0.25d0*(s2sq+f*f)*rhusq+(1.d0/3.d0)*f*tu*tu*tu
        jj12u=jj12u+0.125d0*rhusq*rhusq
      jj12=jj12+(jj12u-jj12l)
!
! 1A-2C: J12=-h1*(z-h2)*rh^2+(2/3)*h1*t^3      
!        I12=-2*h1*(z-h2)-2*h1*t
      rhl=rhmn
      tl=dsqrt(s2sq-rhl*rhl)
      jj12l=-h1*(z-h2)*rhl*rhl+(2.d0/3.d0)*h1*tl*tl*tl
      rhu=min(rh1,rh4)
      tu=dsqrt(s2sq-rhu*rhu)
      jj12u=-h1*(z-h2)*rhu*rhu+(2.d0/3.d0)*h1*tu*tu*tu
      jj12=jj12+(jj12u-jj12l)
!
! 1B-2B: J12=- s^2*rh^2 + 1/2*rh^4
!        I12=-2*s2sq+2*rhsq
      rhl=max(rh1,rh4)
      tl=dsqrt(s2sq-rhl*rhl)
      jj12l=-s2sq*rhl*rhl+(0.5d0)*rhl*rhl*rhl*rhl
      rhu=rh3
      if(rhu.lt.rhl) rhu=rhl
      tu=dsqrt(s2sq-rhu*rhu)
      jj12u=-s2sq*rhu*rhu+(0.5d0)*rhu*rhu*rhu*rhu
      jj12=jj12+(jj12u-jj12l)
!
! 1B-2C: J12=- (1/4) (3 s2^2 - e^2) rh^2 - (1/3) e t^3 + (3/8) rh^4
!        I12=0.5*e*e+e*t-1.5*s2sq+1.5*rhsq      
      rhl=rh1
      tl=dsqrt(s2sq-rhl*rhl)
      jj12l=-0.25d0*(3.d0*s2sq-e*e)*rhl*rhl-(1.d0/3.d0)*e*tl*tl*tl
        jj12l=jj12l+(3.d0/8.d0)*rhl*rhl*rhl*rhl
      rhu=max(rhl,rh4)
      tu=dsqrt(s2sq-rhu*rhu)
      jj12u=-0.25d0*(3.d0*s2sq-e*e)*rhu*rhu-(1.d0/3.d0)*e*tu*tu*tu
        jj12u=jj12u+(3.d0/8.d0)*rhu*rhu*rhu*rhu
      jj12=jj12+(jj12u-jj12l)
!
!---------------------------
! Compute s34:
!
! 3A-4C: J34=- h1*(z + h2)*rh^2 - 2/3*h1*t^3
!        I34=-(z+h2-t)*2*h1      
      rhl=rhmn
      tl=dsqrt(s2sq-rhl*rhl)
      jj34l=-h1*(z+h2)*rhl*rhl-(2.d0/3.d0)*h1*tl*tl*tl
      rhu=max(rhl,rh5)
      tu=dsqrt(s2sq-rhu*rhu)
      jj34u=-h1*(z+h2)*rhu*rhu-(2.d0/3.d0)*h1*tu*tu*tu
      jj34=jj34+(jj34u-jj34l)
!
! 3B-4C: J34=(1/4) (u^2 + s2^2) rh^2 + 1/3 u t^3
!        I34=0.5*u*u-u*t+0.5*s2sq-0.5*t^2     
      rhl=rh5
      tl=dsqrt(s2sq-rhl*rhl)
      jj34l=0.25d0*(u*u+s2sq)*rhl*rhl+(1.d0/3.d0)*u*tl*tl*tl
        jj34l=jj34l-0.125d0*rhl*rhl*rhl*rhl
      rhu=max(rhl,rh6)
      tu=dsqrt(s2sq-rhu*rhu)
      jj34u=0.25d0*(u*u+s2sq)*rhu*rhu+(1.d0/3.d0)*u*tu*tu*tu
        jj34u=jj34u-0.125d0*rhu*rhu*rhu*rhu
      jj34=jj34+(jj34u-jj34l)
!
      jj=-2.d0*pi*(jj0+jj12-jj34) ! area integral of rho*u(rho,z) 
      jj=0.25d0*jj/h1/h2          ! normalize to particle lengths
!
      return
      end
!----------------------------------------------------------------
!
!======  SR dfourier ======
!
      subroutine dfourier(nz,dz,fk,fnr,fni)
!----------------------------------------------------------------
! Do the Discrete Fourier Transform of a table fk(k=0;...;N-1)
! Fn = sum_k=0^N-1 fk e^2*pi*i*k*n/N  with N even ;      and n=0;+-1;:::;+-N/2
! For definitions see:
! http://homepage.univie.ac.at/franz.vesely/ves/cp_tut/nol2h/new/c9ap_s2fft.html
!----------------------------------------------------------------
      implicit real*8(a-h,o-z)
      dimension fk(0:1200),fkt(0:1200),ffr(0:1200),ffi(0:1200)
     1 ,fnr(-600:600),fni(-600:600)
!
! fk .... table of J(z_k)
! fkt ... table of J(z_k), reconstructed from discr. FT
! ffr,ffi ... discrete Fourier coefficients Fn, real and imaginary
! fnr,fni ... discrete Fourier coefficients Fn, ordered from -nz/2 to nz/2
!    
      pi=4.d0*datan(1.d0)         
      pi2=2.d0*pi
!
      do 122 ng=-nz/2,nz/2
        fnr(ng)=0.d0
        fni(ng)=0.d0
122   continue
!
      do 12 ng=-nz/2,nz/2
        do 24 k=0,nz-1
          arg=pi2*float(k)*float(ng)/float(nz)
          ckn=cos(arg)
          skn=sin(arg)
          fnr(ng)=fnr(ng)+fk(k)*ckn
          fni(ng)=fni(ng)+fk(k)*skn
24      continue
!        write(12,875) ng,fnr(ng),fni(ng)
875     format(1x,' ng,fnr(ng),fni(ng)=',i5,7f10.5)
12    continue 
!---------------------------------------------
! As a test, reproduce fk from the DFTransform
!
! First, renumber the gg(ng) --> ff(n):
      do 25 n=0,nz-1
        ffr(n)=fnr(n)
        ffi(n)=fni(n)
        if(n.gt.nz/2) then
          ffr(n)=fnr(n-nz)
          ffi(n)=fni(n-nz)
        endif
        write(12,976) n,ffr(n),ffi(n)
976     format(1x,' n,ffr(n),ffi(n)=',i5,7f10.5)
25    continue
!
! Now do the reconstruction:
      do 13 k=0,nz-1
        fkt(k)=0.d0
        do 26 n=0,nz-1
          arg=pi2*float(k)*float(n)/float(nz)
          ckn=cos(arg)
          skn=sin(arg)
          fkt(k)=fkt(k)+(ffr(n)*ckn+ffi(n)*skn)
26      continue
        fkt(k)=fkt(k)/float(nz)
        write(99,874) k,fkt(k),fk(k) ! test: fkt=fk at all grid pts? 
874     format(1x,' k,frr,fk=',i5,7f10.5)
13    continue    
!
      return
      end
