anchor.diff:    case(dist_fmtPX)
       yv1(1:napx) = (dist_partCol2(1:napx)/ejfv(1:napx))*c1e3
     case(dist_fmtPXP0)
  -    yv1(1:napx) = (dist_partCol2(1:napx)/e0f)*c1e3
  +    yv1(1:napx) = ((dist_partCol2(1:napx)*e0f)/ejfv(1:napx))*c1e3
     case(dist_fmtPXN)
       doNormal = .true.
     case(dist_fmtPhiX)
  
positive.code:   case(dist_fmtNONE)
      sigmv(1:napx) = zero
    case(dist_fmtZETA)
      sigmv(1:napx) = dist_partCol5(1:napx)*rvv(1:napx)
    case(dist_fmtSIGMA)
      sigmv(1:napx) = dist_partCol5(1:napx)
    case(dist_fmtDT)
      sigmv(1:napx) = -(e0f/e0)*(dist_partCol5(1:napx)*clight)
    case(dist_fmtZN)
      doNormal = .true.
    case(dist_fmtJZ)
  
negative.code:   use numerical_constants
  
    logical doAction, doNormal
    real(kind=fPrec) beta0
  
    doAction = .false.
    doNormal = .false.
  
    beta0 = e0f/e0
  
    ! Forward energy/momentum must be calculated first
    select case(dist_partFmt(6))
    case(dist_fmtNONE)
      ejv(1:napx) = e0
      call part_updatePartEnergy(1,.false.)
    case(dist_fmtE)
      ejv(1:napx) = dist_partCol6(1:napx)
      call part_updatePartEnergy(1,.false.)
    case(dist_fmtP)
      ejfv(1:napx) = dist_partCol6(1:napx)
      call part_updatePartEnergy(2,.false.)
    case(dist_fmtDEE0)
      ejv(1:napx) = (one + dist_partCol6(1:napx))*e0
      call part_updatePartEnergy(1,.false.)
    case(dist_fmtDPP0)
      dpsv(1:napx) = dist_partCol6(1:napx)
      call part_updatePartEnergy(3,.false.)
    case(dist_fmtPT)
      ejv(1:napx) = (dist_partCol6(1:napx)*e0f)*beta0 + e0
      call part_updatePartEnergy(1,.false.)
    case(dist_fmtPZN)
      doNormal = .true.
  