This is old stuff with conversion between PTC units and Bmad units.
Current code just tells PTC to use Bmad units so PTC handles the bookkeeping.


!------------------------------------------------------------------------
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!+
! Subroutine taylor_bmad_to_ptc (taylor_bmad, beta0, taylor_ptc, mat6)
!
! Routine to convert a Bmad taylor map to PTC taylor.
!
! Input:
!   taylor_bmad(6) -- type(taylor_struct): Bmad Taylor.
!   beta0          -- real(rp): Reference particle velocity
!
! Output:
!   taylor_ptc(6)  -- type(real_8): PTC coordinates.
!-

subroutine taylor_bmad_to_ptc (taylor_bmad, beta0, taylor_ptc)

use s_fibre_bundle

implicit none

type (taylor_struct) taylor_bmad(:)
type (real_8) taylor_ptc(:), t_ptc(6)

real(rp) beta0

! taylor_ptc(5) = (E - E0) / P0c
! taylor_ptc(6) = c (t - t0)
! 1/beta0 + taylor_ptc(5) == E / P0c

t_ptc = taylor_bmad
taylor_ptc = t_ptc
taylor_ptc(5) = (t_ptc(6)**2 + 2*t_ptc(6)) / (1/beta0 + sqrt(1/beta0**2+t_ptc(6)**2+2*t_ptc(6)))
taylor_ptc(6) = -t_ptc(5) * (1/beta0 + taylor_ptc(5)) / (1 + t_ptc(6))

end subroutine taylor_bmad_to_ptc 

!------------------------------------------------------------------------
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!+
! Subroutine taylor_ptc_to_bmad (taylor_ptc, beta0, taylor_bmad)
!
! Routine to convert a PTC real_8 taylor to a Bmad Taylor.
!
! Input:
!   taylor_ptc(6)  -- real_8: PTC taylor.
!   beta0          -- real(rp): Reference particle velocity
!
! Output:
!   taylor_bmad(6) -- taylor_struct: Bmad Taylor.
!-

subroutine taylor_ptc_to_bmad (taylor_ptc, beta0, taylor_bmad)

use s_fibre_bundle

implicit none

type (taylor_struct) taylor_bmad(:)
type (real_8) taylor_ptc(:), t_ptc(6)
real(rp) beta0

!

t_ptc(6) = (2*taylor_ptc(5)/beta0+taylor_ptc(5)**2)/(sqrt(1+2*taylor_ptc(5)/beta0+taylor_ptc(5)**2)+1)
t_ptc(5) = -taylor_ptc(6) * (1 + t_ptc(6)) / (1/beta0 + taylor_ptc(5))

taylor_bmad = t_ptc

end subroutine taylor_ptc_to_bmad 

!------------------------------------------------------------------------
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!+
! Subroutine vec_bmad_to_ptc (vec_bmad, beta0, vec_ptc, conversion_mat)
!
! Routine to convert a Bmad vector map to PTC vector,
!
! Input:
!   vec_bmad(6) -- real(rp): Bmad coordinates.
!   beta0       -- real(rp): Reference particle velocity
!
! Output:
!   vec_ptc(6)     -- real(rp): PTC coordinates.
!   conversion_mat -- real(rp), optional: Jacobian matrix of Bmad -> PTC conversion map.
!-

subroutine vec_bmad_to_ptc (vec_bmad, beta0, vec_ptc, conversion_mat)

implicit none

real(rp) vec_bmad(:), vec_ptc(:)
real(rp) beta0, vec_temp(6)
real(rp), optional :: conversion_mat(6,6)
real(rp) factor1, factor2

! vec_ptc(5) = (E - E0) / P0c
! vec_ptc(6) = c (t - t0)
! 1/beta0 + vec_ptc(5) == E / P0c

vec_temp = vec_bmad
vec_temp(5) = (vec_bmad(6)**2 + 2*vec_bmad(6)) / (1/beta0 + sqrt(1/beta0**2+vec_bmad(6)**2+2*vec_bmad(6)) )
vec_temp(6) = -vec_bmad(5) * (1/beta0 + vec_temp(5)) / (1 + vec_bmad(6))

if (present(conversion_mat)) then
  call mat_make_unit(conversion_mat)
  factor1 = sqrt(1/beta0**2+vec_bmad(6)**2+2*vec_bmad(6))
  factor2 = 1+beta0**2*vec_bmad(6)*(2+vec_bmad(6))
  conversion_mat(5,5) = 0
  conversion_mat(5,6) = beta0**2*(1+vec_bmad(6))*factor1/factor2
  conversion_mat(6,5) = -(1/beta0+beta0*vec_bmad(6)*(2+vec_bmad(6))/(1+beta0*factor1))/(1+vec_bmad(6))
  conversion_mat(6,6) = -((beta0**2-1)*vec_bmad(5)*factor1)/((1+vec_bmad(6))**2*factor2)
end if

vec_ptc = vec_temp

end subroutine vec_bmad_to_ptc 

!------------------------------------------------------------------------
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!+
! Subroutine vec_ptc_to_bmad (vec_ptc, beta0, vec_bmad, conversion_mat, state)
!
! Routine to convert a PTC orbit vector to a Bmad orbit vector.
!
! Input:
!   vec_ptc(6)  -- real(rp): PTC coordinates.
!   beta0       -- real(rp): Reference particle velocity
!
! Output:
!   vec_bmad(6)    -- real(rp): Bmad coordinates.
!   conversion_mat -- real(rp), optional: Jacobian matrix of PTC -> Bmad conversion map.
!   state          -- integer, optional: Set to lost_pz_aperture$ if energy is too low. Set to alive$ otherwise.
!-

subroutine vec_ptc_to_bmad (vec_ptc, beta0, vec_bmad, conversion_mat, state)

implicit none

real(rp) vec_bmad(:), vec_ptc(:)
real(rp) beta0, vec_temp(6)
real(rp), optional :: conversion_mat(6,6)
real(rp) factor1, factor2
integer, optional :: state

!

if (present(state)) state = alive$

factor1 = 1+2*vec_ptc(5)/beta0+vec_ptc(5)**2
if (factor1 <= 0) then
  if (present(state)) state = lost_pz_aperture$
  return
endif

vec_temp = vec_ptc
vec_temp(6) = (2*vec_ptc(5)/beta0+vec_ptc(5)**2)/(sqrt(factor1)+1)
vec_temp(5) = -vec_ptc(6) * (1 + vec_temp(6)) / (1/beta0 + vec_ptc(5))

if (present(conversion_mat)) then
  call mat_make_unit(conversion_mat)
  factor1 = sqrt(1+2*vec_ptc(5)/beta0+vec_ptc(5)**2)
  factor2 = beta0+2*vec_ptc(5)+beta0*vec_ptc(5)**2 
  conversion_mat(5,5) = beta0*(beta0**2-1)*factor1*vec_ptc(6)/((1+beta0*vec_ptc(5))**2*factor2)
  conversion_mat(5,6) = -(1+vec_ptc(5)*(2+beta0*vec_ptc(5))/(beta0*(1+factor1)))/(1/beta0+vec_ptc(5))
  conversion_mat(6,5) = (1+beta0*vec_ptc(5))*factor1/factor2
  conversion_mat(6,6) = 0
end if

vec_bmad = vec_temp

end subroutine vec_ptc_to_bmad 

!------------------------------------------------------------------------
!------------------------------------------------------------------------
!------------------------------------------------------------------------
!+
! Subroutine real_8_to_taylor (y8, beta0, beta1, bmad_taylor)
!
! Routine to convert a PTC real_8 map to a Bmad Taylor map.
! The conversion includes the conversion between Bmad and PTC time coordinate systems.
!
! Input:
!   y8(6)           -- real_8: PTC Taylor map. NOTE: y8 is used as scratch space and therefore trashed.
!   beta0           -- real(rp): Reference particle velocity at beginning of map
!   beta1           -- real(rp): Reference particle velocity at end of map
!   bmad_taylor(6)  -- Taylor_struct: Only %ref is used at input.
!     %ref            -- Reference orbit
!
! Output:
!   bmad_taylor(6) -- Taylor_struct: Bmad Taylor map.
!-

subroutine real_8_to_taylor (y8, beta0, beta1, bmad_taylor)

use s_fibre_bundle

implicit none

type (taylor_struct) :: bmad_taylor(:)
type (real_8) y8(:), rr(6), bet, ss(6)
type (damap) bm, id, si

real(rp) beta0, beta1, fix0(6)

!

call alloc (bm, id, si)
call alloc (rr)
call alloc (bet)
call alloc (ss)

bm = y8

fix0 = bm
id = 1

rr = id + bmad_taylor%ref

ss = rr 
ss(5) = (rr(6)**2+2.d0*rr(6))/(1.d0/beta0 + sqrt(1.d0/beta0**2+rr(6)**2+2.d0*rr(6)))
bet = (1.d0+rr(6))/(1.d0/beta0+ss(5))
ss(6) = -rr(5)/bet

si=ss  ! bmad to ptc map

bm = bm * si
bm = fix0

rr = bm
ss = rr
ss(6) = (2.d0*rr(5)/beta1+rr(5)**2)/(sqrt(1.d0+2.d0*rr(5)/beta1+rr(5)**2)+1.d0)
bet = (1.d0+ss(6))/(1.d0/beta1+rr(5))
ss(5) = -bet*rr(6)

bmad_taylor = ss

call kill (rr)
call kill (ss)
call kill (bet)
call kill (bm, id, si)

end subroutine real_8_to_taylor