program adiabatic_rotor
c     written by G. Ballentine, G. Bertsch, N. Onishi and K. Yabana
c     based on the theory of Bertsch, Onishi and Yabana, Zeitschrift
c     fuer Physik D34 (1995) 213.  This program may be freely used
c     providing such use cites the source,
c     G. Ballentine, G. Bertsch, N. Onishi and K. Yabana, (2005) to be
c     published.
c     program should output <ubar> = 0.192
      implicit none
      integer Imax,Mzmax,M3max,I,Mz,M3,KECH,iter,iet,intu,n
      real*8 pi,eps,rImax,dI,dMz,dM3,betamu0B,J1J3,ht,zsum,pol
      real*8 rI,rMz,rM3,weight,Jtheta0,et,a,b,c,d,u0,u1,u2,Jtheta
      real*8 nu,taut,ellp1,ellp2,Jint,ubar,uhist,uu,emin,umin,usum

c note that Mzmax =/ M3max, to remove multiple solutions
      parameter(Imax=200,Mzmax=200,M3max=201)
      parameter(rImax=4.d0,eps=5d-3,pi=3.141592653589793d0)
      dimension uhist(100)

c there are two basic inputs, x=beta*mu_0*B, and J1/J3
c beta*mu_0*B
      betamu0B=1.0d0
      if(betamu0B.lt.0.001d0) stop 'parameter out of range: x < 0.001'
c J1/J3
      J1J3=1.d0
      do n=1,100
        uhist(n)=0.d0
      enddo

      dI=rImax/dfloat(Imax)
      dMz=2.d0/dfloat(Mzmax)
      dM3=2.d0/dfloat(M3max)
      ht=betamu0B

      zsum=0.d0
      pol=0.d0
      do I=1,Imax
c        write(*,*) 'I=',I
        rI=I*dI
c set Mzmax even to limit Mz < 0
      do Mz=1,Mzmax/2
        rMz=-1.d0+(Mz-0.5d0)*dMz
      do M3=1,M3max
        rM3=-1.d0+(M3-0.5d0)*dM3

c rMz =/ rM3
        if(abs(rMz-rM3).lt.1d-5) then
          write(*,*) 'abs(rMz-rM3) is too small'
          stop
        endif

        weight=rI**2*exp(-rI**2*((1-rM3**2)+J1J3*rM3**2))
        zsum=zsum+weight
        Jtheta0=2*pi*rI*(1-abs(rMz-rM3)/2.d0-abs(rMz+rM3)/2.d0)

c adiabatic connection between E and I

c set et at minimum potential energy 
        et=1d10
        do n=1,100
         uu=-1+(n-0.5d0)/50.d0
         emin=rI**2*(rMz-rM3*uu)**2/(1-uu*uu)-ht*uu
         if(emin.le.et) then
            et=emin
            umin=uu
         endif
        enddo
        et=et+1d-3

        iet=0
   10   continue
        a=-ht
        b=-et-rI**2*rM3**2
        c=ht+2*rI**2*rMz*rM3
        d=et-rI**2*rMz**2
        call CBCEQ(A,B,C,D,u0,u1,u2,KECH)
c if not 3 roots, increase energy
        if(KECH.eq.1) then
          write(*,*) 'warning: no three real solutions',umin,et
          et=et+0.01d0
          iet=iet+1
          if(iet.eq.5) stop 'failto find three real solutions'
          go to 10
        endif

c newton method to find energy
        do iter=1,10
          call Jint_sub(u0,u1,u2,Jint)
          Jtheta=2*sqrt(ht)*Jint
          nu=(u2-u1)/(u2-u0)
          taut=ellp1(nu)*2/sqrt(ht*(u2-u0))
          if(abs(Jtheta0-Jtheta).lt.eps) exit
          et=et+(Jtheta0-Jtheta)/taut
          iet=0
   20     continue
          b=-et-rI**2*rM3**2
          d=et-rI**2*rMz**2
          call CBCEQ(A,B,C,D,u0,u1,u2,KECH)
          if(KECH.eq.1) then
           write(*,*) 'warning: no three real solutions,iet=',iet,et
           et=et+1d-4
           iet=iet+1
           if(iet.eq.5) stop 'failto find three real solutions'
           go to 20
          endif
          if(iter.eq.10) then
            write(*,120) I,Mz,M3,abs(Jtheta-Jtheta0)
 120        format(1x,'et not converged',3i5,f15.7)
          endif
        enddo
c u0,u1,u2 are obtained.

        ubar=u0+(u2-u0)*ellp2(nu)/ellp1(nu)
        pol=pol+weight*ubar
        intu=(ubar+1)*50.d0+1
        if(intu.le.0.or.intu.gt.100) write(*,*) 'intu out of range'
c histgram
        uhist(intu)=uhist(intu)+weight

      end do
      end do
      end do

      zsum=zsum*dI*dMz*dM3
c average polarization
      pol=pol*dI*dMz*dM3/zsum 
      write(*,15) betamu0B,J1J3,pol
 15   format(2f8.3,f10.5)
      write(*,*)
      do n=1,100
        usum=usum+uhist(n)
      enddo
      open(7,file='udist.dat',form='formatted')
      do n=1,100
       write(7,100) (n-0.5d0)*0.02d0-1.d0, uhist(n)/usum/0.02d0
 100   format(1x,2f15.6)
      enddo
      close(7)

      stop
      end

c complete elliptic integral of 1st kind, K(k)

       FUNCTION ELLP1(FM)
        IMPLICIT REAL*8 (A-H,O-Z)
        DATA A0,A1/ 1.38629436112D0, 0.09666344259D0/
        DATA A2,A3/ 0.03590092383D0, 0.03742563713D0/
        DATA A4,B0/ 0.01451196212D0, 0.5D0          /
        DATA B1,B2/ 0.12498593597D0, 0.06880248576D0/
        DATA B3,B4/ 0.03328355346D0, 0.00441787012D0/
        FM1=1.0D0-FM
        F1IN=1.0D0/FM1
        ELLP1=A0+FM1*(A1+FM1*(A2+FM1*(A3+A4*FM1)))
     1      +(B0+FM1*(B1+FM1*(B2+FM1*(B3+B4*FM1))))*DLOG(F1IN)
       RETURN
       END

c complete elliptic integral of 2nd kind, E(k)

       FUNCTION ELLP2(FM)
        IMPLICIT REAL*8 (A-H,O-Z)
        DATA A1/    0.44325141463D0/
        DATA A2,A3/ 0.06260601220D0, 0.04757383546D0/
        DATA A4   / 0.01736506451D0/
        DATA B1,B2/ 0.24998368310D0, 0.09200180037D0/
        DATA B3,B4/ 0.04069697526D0, 0.00526449639D0/
        FM1=1.0D0-FM
        F1IN=1.0D0/FM1
        ELLP2=1.0D0+FM1*(A1+FM1*(A2+FM1*(A3+A4*FM1)))
     1      + FM1*(B1+FM1*(B2+FM1*(B3+B4*FM1)))*DLOG(F1IN)
       RETURN
       END

c integral related to the action J_theta
      subroutine Jint_sub(u0,u1,u2,Jint)
      implicit none
      integer i,Nmax
      real*8 u0,u1,u2,Jint,pi,x1,x2,dx,s,x,u
      real*8 ubar,udif,theta
      parameter(Nmax=200,pi=3.141592653589793d0)
      if(u1.le.-1.d0.or.u2.ge.1.d0) stop 'u1 or u2 out of range'
      if(u1.lt.-0.9d0.and.u2.gt.0.9d0) then
        x1=acos(u2) 
        x2=acos(u1)
        if(x2.lt.x1) stop
        dx=(x2-x1)/Nmax
        s=0.d0
        do i=1,Nmax
          x=x1+(i-0.5d0)*dx
          u=cos(x)
          s=s+sqrt((u0-u)*(u1-u)*(u2-u))/sin(x)
        enddo
        Jint=s*dx
      else if(u1.lt.-0.9d0) then
        x1=sqrt(1+u1)
        x2=sqrt(1+u2)
        dx=(x2-x1)/Nmax
        s=0.d0
        do i=1,Nmax
         x=x1+(i-0.5d0)*dx
         u=x*x-1
         s=s+sqrt((u0-u)*(u1-u)*(u2-u))/(1-u*u)*2*x
        enddo
        Jint=s*dx
      else if(u2.gt.0.9d0) then
        x1=sqrt(1-u2)
        x2=sqrt(1-u1)
        dx=(x2-x1)/Nmax
        s=0.d0
        do i=1,Nmax
         x=x1+(i-0.5d0)*dx
         u=1-x*x
         s=s+sqrt((u0-u)*(u1-u)*(u2-u))/(1-u*u)*2*x
        enddo
        Jint=s*dx
      else
        ubar=(u2+u1)/2.d0
        udif=(u2-u1)/2.d0
        s=0.d0
        do i=1,Nmax
          theta=(i-0.5d0)*pi/Nmax-(pi/2)
          u=ubar+udif*sin(theta)
          s=s+cos(theta)**2*sqrt(u-u0)/(1-u*u)
        enddo
        Jint=s*(pi/Nmax)*udif**2
      endif
      return
      end

c solve 3rd order algebraic equation ax^3+bx^2+cx+d=0
c if a=0 stop
c if three real solutions, KECH=0 and x1<x2<x3
c if less than three real solutions, KECH=1

      SUBROUTINE CBCEQ(A,B,C,D,X1,X2,X3,KECH)
      implicit real*8(a-h,o-z)
      data pi/3.141592653589793d0/
      if(abs(a).lt.1d-10) stop 'a is too small'
      p  = c/a - b*b/a/a/3.d0
      q  = (2.d0*b*b*b/a/a/a - 9.d0*b*c/a/a + 27.d0*d/a) / 27.d0
      DD = p*p*p/27.d0 + q*q/4.d0
      if(DD.ge.0.d0) then
       KECH=1
       return
      endif
      KECH=0
      phi = acos(-q/2.d0/sqrt(abs(p*p*p)/27.d0))
      temp1=2.d0*sqrt(abs(p)/3.d0)
      x1 =  temp1*cos(phi/3.d0)-b/a/3.d0
      x2 = -temp1*cos((phi+pi)/3.d0)-b/a/3.d0
      x3 = -temp1*cos((phi-pi)/3.d0)-b/a/3.d0
      if(x1.gt.x2) then
       x0=x1
       x1=x2
       x2=x0
      endif
      if(x1.gt.x3) then
       x0=x1
       x1=x3
       x3=x0
      endif
      if(x2.gt.x3) then
       x0=x2
       x2=x3
       x3=x0
      endif
      if(x1.ge.-1.d0) write(*,*) 'warning u0.ge.-1',x1
      if(x2.le.-1.d0.or.x3.ge.1.d0) then
      write(*,*) x1,x2,x3
      pause
      endif
      return
      end