Subroutine jacobi_grid(ichanl)
Use Numeric_Kinds_Module
Use CommonInfo_Module
USE PiFactrs_Module
Use APH_Module
Use Jacobi_Module
Use ChiAng_Module
Use QuantumNumber_Module
Implicit None
Save
!=========================================================================================
!  Routine to transform from APH to Jacobi coordinates.
! Assignment Scheme: (s1,S1,thetacap1) <== (rho1,theta1,chi1)
!                    (s2,S2,thetacap2) <== (rho1,theta1,chi2)
! ichanl = 1 => A ; 2 => B ; 3 => C
!
! Need to add in assignment for different Lamda values if Jtot > 0.
!   This system repeats for each lamda.

!=========================================================================================
!    I N P U T
INTEGER,INTENT(IN) :: ichanl

!=========================================================================================
!    I N T E R N A L S
CHARACTER(LEN=100) :: filename
LOGICAL  :: first=.true., write_chanl=.true.
INTEGER  :: irho, itheta, ichi, ispt, i
REAL(dp) :: smallsmin, smallsmax, largesmin, largesmax, thetacapmin, thetacapmax
REAL(dp) :: rhoterm, sintheta, coschi, sinchi, cos2chi, sin2chi, sinctheta, acosarg
REAL(dp) :: numerator, denominator, tanarg

!=========================================================================================
!    F U N C T I O N S

INTEGER :: aphindex

!=========================================================================================
! Allocate coordinate arrays

IF (.not.ALLOCATED(smalls)) THEN
  ALLOCATE(smalls(ndim), larges(ndim), thetacap(ndim), betaqf(ndim))
ENDIF

!=========================================================================================

Do irho=1,nrho
     rhoterm=rho_val(irho)/Sqrt(2.d0) 
   
     Do itheta=1,ntheta
        sintheta=Sin(theta_val(itheta))

        Do ichi=1,nchi

          sin2chi=Sin(2.d0*(chi_val(ichi)-chif1(ichanl)))
          cos2chi=Cos(2.d0*(chi_val(ichi)-chif1(ichanl)))

          ispt=nchi*ntheta*(irho-1)+nchi*(itheta-1)+ichi
          smalls(ispt)=rhoterm*Sqrt(1.d0-sintheta*cos2chi)
          larges(ispt)=rhoterm*Sqrt(1.d0+sintheta*cos2chi)
          numerator=sintheta*sin2chi
          denominator=(Sqrt(1.d0-(sintheta*cos2chi)**2))
          acosarg=numerator/denominator

!         This is to overcome division by zero when calculating thetacap.

          If(denominator.eq.0.d0) Then
            acosarg=(numerator/Abs(numerator))
          End If

!         These are to overcome machine precision and make sure Abs(acosarg).le.1
        
          If(Abs(acosarg).lt.1.d0.and.1.d0-Abs(acosarg).lt.5.d-14) Then
            acosarg=(numerator/Abs(numerator))
          End If

          If(Abs(acosarg).gt.1.d0.and.Abs(acosarg)-1.d0.lt.5.d-14) Then
            acosarg=(numerator/Abs(numerator))
          End If

          If(Abs(acosarg).gt.1.d0.and.Abs(acosarg)-1.d0.gt.5.d-14) Then
              Print*,'ERROR: Argument of acos outside -1 to 1'
              Write(Out_Unit,*) '========================================'
              Write(Out_Unit,*) 'ERROR: Argument of acos outside -1 to 1' 
              Stop 'jacobi_grid'          
          End If

          thetacap(ispt)=acos(acosarg)

          ! Calculate Betaqf
          sinchi=Sin(chi_val(ichi)-chif1(ichanl))
          coschi=Cos(chi_val(ichi)-chif1(ichanl))
          sinctheta=Sin(thetacap(ispt))
         
          numerator=smalls(ispt)*sinchi*sinctheta
          denominator=(larges(ispt)*coschi) + (smalls(ispt)*sinchi*acosarg)
          tanarg=numerator/denominator
          betaqf(ispt)=ATan(tanarg)

          If(Abs(betaqf(ispt)).gt.pi) Then
            Print*,'ERROR: betaqf out-of-bounds'
            Print*,'betaqf =',betaqf
            STOP 'jacobi_grid'
          EndIf

          End Do
     End Do
End Do

!=========================================================================
! Find minimum and maximum values for each coordinate

IF (first) THEN

! Find boundary values along rho infinity

largesmin=1.d10
largesmax=0.d0
smallsmin=1.d10
smallsmax=0.d0
thetacapmin=1.d10
thetacapmax=0.d0

irho=ind_rho_infty
DO itheta=1,ntheta
  DO ichi=1,nchi

    i=aphindex(ntheta,nchi,irho,itheta,ichi,.false.)

    IF (larges(i).gt.largesmax) THEN
      largesmax=larges(i)
    ENDIF

    IF (larges(i).lt.largesmin) THEN
      largesmin=larges(i)
    ENDIF

    IF (smalls(i).gt.smallsmax) THEN
      smallsmax=smalls(i)
    ENDIF

   IF (smalls(i).lt.smallsmin) THEN
     smallsmin=smalls(i)
   ENDIF

   IF (thetacap(i).gt.thetacapmax) THEN
     thetacapmax=thetacap(i)
   ENDIF

   IF (thetacap(i).lt.thetacapmin) THEN
     thetacapmin=thetacap(i)
   ENDIF

  ENDDO
ENDDO

WRITE(Out_Unit,41) 'Jacobi Grid Info at rho_infty = ', rho_val(ind_rho_infty)
WRITE(Out_Unit,42) 'min, au', 'max, au'
WRITE(Out_Unit,43) 'large S  ',largesmin,largesmax
WRITE(Out_Unit,43) 'small s  ',smallsmin,smallsmax
WRITE(Out_Unit,43) 'cap theta',thetacapmin,thetacapmax

first=.false.

ENDIF

IF (write_chanl) THEN

  IF (ichanl.eq.req_chanls) write_chanl=.false.

largesmin=1.d10
largesmax=0.d0
smallsmin=1.d10
smallsmax=0.d0
thetacapmin=1.d10
thetacapmax=0.d0

DO i=1,ndim

 IF (larges(i).gt.largesmax) THEN
   largesmax=larges(i)
 ENDIF

 IF (larges(i).lt.largesmin) THEN
   largesmin=larges(i)
 ENDIF

 IF (smalls(i).gt.smallsmax) THEN
   smallsmax=smalls(i)
 ENDIF

 IF (smalls(i).lt.smallsmin) THEN
   smallsmin=smalls(i)
 ENDIF

 IF (thetacap(i).gt.thetacapmax) THEN
   thetacapmax=thetacap(i)
 ENDIF

 IF (thetacap(i).lt.thetacapmin) THEN
   thetacapmin=thetacap(i)
 ENDIF

ENDDO

WRITE(Out_Unit,40) 'Jacob Grid Info for Channel: ',ichanl
WRITE(Out_Unit,42) 'min, au', 'max, au'
WRITE(Out_Unit,43) 'large S  ', largesmin, largesmax
WRITE(Out_Unit,43) 'small s  ', smallsmin, smallsmax
WRITE(Out_Unit,43) 'cap theta', thetacapmin, thetacapmax

ENDIF

40 FORMAT(/1X,a,I1)
41 FORMAT(/1X,a,F7.3)
42 FORMAT(1X,T18,a,T28,a,/1X,50('-'))
43 FORMAT(1X,a,6('.'),F7.3,3X,F7.3)

Return
End Subroutine jacobi_grid