subroutine osc_parameters(lat, ix_start,ix_end, co,   osc_param, verbose)
  use bmad
  use da_mod
  use da_interface, exclude => osc_parameters
  ! use mode3_mod

 implicit none
 type (lat_struct) lat
 type (coord_struct), allocatable :: co(:)
! type (normal_modes_struct) mode
 type (taylor_struct) t_map(6)
 type (osc_param_struct) osc_param
 character*16, save :: osc_start, osc_end

 integer ix_start, ix_end
 integer i

 real(rp), allocatable:: ft_matrix(:,:)
 real(rp) mat(6,6), detmn_dag
 real(rp) beta_p, alpha_p, gamma_p , coordinates(4) !, sig2dsp,sig2dse,sig2ds
 real(rp) beta_py, alpha_py, gamma_py
 real(rp)wavelength
 real(rp), save ::  k, mu0
 real(rp) emitx, emity, emitz, det
 real(rp) emitx_0/2.e-9/, emity_0/.5e-9/, emitz_0/1.e-6/
 real(rp) sum_theta_x2, sum_theta_y2
 real(rp) T(6,6)
 logical first/.true./, verbose
 logical err_flag/.false./
! locate start and end elements and compute matrix from start to end
  
if(first)then
   osc_start = lat%ele(ix_start)%name
   osc_end = lat%ele(ix_end)%name
  print *,' osc_start = ',osc_start,' osc_end = ', osc_end
   k=twopi/osc_param%wavelength
   mu0=2.405
  print '(2(a16,i10),a16,es12.4)',' ix_start = ',ix_start,' ix_end = ',ix_end, ' wavelength = ', osc_param%wavelength

   coordinates(1:3) = lat%ele(ix_end)%floor%r(1:3)-lat%ele(ix_start)%floor%r(1:3)
   coordinates(4) =   lat%ele(ix_end)%s-lat%ele(ix_start)%s
  first=.false.
endif

   if(osc_param%linear)call transfer_line_matrix_calc(lat,osc_start,osc_end, mat, detmn_dag)

   beta_p = lat%ele(ix_start)%a%beta
   alpha_p = lat%ele(ix_start)%a%alpha
   gamma_p = lat%ele(ix_start)%a%gamma
   beta_py = lat%ele(ix_start)%b%beta
   alpha_py = lat%ele(ix_start)%b%alpha
   gamma_py = lat%ele(ix_start)%b%gamma
   osc_param%tmat(:,:) = mat(:,:)
   osc_param%tilde_M56 = mat(5,1)*lat%ele(ix_start)%x%eta +mat(5,2)*lat%ele(ix_start)%x%etap +mat(5,6)

   if(osc_param%t2mat)then
      call transfer_map_calc(lat, t_map, err_flag, ix1= ix_start, ix2=ix_end, concat_if_possible = .true.) 
      osc_param%t5mat(1,1) = taylor_coef(t_map(5),taylor_expn([1,1]))
      osc_param%t5mat(1,2) = taylor_coef(t_map(5),taylor_expn([1,2]))
      osc_param%t5mat(1,6) = taylor_coef(t_map(5),taylor_expn([1,6]))
      osc_param%t5mat(2,2) = taylor_coef(t_map(5),taylor_expn([2,2]))
      osc_param%t5mat(2,6) = taylor_coef(t_map(5),taylor_expn([2,6]))
      osc_param%t5mat(3,3) = taylor_coef(t_map(5),taylor_expn([3,3]))
      osc_param%t5mat(3,4) = taylor_coef(t_map(5),taylor_expn([3,4]))
      osc_param%t5mat(3,6) = taylor_coef(t_map(5),taylor_expn([3,6]))
      osc_param%t5mat(4,4) = taylor_coef(t_map(5),taylor_expn([4,4]))
      osc_param%t5mat(4,6) = taylor_coef(t_map(5),taylor_expn([4,6]))
      osc_param%t5mat(6,6) = taylor_coef(t_map(5),taylor_expn([6,6]))
      t(:,:) = osc_param%t5mat(:,:)
      osc_param%second_order_x = 2*(T(1,1)*beta_p - T(1,2)*alpha_p + T(2,2)*gamma_p)**2 +(T(1,2)**2-4*T(1,1)*T(2,2))
      osc_param%second_order_y = 2*(T(3,3)*beta_py - T(3,4)*alpha_py + T(4,4)*gamma_py)**2 +(T(3,4)**2-4*T(3,3)*T(4,4))
   endif
   

   if(osc_param%det)then
    emity=emity_0
    emitz=emitz_0
    do i=0,10
       emitx = emitx_0*i**2
       if(i == 0)emitx = emitx_0/100.
       call nonlin_map_calc(lat,ix_start, ix_end, emitx,emity,emitz, co, det, ft_matrix)
            print '(4(a9,es12.4))',' emitx = ', emitx,' M(5,1) =', ft_matrix(5,1),' M(5,2) =', ft_matrix(5,2),' M(5,6) =', ft_matrix(5,6)
     end do

   endif
   if(osc_param%theta2)then
      call int_theta2(lat, ix_start, ix_end, sum_theta_x2, sum_theta_y2)
      osc_param%int_theta_x2 = sum_theta_x2
      osc_param%int_theta_y2 = sum_theta_y2
   endif

   osc_param%osc_dp_max = mu0/abs(k*osc_param%tilde_M56)
!   sig2dsp = (mode%sige_e)**2 * tilde_M56**2
   osc_param%osc_emit_max = mu0**2/(k**2 * (beta_p*mat(5,1)**2-2*alpha_p*mat(5,1)*mat(5,2)+gamma_p*mat(5,2)**2))
!   sig2ds = sig2dsp + sig2dse
   osc_param%delta_s =     lat%ele(ix_end)%s-lat%ele(ix_start)%s -abs(lat%ele(ix_end)%floor%r(3)-lat%ele(ix_start)%floor%r(3))
  if(verbose)then
   print *
   print *, 'Optical stochastic cooling insert parameters'

   print *,' osc_start = ',osc_start,' osc_end = ', osc_end
   print '(a,3es12.4)',' beta/alpha/gamma = ', beta_p,alpha_p,gamma_p
!   print '(a,es12.4)',' Horizontal emittance = ', mode%a%emittance
!   print '(a,es12.4)',' Fractional energy spread = ', mode%sige_e
   print '(4(a10,1x,es12.4,1x))', 'M51 = ',mat(5,1),'M52 = ',mat(5,2),'M56 = ',mat(5,6),'tilde_M56 = ', osc_param%tilde_M56
!   print '(2(a10,1x,es12.4,1x))', 'sig2_DS = ',sig2ds
   print '(a,4es12.4)',' coordinates  Delta x, Delta y, Delta z, Delta s = ',coordinates(1:4)
   print '(a,es12.4)',' Delta s = ',osc_param%delta_s
   print '(2(a16,es12.4,10x))',' emit_max = ', osc_param%osc_emit_max,' dp/p_max = ',osc_param%osc_dp_max

   print '(3(a10,es12.4))','T51 = ',taylor_coef(t_map(5),[1,0,0,0,0,0]),'T52 = ',taylor_coef(t_map(5),[0,1,0,0,0,0]),'T56 = ',taylor_coef(t_map(5),[0,0,0,0,0,1])
   print '(3(a10,es12.4))','T511= ',taylor_coef(t_map(5),taylor_expn([1,1])),'T512= ',taylor_coef(t_map(5),taylor_expn([1,2])),'T516= ',taylor_coef(t_map(5),taylor_expn([1,6]))
   print '(3(a10,es12.4))','T522= ',taylor_coef(t_map(5),taylor_expn([2,2])),'T526= ',taylor_coef(t_map(5),taylor_expn([2,6])),'T566= ',taylor_coef(t_map(5),taylor_expn([6,6]))
   print '(2(a10,es12.4))','T53 = ',taylor_coef(t_map(5),[0,0,1,0,0,0]),'T54 = ',taylor_coef(t_map(5),[0,0,0,1,0,0])
   print '(3(a10,es12.4))','T533= ',taylor_coef(t_map(5),taylor_expn([3,3])),'T534= ',taylor_coef(t_map(5),taylor_expn([3,4])),'T536= ',taylor_coef(t_map(5),taylor_expn([3,6]))
   print '(2(a10,es12.4))','T544= ',taylor_coef(t_map(5),taylor_expn([4,4])),'T546= ',taylor_coef(t_map(5),taylor_expn([4,6]))
   print *
   print '(a,es12.4,a, es12.4)',' sum theta_x2 = ', sum_theta_x2,' sum theta_y2 = ', sum_theta_y2
   print '(2(a,es12.4))',' second_order_x = ' , osc_param%second_order_x,' second_order_y = ' , osc_param%second_order_y
endif
 return 
 end