cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c...File oblshk.f                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      program oblshk

c        Iterates the oblique shock relations to determine the 
c        properties behind an oblique shock. Assumptions: perfect gas.
c        Written by John W. Slater, NASA Lewis Research Center.
c
c-----------------------------------------------------------------------

      Implicit none

      Real mach1, delta, gamma, tol
      Real mach2, theta, p2p1, t2t1, r2r1
      Real pi, dtr

c-----------------------------------------------------------------------

c...Header. 

      write(*,*) ' '
      write(*,*) '--- Oblshk ---'
                       
      write (*,*) 'Enter Mach number, wedge angle, gamma, tol'
      read (*,*) mach1, delta, gamma, tol

      pi    = 4.0 * atan(1.0)
      dtr   = pi / 180.
      delta = delta * dtr

      call sckobl ( mach1, delta, gamma, tol,  
     &              mach2, theta, p2p1, t2t1, r2r1 )

      write (*,*) '    '
      write (*,*) 'Oblique Shock Conditions'
      write (*,*) 'Inflow Mach number = ', mach1
      write (*,*) 'Wedge angle (deg)  = ', delta / dtr
      write (*,*) '     '
      write (*,*) 'Perfect gas properties behind shock'
      write (*,*) 'Outflow Mach       =  ', mach2
      write (*,*) 'Shock Angle (deg)  =  ', theta / dtr
      write (*,*) 'Pressure ratio     =  ', p2p1
      write (*,*) 'Temperature ratio  =  ', t2t1
      write (*,*) 'Density ratio      =  ', r2r1

c...End Statement.

      write (*,*) '     '
      write (*,*) '--- end of Oblshk ---'
      write (*,*) '     '

c-----------------------------------------------------------------------

      stop
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine sckobl ( mach1, delta, gamma, tol,  
     &                    mach2, theta, p2p1, t2t1, r2r1 )
c
c        Iterates the oblique shock relations to determine properties 
c        behind an oblique shock.  Assumptions: perfect gas.
c
c-----------------------------------------------------------------------

      Implicit none

      Real mach1, delta, gamma, tol
      Real mach2, theta, p2p1, t2t1, r2r1
      Real tbm
      Real delt1, delt2, thet1, thet2, thet
      Real rmn1, rmn2, gp1, gm1, t1, t2

c-----------------------------------------------------------------------

c...Choose an initial shock angle to be increment of wedge.        

      thet1 = 1.2 * delta
      delt1 = tbm ( mach1, thet1 )
      thet2 = 1.02 * thet1

c...Start of iteration.

   20 continue

        delt2 = tbm ( mach1, thet2, gamma )

c...    use secant iteration

        if ( abs((delt2-tan(delta))/tan(delta)) .gt. tol ) then
          thet  = thet2 - (tan(delta)-delt2)*(thet2-thet1)/(delt1-delt2)
          thet1 = thet2
          delt1 = delt2
          thet2 = thet
          go to 20
        endif

c...End of iterations.
                                                      
      theta = thet2

c...Compute properties behind shock based on perfect gas.               .

      rmn1  = mach1 * sin( theta )
              
      gp1   = gamma + 1.0
      gm1   = gamma - 1.0                            

      r2r1  = gp1 * rmn1**2 / ( gm1 * rmn1**2 + 2.0 )

      p2p1  = 1.0 + 2.0 * gamma * ( rmn1**2 - 1.0 ) / gp1

      t1    = ( rmn1**2 + 2.0 / gm1 ) 
      t2    = ( 2.0 * gamma * rmn1**2 / gm1 - 1.0 )
      rmn2  = sqrt( t1 / t2 )

      t2t1  = p2p1 / r2r1

      mach2 = rmn2 / sin( theta - delta )

c-----------------------------------------------------------------------

      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      function tbm ( mach, beta, gamma )

c        Evaluates theta-beta-mach relation for an oblique shock.
c
c-----------------------------------------------------------------------              

      Implicit none

      Real  mach, beta, gamma, c1, c2, tbm

c-----------------------------------------------------------------------              

      c1  = mach * mach * sin(beta) * sin(beta) - 1.0
      c2  = mach * mach * ( gamma + cos(2.0*beta) ) + 2.0
      tbm = 2.0 * (1.0/tan(beta)) * ( c1 / c2 )

c-----------------------------------------------------------------------

      return
      end