anchor.diff:    call alloc(dist_partCol5, napx, zero, "dist_partCol5")
     call alloc(dist_partCol6, napx, zero, "dist_partCol6")
   
  -  dist_tMat(1:6,1:6) = tas(1:6,1:6)
  -  dist_tMat(1:5,6)   = dist_tMat(1:5,6)*c1m3
  -  dist_tMat(6,1:5)   = dist_tMat(6,1:5)*c1e3
  +  ! If the T matrix hasn't been set from input, we set it here
  +  if(dist_normMethod == dist_normSIXTRACK) then
  +    dist_tMat(1:6,1:6) = tas(1:6,1:6)
  +    dist_tMat(1:5,6)   = dist_tMat(1:5,6)*c1m3
  +    dist_tMat(6,1:5)   = dist_tMat(6,1:5)*c1e3
  +  end if
  +  if(dist_normMethod == dist_normINPUT .and. dist_tMatRow /= 6) then
  +    write(lerr,"(a,i0,a)") "DIST> ERROR The T-matrix must be 6 rows, but only ",dist_tMatRow," rows were given"
  +    call prror
  +  end if
  +
  +  ! Scale longitudinal emittance
  +  select case(dist_emit3Unit)
  +  case(1) ! Micrometre
  +    dist_beamEmit(3) = dist_beamEmit(3)*c1m6
  +  case(2) ! eV s
  +    dist_beamEmit(3) = dist_beamEmit(3)*((clight*beta0)/((four*pi)*e0))
  +  end select
  +  write(lout,"(a,1pe13.6,a)") "DIST> EMITTANCE X = ",dist_beamEmit(1)," m"
  +  write(lout,"(a,1pe13.6,a)") "DIST> EMITTANCE Y = ",dist_beamEmit(2)," m"
  +  write(lout,"(a,1pe13.6,a)") "DIST> EMITTANCE Z = ",dist_beamEmit(3)," m"
   
   #ifdef DISTLIB
     if(dist_distLib) then
       call distlib_init(1)
       call distlib_setEnergyMass(e0, nucm0)
       call distlib_setEmittance12(dist_beamEmit(1), dist_beamEmit(2))
  -  ! call distlib_setEmittance3(dist_beamEmit(3))
  -    flatTas = pack(transpose(dist_tMat),.true.)
  -    call distlib_setTasMatrix(flatTas)
  +    call distlib_setEmittance3(dist_beamEmit(3))
  +    select case(dist_normMethod)
  +    case(dist_normSIXTRACK, dist_normINPUT)
  +      do i=1,6
  +        write(lout,"(a,i0,a,6(1x,1pe13.6))") "DIST> TMATRIX(",i,",1:6) =",dist_tMat(i,1:6)
  +      end do
  +      flatTas = pack(transpose(dist_tMat),.true.)
  +      call distlib_setTasMatrix(flatTas)
  +    case(dist_normTWISS)
  +      call distlib_setTwiss(dist_twBeta(1), dist_twAlpha(1), dist_twBeta(2), dist_twAlpha(2))
  +      call distlib_setDispersion(dist_bDisp(1), dist_bDisp(2), dist_bDisp(3), dist_bDisp(4))
  +      write(lout,"(a,2(1x,1pe13.6))") "DIST> BETA  X/Y =",dist_twBeta
  +      write(lout,"(a,2(1x,1pe13.6))") "DIST> ALPHA X/Y =",dist_twAlpha
  +    end select
     end if
   #endif
   
  
positive.code:   use crcoall
    use floatPrecision
    use parpro, only : mFileName
    use numerical_constants, only : zero
  
    implicit none
  
  
negative.code:   integer, parameter :: dist_fillUNIFORM    = 4  ! Uniform distribution
    integer, parameter :: dist_fillLINEAR     = 5  ! Linear fill
    integer, parameter :: dist_fillCOUNT      = 6  ! Integer range
  
    !  Normalisation Methods
    ! =======================
    integer, parameter :: dist_normNONE       = 0  ! No normalisation
    integer, parameter :: dist_normSIXTRACK   = 1  ! Use SixTrack tas matrix
    integer, parameter :: dist_normINPUT      = 2  ! Use matrix in DIST block
    integer, parameter :: dist_normTWISS      = 3  ! Use twiss in DIST block
  
    integer,          private, save :: dist_normMethod  = dist_normSIXTRACK ! Flag for source of normalisation matrix
    integer,          private, save :: dist_tMatRow     = 0                 ! How many T-matrix rows we've read
    real(kind=fPrec), private, save :: dist_tMat(6,6)   = zero              ! T-matrix to use for normalisation
    real(kind=fPrec), private, save :: dist_twBeta(2)   = one               ! Twiss beta
    real(kind=fPrec), private, save :: dist_twAlpha(2)  = zero              ! Twiss alpha
    real(kind=fPrec), private, save :: dist_bDisp(4)    = zero              ! Dispersion
    real(kind=fPrec), private, save :: dist_beamEmit(3) = zero              ! Beam emittance
    integer,          private, save :: dist_emit3Unit   = 0                 ! Flag for converting longitudinal emittance
  
    !  Column Formats
    ! ================
  