
c
c
c     =====================================================
      subroutine rp1(maxm,meqn,mwaves,mbc,mx,ql,qr,auxl,auxr,
     &                  wave,s,amdq,apdq)
c     =====================================================
c
c     # begin_html [use: jsMath]
c       [www.clawpack.org CLAWPACK] Riemann solver
c       in 1 space dimension.
c
c       Solve Riemann problems for the 1D acoustics equations with
c       variable coefficients (heterogeneous media)
c          \[   q_t + A(x) q_x = 0  \]
c       where
c          \[   q(x,t) = \vector{ p(x,t)\\ u(x,t)}           \]
c       and the coefficient matrix is          
c          \[   A(x) = \begin{matrix}
c                        0         & K(x)\\
c                        1/\rho(x) & 0 
c                        \end{matrix}.
c          \]
c
c       Here $p$ is the pressure  perturbation and $u$ the velocity.
c       The parameters $K$ (bulk modulus) and $\rho$ (density) are
c       used in the setaux routine to set the impedance and sound speed
c       in each cell.
c          
c       <b>On input:</b>
c         <ul>
c         <li> ql contains the state vector at the left edge of each cell,
c         <li> qr contains the state vector at the right edge of each cell,
c         <li> auxl(i,1) should contain the impedance Z in cell i:
c                 $Z = \sqrt{K\rho}$. 
c         <li> auxl(i,2) should contain the sound speed c in cell i:
c                 $c = \sqrt{K/\rho}$.
c         </ul> 
c
c       Note that the i'th Riemann problem has left state qr(i-1,:)
c                                          and right state ql(i,:).
c       From the basic clawpack routine step1, rp1 is called with ql = qr = q
c            and auxl = auxr = aux.
c
c       <b>On output:</b>
c         <ul>
c         <li> wave contains the waves, 
c         <li> s the speeds, 
c         <li> amdq the  left-going flux difference  ${\cal A}^-\Delta Q$, 
c         <li> apdq the right-going flux difference  ${\cal A}^+\Delta Q$ 
c       </ul>
c
c       For additional documentation on Riemann solvers rp1, see
c       [www.clawpack.org/rp1.html]
c
c       For details on solution of the Riemann problem for 
c       variable coefficient acoustics, see Chapter 9 of
c       [www.clawpack.org/book.html FVMHP].
c
c     # end_html
c
c
      implicit double precision (a-h,o-z)
c
      dimension auxl(1-mbc:maxm+mbc, 2)
      dimension auxr(1-mbc:maxm+mbc, 2)
      dimension wave(1-mbc:maxm+mbc, meqn, mwaves)
      dimension    s(1-mbc:maxm+mbc, mwaves)
      dimension   ql(1-mbc:maxm+mbc, meqn)
      dimension   qr(1-mbc:maxm+mbc, meqn)
      dimension apdq(1-mbc:maxm+mbc, meqn)
      dimension amdq(1-mbc:maxm+mbc, meqn)
      common /comlim/ mylim,mrplim(2)
c
c
c
c
c     # begin_html
c       Split the jump in $Q$ at each interface into waves.
c       First find $\alpha^1$ and $\alpha^2$, the coefficients of the 
c       2 eigenvectors.
c       
c       For the left-going wave we use the eigenvector for a 1-wave
c       propagating with speed $-c_{i-1}$ into the cell on the left, which is
c         \[ r^1_{i-1/2} = \vector{ -Z_{i-1} \\ 1} \] 
c
c       For the right-going wave we use the eigenvector for a 2-wave
c       propagating with speed $+c_i$ into the cell on the right, which is
c         \[ r^2_{i-1/2} = \vector{ Z_i \\ 1} \] 
c
c       Note that this corresponds to solving the Riemann problem
c         \[  q_t + A_{i-1/2} q_x = 0  \]
c       where the eigendecomposition of $A_{i-1/2}$ is 
c       $A_{i-1/2} = R_{i-1/2} \Lambda_{i-1/2} R^{-1}_{i-1/2}$, with
c          \[   R_{i-1/2} = \begin{matrix}
c                        -Z_{i-1} & Z_i \\
c                         1 & 1
c                        \end{matrix}, \quad
c               \Lambda_{i-1/2} = \begin{matrix}
c                        -c_{i-1} & 0 \\
c                         0 & c_i
c                        \end{matrix}, 
c          \]
c          \[
c               R^{-1}_{i-1/2} = \frac{1}{Z_{i-1}+Z_i}
c                        \begin{matrix}
c                        -1 & Z_i \\
c                         1 & Z_{i-1} 
c                        \end{matrix}, \quad
c          \]
c     # end_html
c
      do 20 i = 2-mbc, mx+mbc

c        # loop over all cells, including ghost cells

c        # impedances:
         zi = auxl(i,1)
         zim = auxl(i-1,1)

c        # begin_html
c        # Set $\delta = Q_i - Q_{i-1}$
c        # end_html
         delta1 = ql(i,1) - qr(i-1,1)
         delta2 = ql(i,2) - qr(i-1,2)

c        # begin_html
c        # Set $\alpha = R_{i-1/2}^{-1} \delta$
c        # end_html
         a1 = (-delta1 + zi*delta2) / (zim + zi)
         a2 =  (delta1 + zim*delta2) / (zim + zi)
c
c        # begin_html
c          Compute the waves and speeds:
c           \[ {\cal W}^1 = \alpha^1 r_{i-1/2}^1,~~~~
c                               {\cal W}^2 = \alpha^2 r_{i-1/2}^2,
c           \]
c           \[ s^1 = -c_{i-1},~~~~ s^2 = +c_i.  \]
c        # end_html
c
         wave(i,1,1) = -a1*zim
         wave(i,2,1) = a1
         s(i,1) = -auxl(i-1,2)
c
         wave(i,1,2) = a2*zi
         wave(i,2,2) = a2
         s(i,2) = auxl(i,2)
c
   20    continue
c
c
c     # begin_html
c       Compute the leftgoing and rightgoing fluctuations:
c
c       For this problem we always have $s^1 <0$ and $s^2 > 0$, so
c       \[
c          {\cal A}^-\Delta Q = s^1{\cal W}^1, \quad 
c          {\cal A}^+\Delta Q = s^2{\cal W}^2.
c       \]
c     # end_html
c
      do 220 m=1,meqn
         do 220 i = 2-mbc, mx+mbc
            amdq(i,m) = s(i,1)*wave(i,m,1)
            apdq(i,m) = s(i,2)*wave(i,m,2)
  220       continue
c
      return
      end
