Module m_Orthogonality

USE m_Parameters

USE m_Arrays, Only: ort, ans, bns, cns, uh, vh, zh, wg

IMPLICIT NONE

PRIVATE

PUBLIC :: Orthogonality_Symmetric, Orthogonality_Asymmetric, &
          Orthogonality_Test, Print_Coefficients


CONTAINS


SUBROUTINE Orthogonality_Symmetric (s)
 
  ! Computes the Orthogonality Characteristics
  ! for the Symmetric Hough Functions
  ! Using the Coefficients of Each Mode.

  IMPLICIT NONE

  INTEGER (KIND=ki8), INTENT (IN) :: s

  INTEGER (KIND=ki8) :: m, l, n

  REAL (KIND=kr8) :: akea, akeb, ape, tol, ek, rm, re

  tol=2.0**(-29_kr8)
  DO m=1,no
    DO l=1,no
      akea=0.0
      akeb=0.0
      ape=0.0
      DO n=1,na
        rm=REAL(s+2*n-2,kr8)
        re=rm*(rm+1.0)/epsl
        akea=akea+re*ans(n,m)*ans(n,l)
        re=(rm+1.0_kr8)*(rm+2.0_kr8)/epsl
        akeb=akeb+re*bns(n,m)*bns(n,l)
        ape=ape+cns(n,m)*cns(n,l)
      END DO
      ek=akea+akeb+ape
      IF (ABS(ek) < tol) ek=0.0_kr8
      IF (ABS(ek-1.0_kr8) < tol) ek=1.0_kr8
      ort(m,l)=ek
    END DO
  END DO

END SUBROUTINE Orthogonality_Symmetric


SUBROUTINE Orthogonality_Asymmetric (s)
 
  ! Computes the Orthogonality Characteristics
  ! for the Asymmetric Hough Functions
  ! Using the Coefficients of Each Mode.

  IMPLICIT NONE

  INTEGER (KIND=ki8), INTENT (IN) :: s

  INTEGER (KIND=ki8) :: m, l, n

  REAL (KIND=kr8) :: akea, akeb, ape, tol, ek, rm, re

  tol=2.0**(-29_kr8)
  DO m=1,no
    DO l=1,no
      akea=0.0
      akeb=0.0
      ape=0.0
      DO n=1,na
        rm=REAL(s+2*n-2,kr8)
        re=(rm+1.0_kr8)*(rm+2.0_kr8)/epsl
        akea=akea+re*ans(n,m)*ans(n,l)
        re=rm*(rm+1.0_kr8)/epsl
        akeb=akeb+re*bns(n,m)*bns(n,l)
        ape=ape+cns(n,m)*cns(n,l)
      END DO
      ek=akea+akeb+ape
      IF (ABS(ek) < tol) ek=0.0_kr8
      IF (ABS(ek-1.0_kr8) < tol) ek=1.0_kr8
      ort(m,l)=ek
    END DO
  END DO

END SUBROUTINE Orthogonality_Asymmetric


SUBROUTINE Orthogonality_Test (s, nome)
 
  ! Performs Test for Hough Function Orthogonality.
  ! This is Essentially a Convergence Test.

  IMPLICIT NONE

  INTEGER (KIND=ki8), INTENT (IN) :: s

  CHARACTER (LEN=6), INTENT (IN) :: nome

  INTEGER (KIND=ki8) :: j, ni, nj, n1, n2, n3

  REAL (KIND=kr8) :: tol, orth

  CHARACTER (LEN= 5) :: orn
  CHARACTER (LEN= 8) :: orm
  CHARACTER (LEN=11) :: fmc
  CHARACTER (LEN=19) :: fms

  CHARACTER (LEN=5), DIMENSION (no,no) :: orc

  IF (ipr < 0) RETURN
  n1=nb+1
  n2=2*nb
  n3=n2+1
  tol=2.0_kr8**(-24)
  fms='(/,A15,/,12(1X,I5))'
  fmc='(12(1X,A5))'
  DO ni=1,no
    DO nj=1,no
      orth=0.0_kr8
      DO j=1,jg
        orth=orth+(uh(ni,j)*uh(nj,j)+vh(ni,j)*vh(nj,j)+zh(ni,j)*zh(nj,j))*wg(j)
      END DO
      orth=2.0_kr8*orth
      IF (ABS(orth) <= tol) orth=0.0_kr8
      WRITE (UNIT=orm, FMT='(1PE8.1)') orth
      orn=orm(2:4)//orm(6:6)//orm(8:8)
      IF (orn == '0.0+0') orn=' Zero'
      IF (orn == '1.0+0') orn='  One'
      orc(ni,nj)=orn
    END DO
  END DO

  WRITE (UNIT=86, FMT='(/,1X,A6,1X,A11,I3,A15)') nome, ' orth  s = ', s, ' - Numerical  :'
  WRITE (UNIT=86, FMT=fms) ' gw x gw - n : ', (ni,ni=1,nb)
  DO ni=1,nb
    WRITE (UNIT=86, FMT=fmc) (orc(ni,nj),nj=1,nb)
  END DO
  WRITE (UNIT=86, FMT=fms) ' gw x ge - n : ', (ni,ni=n1,n2)
  DO ni=1,nb
    WRITE (UNIT=86, FMT=fmc) (orc(ni,nj),nj=n1,n2)
  END DO
  WRITE (UNIT=86, FMT=fms) ' gw x rb - n : ', (ni,ni=n3,no)
  DO ni=1,nb
    WRITE (UNIT=86, FMT=fmc) (orc(ni,nj),nj=n3,no)
  END DO
  WRITE (UNIT=86, FMT=fms) ' ge x ge - n : ', (ni,ni=n1,n2)
  DO ni=n1,n2
    WRITE (UNIT=86, FMT=fmc) (orc(ni,nj),nj=n1,n2)
  END DO
  WRITE (UNIT=86, FMT=fms) ' ge x rb - n : ', (ni,ni=n3,no)
  DO ni=n1,n2
    WRITE (UNIT=86, FMT=fmc) (orc(ni,nj),nj=n3,no)
  END DO
  WRITE (UNIT=86, FMT=fms) ' rb x rb - n : ', (ni,ni=n3,no)
  DO ni=n3,no
    WRITE (UNIT=86, FMT=fmc) (orc(ni,nj),nj=n3,no)
  END DO

  DO ni=1,no
    DO nj=1,no
      orth=ort(ni,nj)
      IF (ABS(orth) <= tol) orth=0.0_kr8
      WRITE (UNIT=orm, FMT='(1PE8.1)')orth
      orn=orm(2:4)//orm(6:6)//orm(8:8)
      IF (orn == '0.0+0') orn=' Zero'
      IF (orn == '1.0+0') orn='  One'
      orc(ni,nj)=orn
    END DO
  END DO

  WRITE (UNIT=87, FMT='(/,1X,A6,1X,A11,I3,A15)') nome, ' orth  s = ', s, ' - Analytical :'
  WRITE (UNIT=87, FMT=fms) ' gw x gw - n : ', (ni,ni=1,nb)
  DO ni=1,nb
    WRITE (UNIT=87, FMT=fmc) (orc(ni,nj),nj=1,nb)
  END DO
  WRITE (UNIT=87, FMT=fms) ' gw x ge - n : ', (ni,ni=n1,n2)
  DO ni=1,nb
    WRITE (UNIT=87, FMT=fmc) (orc(ni,nj),nj=n1,n2)
  END DO
  WRITE (UNIT=87, FMT=fms) ' gw x rb - n : ', (ni,ni=n3,no)
  DO ni=1,nb
    WRITE (UNIT=87, FMT=fmc) (orc(ni,nj),nj=n3,no)
  END DO
  WRITE (UNIT=87, FMT=fms) ' ge x ge - n : ', (ni,ni=n1,n2)
  DO ni=n1,n2
    WRITE (UNIT=87, FMT=fmc) (orc(ni,nj),nj=n1,n2)
  END DO
  WRITE (UNIT=87, FMT=fms) ' ge x rb - n : ', (ni,ni=n3,no)
  DO ni=n1,n2
    WRITE(UNIT=87, FMT=fmc) (orc(ni,nj),nj=n3,no)
  END DO
  WRITE (UNIT=87, FMT=fms) ' rb x rb - n : ', (ni,ni=n3,no)
  DO ni=n3,no
    WRITE (UNIT=87, FMT=fmc) (orc(ni,nj),nj=n3,no)
  END DO

END SUBROUTINE Orthogonality_Test


SUBROUTINE Print_Coefficients (cft, s, name6, name3)
 
  ! Printout the Coefficients Matrix CFT

  IMPLICIT NONE

  INTEGER (KIND=ki8), INTENT (IN) :: s

  CHARACTER (LEN=3), INTENT (IN) :: name3
  CHARACTER (LEN=6), INTENT (IN) :: name6

  REAL (KIND=kr8), DIMENSION (nm,no), INTENT (IN) :: cft

  INTEGER (KIND=ki8) :: n, m, mx, mm, mn

  REAL (KIND=kr8) :: cfx

  IF (ipr <= 0) RETURN
  WRITE (UNIT=88, FMT='(A)') ' '
  WRITE (UNIT=88, FMT='(1X,A,1X,2A,I3)') name3, name6, ' s = ', s
  WRITE (UNIT=88, FMT='(A)') ' '
  mm=0
  DO n=1,no
    mn=0
    cfx=-1.0_kr8
    mx=0
    DO m=1,na
      cfx=MAX(cfx,ABS(cft(m,n)))
      IF (cfx == ABS(cft(m,n))) mx=m
    END DO
    IF (cfx /= 0.0_kr8 .AND. mx /= 0) THEN
      DO m=mx,na
        IF (ABS(cft(m,n))/cfx < cutcft) THEN
          mn=m
          EXIT
        END IF
      END DO
    END IF
    mm=MAX(mm,mn)
    WRITE (UNIT=88, FMT='(3(A,I3),A,1PG12.5)')  &
                        '  n = ', n, ' mn = ', mn, ' mx = ', mx, ' cfx = ',cfx
    IF (ipr > 1) WRITE (UNIT=88, FMT='(1P6G12.5)') (cft(m,n),m=1,na)
  END DO
  WRITE (UNIT=88, FMT='(A,I3)') ' mm = ', mm

END SUBROUTINE Print_Coefficients


END Module m_Orthogonality