anchor.diff:    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
  +    ejv(1:napx) = e0
       call part_updatePartEnergy(1,.false.)
     case(dist_fmtE)
  -    ejv(1:napx)  = dist_partCol6(1:napx)
  +    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
  +    ejv(1:napx) = (one + dist_partCol6(1:napx))*e0
       call part_updatePartEnergy(1,.false.)
     case(dist_fmtDPP0)
  -    dpsv(1:napx) = dist_partCol6(1:napx)*c1e3
  +    dpsv(1:napx) = dist_partCol6(1:napx)
       call part_updatePartEnergy(3,.false.)
     case(dist_fmtPT)
  -    ejv(1:napx)  = e0 + dist_partCol6(1:napx)*e0f
  +    ejv(1:napx) = (dist_partCol6(1:napx)*e0f)*beta0 + e0
       call part_updatePartEnergy(1,.false.)
     case(dist_fmtPZN)
       doNormal = .true.
  
positive.code:     call dist_appendFormat(dist_fmtP,        uFac, 6)
    case("DE/E0","DEE0")
      call dist_appendFormat(dist_fmtDEE0,     one,  6)
    case("DP/P0","DPP0")
      call dist_appendFormat(dist_fmtDPP0,     one,  6)
    case("DE/P0","DEP0","PT")
      call dist_appendFormat(dist_fmtPT,       one,  6)
  
negative.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) = -beta0*(dist_partCol5(1:napx)*clight)
    case(dist_fmtZN)
      doNormal = .true.
    case(dist_fmtJZ)
  