!+
Module PrandtlMeyerMethods
! ---------------------------------------------------------------------------
! PURPOSE - Compute the Prandtl-Meyer function and its inverse
! AUTHOR  - Ralph L. Carmichael, Public Domain Aeronautical Software
! REVISION HISTORY
!   DATE   VERS PERSON  STATEMENT OF CHANGES
! 03Mar95   0.1   RLC   Original coding
! 03Jun02   0.2   RLC   Total revision
! 24Mar03   0.3   RLC   Restructured into a module
! 28Mar03   0.4   RLC   Defend against invalid input values
! NOTES - Versions of each function are provided for single and double 
!  precision variables with a MODULE PROCEDURE declaration that allows the
!  user to use the same function name regardless of precision.
!  Error values are given if invalid input data is supplied.
IMPLICIT NONE
 
  
  CHARACTER(LEN=*),PARAMETER,PUBLIC:: PM_METHODS_VERSION   = '0.3 (24Mar03)'
  INTEGER,PARAMETER,PRIVATE:: SP = SELECTED_REAL_KIND(6)  ! 6 decimal digits
  INTEGER,PARAMETER,PRIVATE:: DP = SELECTED_REAL_KIND(12)  ! 12 decimal digits
  REAL(DP),PARAMETER,PRIVATE:: ZERO=0.0_DP, ONE=1.0_DP, SIX=6.0_DP
  REAL(DP),PARAMETER,PRIVATE:: TWOTHIRDS=2.0_DP/3.0_DP
  REAL(DP),PARAMETER,PRIVATE:: PI = 3.141592653589793238462643_DP
  REAL(DP),PARAMETER,PRIVATE:: NUMAX = 2.27685322704990680000_DP ! (pi/2)*(Sqrt(6)-1)
  
  
  PRIVATE:: EvaluatePolynomial
  PUBLIC:: HallInversePrandtlMeyer
  PUBLIC:: NewtonInversePrandtlMeyer
  PUBLIC:: PrandtlMeyer
  
  INTERFACE HallInversePrandtlMeyer
    MODULE PROCEDURE HallInverseSingle,HallInverseDouble
  END INTERFACE
  INTERFACE NewtonInversePrandtlMeyer
    MODULE PROCEDURE NewtonInverseSingle,NewtonInverseDouble
  END INTERFACE
  INTERFACE PrandtlMeyer
    MODULE PROCEDURE PrandtlMeyerSingle,PrandtlMeyerDouble
  END INTERFACE

  
! ---------------------------------------------------------------------------
      
CONTAINS

!+
FUNCTION EvaluatePolynomial(x,a) RESULT(f)                  ! a PURE function
! ---------------------------------------------------------------------------
! PURPOSE -Evaluate a polynomial with coeff a(1),a(2),... at X
! NOTES - f=a(1)+a(2)*x+a(3)*x**2+....+a(n)*x**(n-1)
!   Straightforward implementation of Horner's rule
  REAL(DP),INTENT(IN):: x   ! evaluation point
  REAL(DP),INTENT(IN),DIMENSION(:):: a   ! coefficients
  REAL(DP):: f   ! value of polynomial at x
  INTEGER:: k,n
!----------------------------------------------------------------------------
  n=SIZE(a)   ! order of polynomial is n-1, not n
  IF (n <= 0) THEN
    f=0.0
    RETURN
  END IF

  f=a(n)
  DO k=n-1,1,-1
    f=a(k)+x*f
  END DO

  RETURN
END FUNCTION EvaluatePolynomial   ! -----------------------------------------

!+
FUNCTION HallInverseDouble(nu) RESULT(mach)                 ! a PURE function
! ---------------------------------------------------------------------------
! PURPOSE - Compute the inverse Prandtl-Meyer function using the rational
!  polynomial approximation of I.M. Hall
  REAL(DP),INTENT(IN):: nu   ! Prandtl-Meyer angle, radians
  REAL(DP):: mach   ! Mach number
  
! constants from I.M. Hall, The Aeronautical Journal, Sept 1975, p.417 
  REAL(DP),PARAMETER,DIMENSION(4):: C1 = &
    (/ 1.0_DP, 1.3604_DP, 0.0962_DP, -0.5127_DP /)
  REAL(DP),PARAMETER,DIMENSION(3):: C2 = (/ 1.0_DP, -0.6722_DP, -0.3278_DP /)
 
  REAL(DP):: y
!----------------------------------------------------------------------------
  IF (nu <= ZERO) THEN
    mach=ONE
    RETURN
  END IF
  IF (nu >= NUMAX) THEN
    mach=HUGE(mach)
    RETURN
  END IF
      
  y=(nu/NUMAX)**TWOTHIRDS
  mach=EvaluatePolynomial(y,C1)/EvaluatePolynomial(y,C2)
  RETURN
END Function HallInverseDouble   !-------------------------------------------

!+
FUNCTION HallInverseSingle(nu) RESULT(mach)                 ! a PURE function
! ---------------------------------------------------------------------------
! PURPOSE -Same as HallInverseDouble but for single precision variables.
!  When machDouble is converted to mach, some precision will be lost.
  REAL(SP),INTENT(IN):: nu
  REAL(SP):: mach
  REAL(DP):: nuDouble,machDouble
! ---------------------------------------------------------------------------
  nuDouble=nu
  machDouble=HallInverseDouble(nuDouble)
  mach=machDouble   ! converts DP to SP
  RETURN
END Function HallInverseSingle   ! ------------------------------------------
  

!+
FUNCTION NewtonInverseDouble(nu) RESULT(mach)               ! a PURE function
! ---------------------------------------------------------------------------
! PURPOSE - Compute the inverse Prandtl-Meyer function to the maximum of
!  double precision. I have done extensive experiments that show that full
!  precision will be obtained in four iterations of Newton's method when you
!  use the excellent approximation of Hall as a starting point.
  REAL(DP),INTENT(IN):: nu   ! Prandtl-Meyer angle, radians
  REAL(DP):: mach
  
  REAL(DP):: beta
  INTEGER:: k
  REAL(DP):: lambda,lambda1
  REAL(DP):: x,f,fp
!----------------------------------------------------------------------------
  IF (nu <= ZERO) THEN
    mach=ONE
    RETURN
  END IF  
  IF (nu >= NUMAX) THEN
    mach=HUGE(mach)
    RETURN
  END IF  

  lambda1=SQRT(SIX)
  lambda=ONE/lambda1
 
! First, we make an excellent guess using Hall
  x=HallInverseDouble(nu)
! Next, we make 4 iterative improvements, using Newton
  DO k=1,4
    beta=SQRT(x*x-ONE)
    f=lambda1*ATAN(lambda*beta)-ATAN(beta)  
    fp=(ONE-lambda**2)*beta/(x*(ONE+(lambda*beta)**2))
    x=x-(f-nu)/fp
  END DO  
  mach=x
  
  RETURN  
END Function NewtonInverseDouble   ! ----------------------------------------

!+
FUNCTION NewtonInverseSingle(nu) RESULT(mach)               ! a PURE function
! ---------------------------------------------------------------------------
! PURPOSE - Same as NewtonInverseDouble but for single precision variables.
!  When machDouble is converted to mach, some precision will be lost. 
  REAL(SP),INTENT(IN):: nu
  REAL(SP):: mach
  REAL(DP):: nuDouble,machDouble
! ---------------------------------------------------------------------------
  nuDouble=nu
  machDouble=NewtonInverseDouble(nuDouble)
  mach=machDouble   ! converts DP to SP
  RETURN
END Function NewtonInverseSingle   ! ----------------------------------------
  

!+
FUNCTION PrandtlMeyerDouble(mach) RESULT(nu)                ! a PURE function
! ---------------------------------------------------------------------------
! PURPOSE - Compute the value of the Prandtl-Meyer function
  REAL(DP),INTENT(IN):: mach   ! input Mach number
  REAL(DP):: nu   ! PM function, radians
  REAL(DP):: beta
  REAL(DP):: lambda,lambda1
!----------------------------------------------------------------------------
  IF (mach <= ONE) THEN
    nu=ZERO
    RETURN
  END IF  
  lambda1=SQRT(SIX)
  lambda=ONE/lambda1 
  beta=SQRT(mach*mach-ONE)
  nu=lambda1*ATAN(lambda*beta)-ATAN(beta)  
  RETURN
END Function PrandtlMeyerDouble   ! -----------------------------------------

!+
FUNCTION PrandtlMeyerSingle(mach) RESULT(nu)                ! a PURE function
! ---------------------------------------------------------------------------
  REAL(SP),INTENT(IN):: mach
  REAL(SP):: nu
  REAL(DP):: nuDouble,machDouble
! ---------------------------------------------------------------------------
  machDouble=mach
  nuDouble=PrandtlMeyerDouble(machDouble)
  nu=nuDouble   ! converts DP to SP
  RETURN
END Function PrandtlMeyerSingle   ! -----------------------------------------
  

END Module PrandtlMeyerMethods   ! ==========================================


! If you compile and run the following test program, you should get the
! following results:

!PRANDTL-MEYER
!   M     nu,deg
!  1.0    0.000000000000
!  2.0   26.379760813416
!  3.0   49.757346744346
!  4.0   65.784819797749
!  5.0   76.920215508539
!  6.0   84.955498177708
!  7.0   90.972732328571
!  8.0   95.624671710244
!  9.0   99.318098653168
! 10.0  102.316253173200
! 11.0  104.795679409621
! 12.0  106.878626078334
! 13.0  108.652159934792
! 14.0  110.179837794091
! 15.0  111.509068612111
! 16.0  112.675895924527
! 17.0  113.708188164695
! 18.0  114.627818327055
! 19.0  115.452185472551
! 20.0  116.195297574933
!INVERSE PRANDTL-MEYER
! nu    M,Hall     M,Newton
!  0   1.000000   1.000000000000000
! 10   1.434988   1.434974500874799
! 20   1.775064   1.774975810099576
! 30   2.134097   2.133905033231848
! 40   2.538077   2.537815452375472
! 50   3.012833   3.012608304504067
! 60   3.594034   3.594038276932490
! 70   4.338480   4.338994794662508
! 80   5.346451   5.347855733027846
! 90   6.816257   6.819035966666232
!100   9.205737   9.210489400662052
!110  13.867147  13.874600564266734
!120  27.325537  27.336595585186844
!130 630.877062 630.901085933426320

! It will be instructive to change the constant in the declaration of WP
! below from 12 to 6 and see the effect on the outcome.

!+
PROGRAM TestPrandtlMeyer
! ---------------------------------------------------------------------------
! PURPOSE - Test the operation of the functions in module PrandtlMeyerMethods
! AUTHOR  - Ralph L. Carmichael, Public Domain Aeronautical Software
! REVISION HISTORY
!   DATE   VERS PERSON  STATEMENT OF CHANGES
! 24Mar03   1.0   RLC   Original coding
!----------------------------------------------------------------------------
USE PrandtlMeyerMethods
IMPLICIT NONE

  INTEGER,PARAMETER:: WP = SELECTED_REAL_KIND(12)   ! double
  REAL(WP),PARAMETER:: PI = 3.141592653589793238462643_WP
  
  INTEGER:: k
  REAL(WP):: mach
  REAL(WP):: mhall,mnewton
  REAL(WP):: nu

!----------------------------------------------------------------------------
! First print the value of the Prandtl-Meyer function for Mach=1,2,..20
  WRITE(*,*) 'PRANDTL-MEYER'
  WRITE(*,*) '   M     nu,deg'
  DO k=1,20
    mach=REAL(k)
    nu=PrandtlMeyer(mach)
    WRITE(*,'(F6.1,F18.12)') mach,180*nu/PI
  END DO  

! Then print the inverse Prandtl-Meyer function from the Hall approximation
!  and from the converged Newton calculations
  WRITE(*,*) 'INVERSE PRANDTL-MEYER'
  WRITE(*,*) ' nu    M,Hall     M,Newton'
  DO k=0,130,10
    nu=(PI/180)*REAL(k)
    mhall=HallInversePrandtlMeyer(nu)
    mnewton=NewtonInversePrandtlMeyer(nu)
    WRITE(*,'(I4,F11.6,F20.15)') k,mhall,mnewton
  END DO
  
  STOP
END Program TestPrandtlMeyer   !=============================================