!! . . This module calculates the Jacobian integrals
!! . . and elliptic Jacobian functions required for
!! . . propagation in the Ellipsoidal geometry
module elliptic

  use sep
  integer,private :: chunk
  real, private     :: tol,pi
  real, allocatable,private :: aaa(:),bbb(:),ccc(:)

contains
  !!----------------------------------------------------------------
  subroutine elliptic_setup(chunk_in,tol_in)
    integer :: chunk_in
    real  :: tol_in
    pi = acos(-1.)
    chunk = chunk_in
    tol = tol_in

    allocate ( aaa(chunk), bbb(chunk), ccc(chunk) )

  end subroutine elliptic_setup
  !!----------------------------------------------------------------

  subroutine elliptic12(u,m,F,E,Z)
    !! ELLIPTIC12 evaluates the value of the Incomplete Elliptic Integrals 
    !! of the First, Second Kind and Jacobi's Zeta Function.
    !!
    !!   call elliptic12(u,m,tol,chunk,F,E,Z) where U is a phase in radians, 0<M<1 is 
    !!   the module and TOL is the tolerance (optional). 
    !!
    !!   ELLIPTIC12 uses the method of the Arithmetic-Geometric Mean 
    !!   and Descending Landen Transformation to determine the value of the 
    !!   Incomplete Elliptic Integrals of the First, Second Kind and Jacobi Zeta Function.
    !!
    !!       F(phi,m) = int(1/sqrt(1-m*sin(t)^2), t=0..phi);
    !!       E(phi,m) = int(sqrt(1-m*sin(t)^2), t=0..phi);
    !!       Z(phi,m) = E(u,m) - E(m)/K(m)*F(phi,m).
    !!
    integer :: signU,ii,n,i,ee,NN
    real    :: F,E,Z,m,u,mu,phin,CC,Cp,Ff,tmp

    !! Initialize input and outputs
    F = 0.; E = 0.; Z = 0.; 

    if (m.gt.0. .and. m.lt.1.) then

       signU = sign(1.,u);

       do ii=1,chunk
          aaa(ii)=0. 
          bbb(ii)=0.
          ccc(ii)=0.
       end do

       aaa(1) = 1. 
       ccc(1) = sqrt(   m)
       bbb(1) = sqrt(1.-m)
       n = 0 

       !! . . Arithmetic-Geometric Mean of A, B and C
       i = 1
       do while (abs(ccc(i)) .gt. tol)  
          i = i + 1;
          aaa(i) =      (aaa(i-1) + bbb(i-1))/2.
          bbb(i) =  sqrt(aaa(i-1) * bbb(i-1))
          ccc(i) =      (aaa(i-1) - bbb(i-1))/2.
          if (ccc(i) .le. tol .or. ccc(i-1) .gt. tol ) then   
             n = (i-1)
          end if
       end do

       phin= 0.
       CC  = 0.
       Cp  = CC
       ee  = floor(CC)
       phin= signU*u
       i   = 0
       !! . . Descending Landen Transformation 
       do while (i.lt.n) 
          i = i + 1
          !if (i .lt. n) then
          phin = phin - atan(bbb(i)/aaa(i)*tan(phin)) + pi*(ceiling(phin/pi-0.5))
          ee  = 2**(i-1)
          CC = CC + dble(ee)*ccc(i)*ccc(i)
          Cp = Cp + ccc(i+1)*sin(phin)  
          !end if
       end do

       Ff = phin / ( 2.0*aaa(n)*dble(ee) );                                                      
       F  = signU*Ff                  !! Incomplete Ell. Int. of the First Kind
       Z  = signU*Cp                  !! Jacobi Zeta Function
       E  = signU*(Cp+(1.-0.5*CC)*Ff) !! Incomplete Ell. Int. of the Second Kind

    else if (m.eq.0.) then !! . . Case where M=0
       F = u 
       E = u 
       Z = 0 
    else if (m.eq.1.) then !! . . Case where M=1

       tmp = (abs(u) + pi/2.) /pi 

       !! some kind of error here ...
       !! . . using NN=tmp instead of floor(tmp)
       !!NN = floor(tmp)

       if ( abs(u).gt.pi/2.) then
          F = log(tan(pi/4. + u/2.))
       else 
          F = huge(u)*sign(1.,u)
       end if

       E = ((-1)**NN * sin(abs(u)) +2.*NN)*sign(1.,u)
       Z =  (-1)**NN * sin(u  )                      
       !       E = ((-1)**NN * sin(abs(u)) +2.*NN)*sign(1.,u)
       !       Z =  (-1)**NN * sin(u  )     
    else 
       call seperr('M must be in the range 0 <= M <= 1.')
    end if
  end subroutine elliptic12

  !!----------------------------------------------------------------

  subroutine ellipj(u,m,sn,cn,dn,am)
    !! ellipj: Jacobi elliptic functions and Jacobi amplitude.
    !! . .  [Sn,Cn,Dn,Am] = ELLIPJ(U,M) returns the values of the Jacobi
    !! . .  elliptic functions sn, cn, dn and am evaluated for corresponding
    !! . .  elements of argument u and parameter m where m is limited to [0.,1.] 
    !!
    !! . .  call ellipj(u,m,tol,chunk,sn,cn,dn,am) computes the elliptic functions to
    !! . .  the accuracy tol.  Some definitions of the Jacobi elliptic functions use 
    !! . .  the modulus k instead of the parameter m.  They are related by m = k^2.
    !!
    integer:: ii,i,n
    real   :: u,m,sn,cn,dn,am,phin

    am = 0.; cn = 0.; sn = 0.; dn = 0.

    if (m .gt. 0. .and. m .lt. 1.) then

       do ii=1,chunk
          aaa(ii)=0.
          bbb(ii)=0.
          ccc(ii)=0.
       end do

       aaa(1) = 1.
       ccc(1) = sqrt(   m)
       bbb(1) = sqrt(1.-m)
       n = 0
       i = 1

       !! . . Arithmetic-Geometric Mean of A, B and C
       do while ( abs(ccc(i)) .gt. tol)    
          i = i + 1;
          aaa(i) =     (aaa(i-1) + bbb(i-1))/2.
          bbb(i) = sqrt(aaa(i-1) * bbb(i-1))
          ccc(i) =     (aaa(i-1) - bbb(i-1))/2.
          if (ccc(i) .ge. tol .or. ccc(i-1) .gt. tol ) then
             n = (i-1)
          end if
       end do

       phin = 0.
       phin = (2**dble(n))*aaa(i)*u
       do while (i .gt. 1)
          i = i - 1
          if ( n .ge. i ) then
             phin = (asin(ccc(i+1)*sin(phin)/aaa(i+1))+phin)/2.
          end if
       end do
       am = phin
       sn = sin(phin)
       cn = cos(phin)
       dn = sqrt(1.-m*sin(phin)**2)
    else if (m.eq.0.) then
       am = u
       sn = sin(u)
       cn = cos(u)
       dn = 1. 
    else if (m.eq.1.) then
       am = asin(tanh(u))
       sn = tanh(u)
       cn = 1./cosh(u) 
       dn = 1./cosh(u) 
    else 
       call seperr('M must be in the range 0 <= M <= 1.')
    end if

  end subroutine ellipj

  !!----------------------------------------------------------------

  subroutine deallocate_elliptic()
    deallocate ( aaa,bbb,ccc )
  end subroutine deallocate_elliptic

end module elliptic
