program twiss_track_test

  use bmad

  implicit none

  type (ring_struct)  ring, ring2
  type (ele_struct) ele
  type (coord_struct), allocatable :: orb_(:)   
  type (modes_struct) mode

  real(rp) chrom_x, chrom_y, delta_e
  real(rp) m6(6,6), vec(6)

  integer i, j, n, n_lines, version

  character(40) lattice
  character(200) lat_file
  character(80), pointer :: lines(:)

!  lat_file = "/home/cesrulib/cesr_libs/current/config/bmad/lat/" // &
!                                                "bmad_L9A18A000-_MOVEREC.lat"

  lat_file = "bmad_L9A18A000-_MOVEREC.lat"
  call bmad_parser (lat_file, ring2)
  call write_digested_bmad_file ('digested.file', ring2, 0)
  call read_digested_bmad_file ('digested.file', ring2, version)
  ring = ring2

  print *, 'Number of elements:', ring%n_ele_max
  print *, 'Precision:', rp

  allocate (orb_(0:ring%n_ele_max))     

  orb_(0)%vec = 0
  call twiss_at_start (ring)
  call closed_orbit_calc (ring, orb_, 4)
  call track_all (ring, orb_)
  call ring_make_mat6 (ring, -1, orb_)
  call twiss_at_start (ring)
  call twiss_propagate_all (ring)

  call type_ele (ring%ele_(96), .false., 6, .true., radians$, .true., ring)
  print *
  call type_coord (orb_(96))

!

  allocate (ring%ele_(96)%descrip)
  ring%ele_(96)%descrip = 'First'

  ele = ring%ele_(96)
  ele%descrip = 'Second'

  print *, 'Descriptions Should be: First, Second.'
  print *, '    ', trim(ring%ele_(96)%descrip)
  print *, '    ', trim(ele%descrip)

  delta_e = 0.0
  call chrom_calc (ring, delta_e, chrom_x, chrom_y)
  print *
  print *, 'Chrom:', chrom_x, chrom_y

!

  call mat_make_unit (m6)
  open (1, file = 'twiss.out')
  do n = 1, ring%n_ele_ring
    write (1, *) '!------------------------------------'
    write (1, *) 'Index:', n
    call type2_ele (ring%ele_(n), lines, n_lines, .false., 0, .false., 0, .false., ring)  
    do i = 1, n_lines
      write (1, '(a)') lines(i)
    enddo
    m6 = matmul (ring%ele_(n)%mat6, m6)
    do i = 1, 6
      write (1, '(6f11.5)') (m6(i, j), j = 1, 6)
    enddo
  enddo

!

  call radiation_integrals (ring, orb_, mode)

  print *
  print '(a,1p,6e18.10)', 'Synch_int1/3:', mode%synch_int
  print '(a,1p,6e18.10)', 'Sig_e/z:     ', mode%sige_e, mode%sig_z, mode%e_loss
  print '(a,1p,6e18.10)', 'a', mode%a%emittance, mode%a%synch_int(4:5)
  print '(a,1p,6e18.10)', 'b', mode%b%emittance, mode%b%synch_int(4:5)
  print '(a,1p,6e18.10)', 'z', mode%z%emittance, mode%z%synch_int(4)

!
 
  print *
  print *, 'TEST OK IF NO MESSAGES BELOW:'

  ele = ring%ele_(96)

  if (ele%name /= 'B11W') print *, 'ERROR: ELEMENT NAME IS NOT CORRECT!'
  if (ele%key /= sbend$)  print *, 'ERROR: ELEMENT KEY IS NOT CORRECT!'

  call check (ele%x%beta,   5.74803d0, 1D-4, 'BETA_X')
  call check (ele%x%alpha,  0.61270d0, 1D-4, 'ALPHA_X')
  call check (ele%x%eta,    0.37872d0, 1D-4, 'ETA_X')
  call check (ele%x%phi,    6.67082d0, 1D-4, 'PHI_X')
  call check (ele%y%beta,  30.04289d0, 1D-4, 'BETA_y')
  call check (ele%y%alpha, -1.70015d0, 1D-4, 'ALPHA_y')
  call check (ele%y%eta,   -0.02685d0, 1D-4, 'ETA_y')
  call check (ele%y%phi,    5.93839d0, 1D-4, 'PHI_y')

  vec = orb_(96)%vec
  call check (vec(1),  3.29560D-3, 1d-7, 'ORB X')
  call check (vec(2), -1.48903D-3, 1d-7, 'ORB P_X')
  call check (vec(3), -0.05019D-3, 1d-7, 'ORB Y')
  call check (vec(4), -0.00322D-3, 1d-7, 'ORB P_Y')
  call check (vec(5), -1.08662D-3, 1d-7, 'ORB Z')
  call check (vec(6), -0.00276D-3, 1d-7, 'ORB P_Z')

  call check (chrom_x, -1.14165_rp, 1d-4, 'CHROM_X')
  call check (chrom_y, -0.96340_rp, 1d-4, 'CHROM_Y')
 
  call check (mode%synch_int(1), 8.8420544D+00, 1d-6, 'SYNCH_INT(1)')
  call check (mode%synch_int(2), 0.1037476D-00, 1d-6, 'SYNCH_INT(2)')
  call check (mode%synch_int(3), 0.0023000D-00, 1d-6, 'SYNCH_INT(3)')

  call check (mode%sige_e,  6.78028304568D-04, 1D-10, 'SIG_E')
  call check (mode%sig_z,  1.8495854827D-02, 1D-8, 'SIG_Z')
  call check (mode%e_loss, 1.1434987872D+06, 1D-1,  'E_LOSS')
  call check (mode%a%emittance,  2.0453812467D-07, 1D-12, 'A%EMITTANCE')
  call check (mode%b%emittance,  1.8021489297D-11, 1D-14, 'B%EMITTANCE')
  call check (mode%z%emittance,  1.2540715133D-05, 1D-11, 'Z%EMITTANCE')

  call check (mode%a%synch_int(4), -1.0619792D-3, 1d-7, 'A%SYNCH_INT(4)')
  call check (mode%a%synch_int(5),  5.2111937D-4, 1d-7, 'A%SYNCH_INT(5)')

  call check (mode%b%synch_int(4), -6.150725623D-4, 1d-7, 'B%SYNCH_INT(4)')
  call check (mode%b%synch_int(5),  4.571914425D-8, 1d-11, 'B%SYNCH_INT(5)')

  call check (mode%z%synch_int(4), -1.676986D-3, 1D-8, 'Z%SYNCH_INT(4)')

!--------------------------------------------------------------------
contains

subroutine check (now, theory, err, what)

  implicit none

  real(rp) now, theory, err
  character(*) what

!

  if (abs(now - theory) > err) then
    print *
    print *, 'ERROR: BAD ', what, '!'
    print *, '       NOW, THEORY:   ', now, theory
    print *, '       RATIO:', now/theory
  endif

end subroutine

end program