anchor.diff:        real(rs), parameter, private :: big_number = 10.e36_rs
   
   !  ** Hide the helper functions and subroutines
  -      private :: inverse, &
  +      private :: cart2incaz, &
  +                 cross_prod, &
  +                 incaz2cart, &
  +                 inverse, &
                    determinant
   
         CONTAINS
  +
  +!===============================================================================
  +   subroutine CIJ_phase_vels(C, az, inc, pol, avs, vp, vs1, vs2, vsmean, &
  +                             azp, incp, azs1, incs1, azs2, incs2)
  +!===============================================================================
  +!
  +!  Calculate phase velocities for a given direction in an elasticity tensor
  +!
  +!  INPUT:
  +!     C(6,6)       : Elasticity tensor, density-normalised (Aij) [m^2.s^-2]
  +!     az           : Azimuth from the +x1 axis towards the -x2 axis
  +!                    (When looking down the x3 axis towards the origin, azi is
  +!                    +ve clockwise.) [degrees]
  +!     inc          : Incidence measured from the x1-x2 axis towards x3 [degrees]
  +!
  +!  OUTPUT (all optional):
  +!     pol          : Orientation of fast shear wave.  This is measured away from
  +!                    the projection of the +ve x3 axis onto the plane normal to
  +!                    the direction of interest.  Positive is clockwise when
  +!                    looking DOWN the ray towards the origin.  (NB: This is in
  +!                    the opposite way to the ray frame convention.) [degrees]
  +!     avs          : Shear wave anisotropy, measured as 200*(Vs1-Vs2)/(Vs1+Vs2).
  +!     vp           : P-wave phase velocity. [m.s^-1]
  +!     vs1          : Shear wave velocity of fast shear wave. [m.s^-1]
  +!     vs2          : Shear wave velocity of slow shear wave. [m.s^-1]
  +!     vsmean       : Harmonic mean of the fast and slow shear wave velocities.
  +!                    [m.s^-1]
  +!     azp, incp    : Azimuth and inclination of P vibration [deg].
  +!     azs1, incs1  : Same, but for fast S wave [deg].
  +!     azs2, incs2  : Same, but for slow S wave [deg].
  +      real(rs), intent(in) :: C(6,6), az, inc
  +      real(rs), intent(out), optional :: pol, avs, vp, vs1, vs2, vsmean, &
  +                                         azp, incp, azs1, incs1, azs2, incs2
  +      real(rs) :: vp_in, vs1_in, vs2_in
  +      real(rs) :: x(3)  ! Unit vector of interest
  +      real(rs) :: T(3,3) ! Christoffel T matrix
  +      real(rs) :: xp(3), xs1(3), xs2(3) ! Vibration polarisation vectors
  +
  +      ! Get direction vector
  +      x = incaz2cart(inc, az)
  +      ! Construct the Christoffel matrix
  +      call make_T
  +      ! Find eigenvalues and -vectors, which yield the velocities and pols
  +      call T_eigs
  +      ! Fill in needed outputs
  +      ! Find polarisation of fast shear wave
  +      if (present(pol)) call get_pol
  +      if (present(vp)) vp = vp_in
  +      if (present(vs1)) vs1 = vs1_in
  +      if (present(vs2)) vs2 = vs2_in
  +      if (present(avs)) avs = 200._rs*(vs1_in - vs2_in)/(vs1_in + vs2_in)
  +      if (present(vsmean)) vsmean = 2._rs/(1._rs/vs1_in + 1._rs/vs2_in)
  +      if (present(azp)) then
  +         if (.not.present(incp)) then
  +            write(0,'(a)') 'anisotropy_ajn: CIJ_phase_vels: Both azp and incp must be present'
  +            stop 1
  +         endif
  +         call cart2incaz(xp, incp, azp)
  +      endif
  +      if (present(azs1)) then
  +         if (.not.present(incs1)) then
  +            write(0,'(a)') 'anisotropy_ajn: CIJ_phase_vels: Both azs1 and incs1 must be present'
  +            stop 1
  +         endif
  +         call cart2incaz(xs1, incs1, azs1)
  +      endif
  +      if (present(azs2)) then
  +         if (.not.present(incs2)) then
  +            write(0,'(a)') 'anisotropy_ajn: CIJ_phase_vels: Both azs2 and incs2 must be present'
  +            stop 1
  +         endif
  +         call cart2incaz(xs2, incs2, azs2)
  +      endif
  +
  +   contains
  +      subroutine make_T()
  +      ! Create the Christoffel T matrix
  +         integer, parameter :: ijkl(3,3) = reshape((/1, 6, 5,&
  +                                                     6, 2, 4,&
  +                                                     5, 4, 3/), (/3, 3/))
  +         integer :: i, j, k, l, m, n
  +         T = 0._rs
  +         do i = 1, 3
  +            do j = 1, 3
  +               do k = 1, 3
  +                  do l = 1, 3
  +                     m = ijkl(i,j)
  +                     n = ijkl(k,l)
  +                     T(i,k) = T(i,k) + C(m,n)*x(j)*x(l)
  +                  enddo
  +               enddo
  +            enddo
  +         enddo
  +      end subroutine make_T
  +
  +      subroutine T_eigs()
  +      ! Compute eigenvalues of T and sort into decreasing order, filling in
  +      ! the velocity and polarisation vectors for P, S1 and S2
  +         real(rs) :: eigval(3), eigvec(3,3)
  +         integer :: kp, ks1, ks2
  +         call eig_jacobi(T, 3, eigval, eigvec)
  +         kp = maxloc(eigval, 1)
  +         ks2 = minloc(eigval, 1)
  +         ks1 = 6 - kp - ks2
  +         if (ks1 < 1 .or. ks1 > 3) then
  +            write(0,'(a)') 'anisotropy_ajn: CIJ_phase_vels: Error: ks1 is not in range 1-3'
  +            stop
  +         endif
  +         vp_in  = sqrt(eigval(kp))
  +         vs1_in = sqrt(eigval(ks1))
  +         vs2_in = sqrt(eigval(ks2))
  +         xp  = eigvec(:,kp)
  +         xs1 = eigvec(:,ks1)
  +         xs2 = eigvec(:,ks2)
  +      end subroutine T_eigs
  +
  +      subroutine get_pol()
  +      ! Convert the S1 polarisation vector into a fast orientation
  +         real(rs), dimension(3) :: u, xs1p, v
  +         ! Project the fast vector onto the wavefront plane
  +         xs1p = cross_prod(x, cross_prod(x, xs1))
  +         xs1p = xs1p/sqrt(sum(xs1p**2))  ! Normalise
  +         ! Local up vector
  +         u = incaz2up(inc, az)
  +         ! Angle between the projected polarisation and the up vector is pol
  +         v = cross_prod(xs1p, u)
  +         ! Use atan2 to get the correct quadrant (sin/cos = |a^b|/a.b)
  +         pol = to_deg*atan2(sqrt(sum(v**2)), dot_product(xs1p, u))
  +         ! If v is codirectional with x, then the angle is correct; otherwise,
  +         ! we're measuring the angle the wrong way round
  +         if (dot_product(x, v) < 0._rs) pol = -pol
  +         pol = modulo(pol + 90._rs, 180._rs) - 90._rs
  +      end subroutine get_pol
  +
  +      function incaz2up(inci, azi) result(u)
  +         ! Return the local 'up' vector which is normal to x and points along
  +         ! the unit sphere towards z
  +         real(rs), intent(in) :: inci, azi
  +         real(rs) :: u(3), raz, rinc, sini
  +         raz = to_rad*azi
  +         rinc = to_rad*inci
  +         sini = sin(rinc)
  +         u(1) =  cos(raz+pi)*sini
  +         u(2) = -sin(raz+pi)*sini
  +         u(3) =  cos(rinc)
  +      end function incaz2up
  +   end subroutine CIJ_phase_vels
  +!-------------------------------------------------------------------------------
   
   !===============================================================================
      subroutine thom(vp,vs,rho,eps,gam,del,c)
  
positive.code:    end function determinant
  !-------------------------------------------------------------------------------
  
     end module anisotropy_ajn
  !=======================================================================================
  
  
negative.code: !  MODULE : anisotropy_ajn
  !
  !  (C) James Wookey, September 2005 - 2008
  !  School of Earth Sciences, University of Bristol
  !  j.wookey@bristol.ac.uk
  !
  !  (C) Andy Nowacki, October 2008 -
  