program injbump

  use bmad
  use cesr_basic_mod

  implicit none

  type (lat_struct) lat
  type (cesr_struct) cesr
  type (ele_struct) ele 

  integer ix_26, ix_28, ix_36, nargs
  integer n_save
  integer ix_cache
  integer iu,iu1
  integer ix
  character*120 current_lat, lattice
  character*80 lat_file
  character*120 line, last_line

  real(rp) k26,k28,k36
  real(rp) x26_36,   xp26_36,   x28_36,   xp28_36

  call mpm_init('INJ')

  call getlat (current_lat) 
  call choose_cesr_lattice (lattice, lat_file, current_lat, lat)

  call twiss_at_start(lat)
  call twiss_propagate_all(lat)

  call element_locator('bump_26w',lat,ix_26)
  call element_locator('bump_28w',lat,ix_28)
  call element_locator('bump_36w',lat,ix_36)

  ! x from 26 at 36 is
  x26_36 = sqrt(lat%ele(ix_26)%a%beta * lat%ele(ix_36)%a%beta) * &
       sin( lat%ele(ix_36)%a%phi - lat%ele(ix_26)%a%phi)
  xp26_36 = sqrt(lat%ele(ix_26)%a%beta / lat%ele(ix_36)%a%beta) * &
       cos( lat%ele(ix_36)%a%phi - lat%ele(ix_26)%a%phi)
  ! x from 28 cancels x from 26
  x28_36 = sqrt(lat%ele(ix_28)%a%beta * lat%ele(ix_36)%a%beta) * &
       sin( lat%ele(ix_36)%a%phi - lat%ele(ix_28)%a%phi)
  k28 = -x26_36/x28_36
  ! xp from 36 cancels xp from 26 and 28
  xp28_36 = k28*sqrt(lat%ele(ix_28)%a%beta / lat%ele(ix_36)%a%beta) * &
       cos( lat%ele(ix_36)%a%phi - lat%ele(ix_28)%a%phi)
  k36 = -(xp26_36 + xp28_36)   
  k26=1.

  print *
  print '(12x,3a12)','  Horz Beta ','   Phase   ',' Bumper coef'
  print '(a12,3e12.4)',' Bumper 26W ',lat%ele(ix_26)%a%beta, lat%ele(ix_26)%a%phi, k26
  print '(a12,3e12.4)',' Bumper 28W ',lat%ele(ix_28)%a%beta, lat%ele(ix_28)%a%phi, k28
  print '(a12,3e12.4)',' Bumper 36W ',lat%ele(ix_36)%a%beta, lat%ele(ix_36)%a%phi, k36

  call element_locator('bump_26e',lat,ix_26)
  call element_locator('bump_28e',lat,ix_28)
  call element_locator('bump_36e',lat,ix_36)

  ! x from 26 at 36 is
  x26_36 = sqrt(lat%ele(ix_26)%a%beta * lat%ele(ix_36)%a%beta) * &
       sin( lat%ele(ix_36)%a%phi - lat%ele(ix_26)%a%phi)
  xp26_36 = sqrt(lat%ele(ix_26)%a%beta / lat%ele(ix_36)%a%beta) * &
       cos( lat%ele(ix_36)%a%phi - lat%ele(ix_26)%a%phi)
  ! x from 28 cancels x from 26
  x28_36 = sqrt(lat%ele(ix_28)%a%beta * lat%ele(ix_36)%a%beta) * &
       sin( lat%ele(ix_36)%a%phi - lat%ele(ix_28)%a%phi)
  k28 = -x26_36/x28_36
  ! xp from 36 cancels xp from 26 and 28
  xp28_36 = k28*sqrt(lat%ele(ix_28)%a%beta / lat%ele(ix_36)%a%beta) * &
       cos( lat%ele(ix_36)%a%phi - lat%ele(ix_28)%a%phi)
  k36 = -(xp26_36 + xp28_36)   
  k26=1.


  print '(/,a12,3e12.4)',' Bumper 26E ',lat%ele(ix_26)%a%beta, lat%ele(ix_26)%a%phi, k26
  print '(a12,3e12.4)',' Bumper 28E ',lat%ele(ix_28)%a%beta, lat%ele(ix_28)%a%phi, k28
  print '(a12,3e12.4)',' Bumper 36E ',lat%ele(ix_36)%a%beta, lat%ele(ix_36)%a%phi, k36

  print '(/,a)',' Polarity of 26 is hardwired positive '
  print '(a)',' Polarity of 36 is hardwired negative '
  print *
  print '(a)', ' To change the polarity of 28W: FFF VXPUTN CSR BUMP POL 1 1 n '
  print '(a)', '   where n=0 for negative and n=1 for positive polarity ' 
  print '(a)', ' To change the polarity of 28E: FFF VXPUTN CSR BUMP POL 2 2 n '
  print *

  call mpm_goodbye
  call exit

end program injbump