! PROGRAM MaximumRampAngleStudy
! PURPOSE - Make a test of subroutine MaxRampAngle
! AUTHOR -  Ralph L. Carmichael, Public Domain Aeronautical Software
! REVISION HISTORY
!   DATE  VERS PERSON  STATEMENT OF CHANGES
! 11Oct02  1.0   RLC   Initial release to the web site
! 25Jan03  1.1   RLC   Some minor restructuring

!+
MODULE ObliqueShockAux
! ---------------------------------------------------------------------------
! PURPOSE - This module holds the utility subroutine MaxRampAngle that 
!  computes the maximum ramp angle for an attached oblique shock for a given
!  value of Mach number and gamma.
IMPLICIT NONE
  INTEGER,PRIVATE,PARAMETER:: WP = SELECTED_REAL_KIND(10,50)
  REAL(WP),PRIVATE:: machRef,gammaRef   ! set by MaxRampAngle; used by Eq138
!----------------------------------------------------------------------------
  PRIVATE:: Eq138  ! not for use outside this module
  PUBLIC:: MaxRampAngle
CONTAINS
!+
FUNCTION Eq138(theta) RESULT(cotDelta)
! ---------------------------------------------------------------------------
! PURPOSE - Special version of Eq138 that picks up mach and gamma as module
!  variables. Used with BrentMin which requires a function of one variable.
  REAL(WP),INTENT(IN):: theta   ! wave angle, radians
  REAL(WP):: cotDelta   ! cotangent of deflection angle.
  
  REAL(WP):: gp1,ms
  REAL(WP),PARAMETER:: ONE=1.0_WP, TWO=2.0_WP
!----------------------------------------------------------------------------
  gp1=gammaRef+ONE                            ! gammaref is a module variable
  ms=(machRef*SIN(theta))**2                   ! machRef is a module variable
  IF (ms <= ONE) THEN
    cotDelta=1E38   ! something really huge 
  ELSE
    cotDelta=TAN(theta)*(gp1*machRef**2/(TWO*(ms-ONE))-ONE)
  END IF
  RETURN
END Function Eq138   ! ------------------------------------------------------

!+
SUBROUTINE MaxRampAngle(mach,gamma,maxRamp,thetaMaxRamp,neval)
! ---------------------------------------------------------------------------
! PURPOSE - Compute the maximum deflection for given values of mach and gamma.
!  Include the wave angle where this maximum occurs as an output and the
!  number of function evaluations used by BrentMin.
! NOTES - Eq138 returns cot(delta). We want a maximum of delta, so we find the
!  minimum of cot(delta).
!  The usage of BrentMin requires that the function to be minimized be a 
!  function of only one variable. Equation 138 expresses cot(delta) as a
!  function of mach, gamma, and theta. Since theta is our independent
!  variable, put mach and gamma into module variables that may be accessed by
!  all procedures in the module.

 USE ForsytheMalcolmMoler,ONLY: BrentMin
  REAL(WP),INTENT(IN):: mach   ! Mach number
  REAL(WP),INTENT(IN):: gamma  ! ratio of specific heats
  REAL(WP),INTENT(OUT):: maxRamp   ! max deflection angle, radians
  REAL(WP),INTENT(OUT):: thetaMaxRamp  ! theta where max occurs, radians
  INTEGER, INTENT(OUT):: neval   ! number of function evaluations in BrentMin

  REAL(WP),PARAMETER:: HALF=0.5_WP, ONE=1.0_WP
  INTEGER, PARAMETER:: MAXITER=20
  REAL(WP),PARAMETER:: PI = 3.141592653589793238462643_WP
  REAL(WP),PARAMETER:: TOL=1E-8
   
  REAL(WP):: ax,bx   ! limits on theta
  INTEGER:: errCode  ! receives error code from BrentMin
  REAL(WP):: fmin    ! minimum value of the objective function
  REAL(WP):: xmin    ! value of the independent variable where minimum occurs
!----------------------------------------------------------------------------
  machRef=mach        ! so Eq138 can access it 
  gammaRef=gamma      ! so Eq138 can access it
  ax=ASIN(ONE/mach)   ! smallest wave angle (Mach angle)
  bx=HALF*PI          ! wave angle cannot exceed pi/2
  CALL BrentMin(ax,bx,Eq138,tol,MAXITER,neval,errCode,xmin,fmin)
  IF (errCode /= 0) WRITE(*,*) 'Non-zero errCode from BrentMin'
  thetaMaxRamp=xmin
  maxRamp=ATAN(ONE/fmin)  ! fmin is cotangent of delta
  RETURN
END Subroutine MaxRampAngle   ! ---------------------------------------------


END Module ObliqueShockAux   ! ==============================================


!+
PROGRAM MaximumRampAngleStudy
! ---------------------------------------------------------------------------
! PURPOSE - Compute the maximum ramp angle for an attached oblique shock wave
!  at given Mach angle and gamma.
! AUTHOR  - Ralph L. Carmichael, Public Domain Aeronautical Software
! REVISION HISTORY
!   DATE   VERS PERSON  STATEMENT OF CHANGES
! 15Jul02   1.0   RLC   Extracted code from other programs
! 25Jan03   1.1   RLC   Some restructuring for clarity
USE ObliqueShockAux,ONLY: MaxRampAngle
IMPLICIT NONE
 
  INTEGER,PARAMETER:: WP = SELECTED_REAL_KIND(10,50)
  REAL(WP),PARAMETER:: PI = 3.141592653589793238462643_WP
  REAL(WP),PARAMETER:: GAMMA = 1.4_WP
  CHARACTER(LEN=*),PARAMETER:: FMT = '(F7.3,2F12.7,2F10.3,I6)'
  INTEGER:: k
  REAL(WP):: deltaMax
  REAL(WP),DIMENSION(15),PARAMETER:: machTable = (/ &
    1.05, 1.1, 1.2, 1.3, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0,4.0,5.0,6.0,8.0,10.0 /)
  INTEGER:: neval  
  REAL(WP):: thetaDeltaMax
! ---------------------------------------------------------------------------
! For each Mach number in machTable, use MaxRampAngle to compute deltaMax and
!  thetaDeltaMax. Print these angles in radians, then in degrees, followed by
!  the number of function evaluations used by BrentMin in MaxRampAngle.
  WRITE(*,*) '  mach   delta,rad   theta,rad  delta,deg theta,deg neval'
  DO k=1,SIZE(machTable)
    CALL MaxRampAngle(machTable(k),GAMMA,deltaMax,thetaDeltaMax,neval)
    WRITE(*,FMT) machTable(k),deltaMax,thetaDeltaMax,   &
    (180/PI)*deltaMax, (180/PI)*thetaDeltaMax, neval
  END DO
  STOP
END Program MaximumRampAngleStudy   ! =======================================