program oblique
c
c*********************************************************************c
c*                                                                   *c
c*  Written by Timothy J. Cowan                                      *c
c*                                                                   *c
c*********************************************************************c
c*
c*  Description: Interactive oblique shock calculator for a double 
c*               wedge.
c*
c*  References:  John, J.E.A., "Gas Dynamics," 2nd Edition, 
c*               1984: pp 123-134.
c*               Ames Research Staff, "Equations, Tables, and Charts
c*               for Compresssible Flow," NACA Report 1135, pp 8-12.
c*
c*  Updated:     9/25/00 - Fixed singularity problems for oblique
c*               shock case.  Solves for singularity and makes an
c*               educated initial guess for the shock angle.
c*
c*               9/26/00 - Added computation of downstream pressure 
c*               and density ratios.
c*
c*********************************************************************c
c*
c*  Variable definitions
c*
c*  pi    => constant value for pi
c*  gam   => constant value for ratio of specific heats
c*  gm1   => constant value, gam - 1.0
c*  gp1   => constant value, gam + 1.0
c*  tol   => convergence tolerance
c*
c*  xm1   => free stream Mach number
c*  delta => wedge half angle in degrees
c*  del   => wedge half angle in radians
c*
c*  q     => shock angle
c*  xm2   => Mach number behind the shock
c*  t12   => temperature ratio across the shock
c*
c*********************************************************************c
c
      implicit real*8 (A-H,O-Z)

      character*20 textread, filen
      real*8 obl_fnc, obl_drv
c
c---------------------------------------------------------------------c
c *** Define some constants :
c---------------------------------------------------------------------c
c
      pi = 4.0d0*datan(1.0d0)
      gam = 1.40d0
      gm1 = gam - 1.0d0
      gp1 = gam + 1.0d0
      tol = 0.000000001d0

      niter = 100
c
c---------------------------------------------------------------------c
c *** Problem definition:
c---------------------------------------------------------------------c
c
      write(*,'(5(/),'' *** Program Oblique Shock Wave ***'',3(/))')

      write(*,*) 
      write(*,'(/,A,$)') 'Oblique shock or expansion fan (1/0) : '
      read(*,*) itype

      write(*,*) 
      write(*,'(/,A,$)') 'Input the freestream Mach number : '
      read(*,*) xm1
      write(*,'(/,A,$)') 'Input the wedge half angle (deg) : '
      read(*,*) delta

      del = delta*pi/180.0d0
      xmm = xm1*xm1
c
c---------------------------------------------------------------------c
c *** Use iteration to solve for shock angle & downstream Mach number 
c---------------------------------------------------------------------c
c
      if (itype .EQ. 0) goto 200

c *** define singular points in model equation ( den=0, or q=90deg )
c
      qs1 = DASIN( DSQRT( 1.0d0/xmm ) )
      qs2 = 90.0d0*pi/180.0d0

c *** make an initial guess which is to the right of qs1
c
      q = qs1 + 0.01d0

c *** Use Newton's method to solve for root
c
      do i = 1,niter

c ***** update solution using Newton's method
c
        f = obl_fnc( q, del, gam, xm1 )
        df = obl_drv( q, del, gam, xm1 )

        q = q - f/df
        
c ***** make sure angle stays away from singularities
c
   10   if ( q .LE. qs1 ) then
          q = q + 0.5d0*f/df
          goto 10
        endif
   20   if ( q .GE. qs2 ) then
          q = q + 0.5d0*f/df
          goto 20
        endif

c ***** check if tolerance has been met
c
        if ( DABS( f ) .LE. tol ) goto 100

      enddo

      write(*,*) ' '
      write(*,*) ' '
      write(*,*) 'Warning!!!!!!'
      write(*,*) '>> Failed to converge within ',niter,' iterations.'

  100 continue
c
c---------------------------------------------------------------------c
c *** Output the results
c---------------------------------------------------------------------c
c
      sq = DSIN(q)
      cq = DCOS(q)
      sqs = sq*sq
      cqs = cq*cq

      xm2 = ( gp1*gp1*xmm*xmm*sqs - 4.0d0*(xmm*sqs - 1.0d0)*
     &      (gam*xmm*sqs + 1.0d0) )/( (2.0d0*gam*xmm*sqs - gm1)*
     &      (gm1*xmm*sqs + 2.0d0) )
      xm2 = DSQRT( xm2 )
      cp = 4.0d0*( xmm*sqs - 1.0d0 )/( gp1*xmm )

      r2  = gp1*xmm*sqs/( gm1*xmm*sqs + 2.0d0 )
      p2  = 1.0d0 + 2.0d0*gam*( xmm*sqs - 1.0d0 )/gp1
      t2  = p2/r2

      write(*,*) ' '
      write(*,*) ' '
      write(*,*) 'Given :'
      write(*,*) '  Freestream Mach  => ', xm1
      write(*,*) '  Wedge half angle => ', delta, 'degrees'
      write(*,*) ' '
      write(*,*) 'Computed :'
      write(*,*) '  Shock angle           => ', q*180.0d0/pi, 'degrees'
      write(*,*) '  Downstream Mach       => ', xm2
      write(*,*) '  Downstream Cp         => ', cp
      write(*,*) '  Press. Ratio (p2/p1)  => ', p2
      write(*,*) '  Temp. Ratio (T2/T1)   => ', t2
      write(*,*) '  Density Ratio (r2/r1) => ', r2

      goto 999
c
c---------------------------------------------------------------------c
c *** Prandtl-Meyer Expansion fan solution
c---------------------------------------------------------------------c
c
  200 continue

c *** Compute the Prandtl-Meyer angle, v1
c
      v1 = exp_fnc( 0.0d0, gam, xm1 )

c *** The total turning angle is v2 = v1 + 2*delta
c
      v2 = v1 + 2.0d0*del

c *** Guess a Mach number behind the expansion
c
      xm2 = xm1 + 10.0d0*xm1/(180.0d0*v2/pi)

c *** Use Newton's method to solve for root
c
      do i = 1,niter

        f = exp_fnc( v2, gam, xm2 ) 
        df = exp_drv( v2, gam, xm2 )
        xm2 = xm2 - f/df
          
c ***** check if tolerance has been met
c
        write(*,*) i, f, df
        if ( DABS( f ) .LE. tol ) goto 300

      enddo

      write(*,*) ' '
      write(*,*) ' '
      write(*,*) 'Warning!!!!!!'
      write(*,*) '>> Failed to converge within ',niter,' iterations.'

  300 continue
c
c---------------------------------------------------------------------c
c *** Output the results
c---------------------------------------------------------------------c
c
c *** Compute the fan angle & pressure coefficient
c
      xmu = 2.0d0*del + DATAN( 1.0d0/DSQRT(xm1*xm1 - 1.0d0) ) 
     &                - DATAN( 1.0d0/DSQRT(xm2*xm2 - 1.0d0) )

      t2  = ( 1.0d0 + 0.5d0*gm1*xmm )/( 1.0d0 + 0.5d0*gm1*xm2*xm2 )
      p2 = t2**(gam/gm1)
      r2  = p2/t2
      cp = 2.0d0*( p2 - 1.0d0 )/(gam*xmm)

      write(*,*) ' '
      write(*,*) ' '
      write(*,*) 'Given :'
      write(*,*) '  Freestream Mach    => ', xm1
      write(*,*) '  Wedge half angle   => ', delta, 'degrees'
      write(*,*) ' '
      write(*,*) 'Computed :'
      write(*,*) '  Prandtl-Meyer Ang. => ', v1*180.0d0/pi, 'degrees'
      write(*,*) '  Exp. fan angle     => ', xmu*180.0d0/pi, 'degrees'
      write(*,*) '  Downstream Mach    => ', xm2
      write(*,*) '  Downstream Cp      => ', cp
      write(*,*) '  Press. Ratio (p2/p1)  => ', p2
      write(*,*) '  Temp. Ratio (T2/T1)   => ', t2
      write(*,*) '  Density Ratio (r2/r1) => ', r2

c
c---------------------------------------------------------------------c
c *** Successful completion of program
c---------------------------------------------------------------------c
c
  999 write (*,*) ' '
      stop ' OK!!'
      end
c*********************************************************************c
c*********************************************************************c







      real*8 function obl_fnc( q, d, g, m )
c*********************************************************************c
c
      implicit none
      real*8 f, q, d, g, m
      real*8 sq, cq, tq, m2
c
c---------------------------------------------------------------------c
c
c *** define sin, cos, & tan values
c
      sq = DSIN(q)
      cq = DCOS(q)
      tq = sq/cq

      m2 = m*m

      obl_fnc = tq*( 0.5d0*(g + 1.0d0)*m2 / 
     &          (m2*sq*sq - 1.0d0) - 1.0d0 ) - 1.0d0/DTAN(d)

      return
c
      end
c*********************************************************************c
c*********************************************************************c






      real*8 function obl_drv( q, d, g, m )
c*********************************************************************c
c
      implicit none
      real*8 q, d, g, m
      real*8 sq, cq, m2, c1
c
c---------------------------------------------------------------------c
c
      sq = DSIN(q)
      cq = DCOS(q)

      m2 = m*m
      c1 = m*m*sq*sq - 1.0d0

      obl_drv = 1.0d0/(cq*cq) * ( 0.5d0*(g + 1.0d0)*m2 / c1 - 1.0d0 )
     &        - ( (g + 1.0d0)*m2*m2*sq*sq / ( c1*c1 ) )

      return
c
      end
c*********************************************************************c
c*********************************************************************c







      real*8 function exp_fnc( v, g, m )
c*********************************************************************c
c
      implicit none
      real*8 v, g, m, ggg, mmm
c
c---------------------------------------------------------------------c
c
      ggg = DSQRT( (g + 1.0d0)/(g - 1.0d0) )
      mmm = DSQRT( m*m - 1.0d0 )

      exp_fnc = ggg*DATAN( mmm/ggg ) - DATAN( mmm ) - v

      return
c
      end
c*********************************************************************c
c*********************************************************************c







      real*8 function exp_drv( v, g, m )
c*********************************************************************c
c
      implicit none
      real*8 v, g, m, ggg, mmm
c
c---------------------------------------------------------------------c
c
      ggg = DSQRT( (g + 1.0d0)/(g - 1.0d0) )
      mmm = DSQRT( m*m - 1.0d0 )

      exp_drv = m/( 1.0d0 + mmm*mmm/(ggg*ggg) )/mmm - 1.0d0/(m*mmm)

      return
c
      end
c*********************************************************************c
c*********************************************************************c







      character*80 function textread( prompt)
c*********************************************************************c
c
      character*(*) prompt
c
c---------------------------------------------------------------------c
c
      write(*,'(/,a,$)') prompt
      read(*,'(a)') textread
      return
c
      end
c*********************************************************************c
c*********************************************************************c






      integer function namlen( filen )
c*********************************************************************c
c
      character*20 filen
c
c---------------------------------------------------------------------c
c
      namlen = 0
      do i = 20,1,-1
        if ( filen(i:i) .ne. ' ' ) then
          namlen = i
          goto 101
        endif
      enddo
  101 return
c
      end
c*********************************************************************c
c*********************************************************************c