anchor.diff:  
   
          ! the basic idea here is that we do a characteristic
  -       ! decomposition.  The jump in primitive variables (Q)
  -       ! can be transformed to a jump in characteristic variables
  -       ! using the left and right eigenvectors.  Then each wave
  -       ! tells us how much of each characteristic quantity reaches
  -       ! the interface over dt/2.  We only add the quantity if
  -       ! it moves toward the interface.
  +       ! decomposition.  The jump in primitive variables (Q) can be
  +       ! transformed to a jump in characteristic variables using the
  +       ! left and right eigenvectors.  Then each wave tells us how
  +       ! much of each characteristic quantity reaches the interface
  +       ! over dt/2.  We only add the quantity if it moves toward the
  +       ! interface.
          !
  -       ! Given a jump dQ, we can express this as:
  +       ! See Miller & Colella for a good discussion of the
  +       ! characteristic form.  The basic form is:
          !
  -       !       3
  -       ! dQ = sum  (l_i dQ) r_i
  -       !      i=1
          !
  -       ! Then 
  -       !         3
  -       ! A dQ = sum  lambda_i (l_i dQ) r_i
  -       !        i=1
  -       !
  -       ! where lambda_i is the eigenvalue.  (See for example LeVeque's
  -       ! book).
  -       !
  +       !  n+1/2      ~               ~
  +       ! q         = q  -  sum l . ( q   - I  ) r
  +       !  i+1/2,L     L     i   i     L     +    i
  +       ! 
  +       !       ~
  +       ! Where q is the reference state.
  +
          
  -       ! define the reference states
  +       ! define the reference states -- Miller & Colella (2002) argue
  +       ! picking the fastest wave speed.  We follow the convention from
  +       ! Colella & Glaz, and Colella (1990) -- this is intended to 
  +       ! minimize the size of the jump that the projection operates on.
          if (.false.) then
             ! CASTRO method
             Qref_xm(:) = Q%data(i,:)
  
positive.code:     !
      !       / u   r   0   \  
      !   A = | 0   u   1/r |  
      !       \ 0  rc^2 u  /
      !
      ! With the rows corresponding to rho, u, and p
      !
  
negative.code:     !  q(xi) = qminus + xi*(qplus - qminux + q6 * (1-xi) )
      !
      ! with xi = (x - xl)/dx, where xl is the interface of the left
      ! edge of the cell.  qminus and qplus are the values of the 
      ! parabola on the left and right edges of the current cell.
  
      ! Limit the left and right values of the parabolic interpolant
      ! (C&W Eq. 1.10).  Here the loop is over cells, and considers the
      ! values on either side of the center of the cell (Qminus and
      ! Qplus).
      do n = 1, nprim
         do i = U%grid%lo-1, U%grid%hi+1
  
  