SUBROUTINE ComputeDeltaPQ( ray, c, gradc, cnn, cmn, cmm, d_phi, d_p_tilde, d_p_hat, d_q_tilde, d_q_hat )
TYPE( ray3DPt ), INTENT( IN ) :: ray
REAL( KIND=8 ), INTENT( IN ) :: c, gradc( 3 ), cnn, cmn, cmm
COMPLEX( KIND=8 ), INTENT( OUT ) :: d_phi
REAL( KIND=8 ), INTENT( OUT ) :: d_p_tilde( 2 ), d_p_hat( 2 ), d_q_tilde( 2 ), d_q_hat( 2 )
REAL( KIND=8 ) :: c_mat( 2, 2 ), csq
d_phi = ( 1.0 / c ) * ray%t( 3 ) * &
( ray%t( 2 ) * gradc( 1 ) - ray%t( 1 ) * gradc( 2 ) ) / &
( ray%t( 1 ) ** 2 + ray%t( 2 ) ** 2 )
csq = c ** 2
c_mat( 1, : ) = -[ cnn, cmn ] / csq
c_mat( 2, : ) = -[ cmn, cmm ] / csq
d_p_tilde = MATMUL( c_mat, ray%q_tilde )
d_p_hat = MATMUL( c_mat, ray%q_hat )
d_q_tilde = c * ray%p_tilde
d_q_hat = c * ray%p_hat
! d_f = ( c0 * ray0%DetP - cnn0 / csq0 * ray0%DetQ )
! d_g = ( c0 * ray0%DetP - cmm0 / csq0 * ray0%DetQ )
! d_h = ( - cmn0 / csq0 * ray0%DetQ )
! d_DetP = 1.0D0 / csq0 * ( -cmm0 * ray0%f - cnn0 * ray0%g + 2.0 * cmn0 * ray0%h )
! d_DetQ = c0 * ( ray0%f + ray0%g )
END SUBROUTINE ComputeDeltaPQ