subroutine cbar_v_e (lat_in, ele_names, slopes)

  use bmad

  implicit none

  type (lat_struct)  lat_in
  type (lat_struct),save :: lat
  type (ele_struct), pointer :: ele1, ele2

  real(rp) slopes(4)
  character*40 ele_names(4)

  real(rp) l_min, l_max, circum

  integer ele_ix(4), i, ix_min, ix_max
  integer ix, k, j
  logical, save :: first_time = .true.
  character*12 start_name, end_name
  real(rp), allocatable :: gamma_fact(:)
  real(rp), allocatable :: C(:,:,:), Cbar(:,:,:)
  real(rp) energy, del_phi_a, del_phi_b, beta_a, beta_b, alpha_a, alpha_b
  real(rp) iden(2,2)
  real(rp), allocatable :: e(:)
  real(rp) mat(6,6)
  real(rp) detmn_dag
  real(rp) t0(4,4),U(4,4), V(4,4)
  real(rp) Ubar(4,4), Vbar(4,4), Ga(4,4), Gb_inv(4,4)
  type (ele_struct) ele
  real(rp) C_star(2,2)
  real(rp) H(2,2)
  real(rp) detH

  integer num_energy, n_param, n
  real(rp) :: start_energy = -0.0002, end_energy = 0.0002, inc_energy = 0.0001

  type(coord_struct), allocatable, save :: coord(:)

  real(rp) a,b,siga,sigb,chi2,q
  
  lat = lat_in
  if(first_time)then
     n_param = 0
     do i=1,size(ele_names)
        call string_trim(ele_names(i), ele_names(i),ix)
        if(ix >0)then
           n_param=n_param+1
        endif
     end do
     
     do n=1,n_param
        ele_ix(n)=0
        do i=1,lat%n_ele_max
           if(lat%ele(i)%name == ele_names(n))then
              ele_ix(n)=i
              exit
           endif
        end do
        if(ele_ix(n) == 0)then
           print *,' can t find element ', ele_names(n)
           stop
        endif
     end do     

     circum = lat%param%total_length
     l_min = circum
     l_max = 0
     do i = 1, n_param
        ix = ele_ix(i)
        if (lat%ele(ix)%n_slave /= 0) then
           ele1 => pointer_to_slave(lat%ele(ix), 1)
           ele2 => pointer_to_slave(lat%ele(ix), lat%ele(ix)%n_slave)
           
           do j = ele1%ix_ele, ele2%ix_ele
              if(lat%ele(j)%s < l_max)cycle
              if(lat%ele(j)%s > l_min)cycle
              if(lat%ele(j)%s < circum/2) then
                 l_max = lat%ele(j)%s
                 ix_max = j
                 end_name = lat%ele(j)%name
              else
                 l_min = lat%ele(j)%s
                 ix_min = j
                 start_name = lat%ele(j)%name
              endif
           end do
        else
           if(lat%ele(ix)%s < l_max)cycle
           if(lat%ele(ix)%s > l_min)cycle
           if(lat%ele(ix)%s < circum/2) then
              l_max = lat%ele(ix)%s
              ix_max = ix
              end_name = lat%ele(ix)%name
           else
              l_min = lat%ele(ix)%s
              ix_min = ix
              start_name = lat%ele(ix)%name
           endif
        endif
     end do

     iden(1,1)=1
     iden(1,2)=0
     iden(2,1)=0
     iden(2,2)=1

  endif
  call reallocate_coord(coord,lat%n_ele_max)
  
  call twiss_at_start(lat)
  coord(0)%vec(:)=0.
  call closed_orbit_calc( lat, coord, 4 )
  call track_all (lat, coord)
  call lat_make_mat6(lat,-1,coord)
  call twiss_at_start(lat)
  call twiss_propagate_all(lat)  
  
  
  del_phi_a = 2*lat%ele(ix_max)%a%phi
  del_phi_b = 2*lat%ele(ix_max)%b%phi
  
  num_energy=int((end_energy - start_energy)/inc_energy) + 1     
  allocate(gamma_fact(num_energy))
  allocate(C(num_energy,2,2))
  allocate(Cbar(num_energy,2,2))      
  allocate(e(num_energy))
  
  i=1
  energy = start_energy
  do while(energy .le. end_energy)
     e(i)=energy
     energy = energy + inc_energy
     i = i + 1
  enddo
    
  k=1
    
  do while(k .le. num_energy)
     
     do i=0,lat%n_ele_max
        coord(i)%vec(6) = e(k)
     end do
     
     call lat_make_mat6(lat,-1,coord)
     
     call transfer_line_matrix_calc(lat,start_name,end_name, mat, detmn_dag)
     
     H(1,1)=mat(1,3)+mat(4,2)
     H(1,2)=mat(1,4)-mat(3,2)
     H(2,1)=mat(2,3)-mat(4,1)
     H(2,2)=mat(2,4)+mat(3,1)
     detH= H(1,1)*H(2,2)-H(1,2)*H(2,1)
     
     if(((mat(1,1)-mat(3,3)+mat(2,2)-mat(4,4))**2+4*detH).lt.0) then 
        print *, "CBAR_V_E: No real solutions exist"
        stop
     end if
     
     gamma_fact(k) = sqrt(0.5 + 0.5 * sqrt((mat(1,1)-mat(3,3)+mat(2,2)-mat(4,4))**2 / &
                            ((mat(1,1)-mat(3,3)+mat(2,2)-mat(4,4))**2+4*detH)))
     
     C(k,:,:)= -H/(gamma_fact(k)*sqrt((mat(1,1)-mat(3,3)+mat(2,2)-mat(4,4))**2+4*detH))
     
     if((mat(1,1)-mat(3,3)+mat(2,2)-mat(4,4)).lt.0) C(k,:,:)=-C(k,:,:)
     C_star(1,1)=C(k,2,2)
     C_star(2,1)=-C(k,2,1)
     C_star(1,2)=-C(k,1,2)
     C_star(2,2)=C(k,1,1)
     
     V(1:2,1:2)=gamma_fact(k)*iden
     V(3:4,3:4)=gamma_fact(k)*iden
     V(3:4,1:2)=-C_star
     V(1:2,3:4)=C(k,:,:)
     
     Vbar(1:2,1:2)=gamma_fact(k)*iden
     Vbar(3:4,3:4)=gamma_fact(k)*iden
     Vbar(3:4,1:2)=C_star
     Vbar(1:2,3:4)=-C(k,:,:)
     
     do i=1,4
        do j=1,4
           t0(i,j)=(Vbar(i,1)*mat(1,j)+Vbar(i,2)*mat(2,j)+Vbar(i,3)*mat(3,j)+Vbar(i,4)*mat(4,j))
        end do
     end do
     
     do i=1,4
        do j=1,4
           U(i,j)=(t0(i,1)*V(1,j)+t0(i,2)*V(2,j)+t0(i,3)*V(3,j)+t0(i,4)*V(4,j))
        end do
     end do
     
     beta_a = U(1,2)/sin(del_phi_a)
     beta_b = U(3,4)/sin(del_phi_b)
     alpha_a = (U(1,1)-cos(del_phi_a))/sin(del_phi_a)
     alpha_b = (U(3,3)-cos(del_phi_b))/sin(del_phi_b)
     
     Ga(1,1)=sqrt(beta_a)
     Ga(1,2)=0
     Ga(2,1)=-alpha_a/sqrt(beta_a)
     Ga(2,2)=-1/sqrt(beta_a)
     
     Gb_inv(2,2)=-sqrt(beta_b)
     Gb_inv(1,2)=0
     Gb_inv(2,1)=-alpha_b/sqrt(beta_b)
     Gb_inv(1,1)=1/sqrt(beta_b)
     
     do i=1,2
        do j=1,2
           t0(i,j)=(C(k,i,1)*Gb_inv(1,j)+C(k,i,2)*Gb_inv(2,j))
        end do
     end do
     
     do i=1,2
        do j=1,2
           Cbar(k,i,j)=(Ga(i,1)*t0(1,j)+Ga(i,2)*t0(2,j))
        end do
     end do

     k=k+1
  end do

  call fit(e,Cbar(:,1,1),a,b,siga,sigb,chi2,q)
  slopes(1)=b
  call fit(e,Cbar(:,1,2),a,b,siga,sigb,chi2,q)
  slopes(2)=b
  call fit(e,Cbar(:,2,1),a,b,siga,sigb,chi2,q)
  slopes(3)=b
  call fit(e,Cbar(:,2,2),a,b,siga,sigb,chi2,q)
  slopes(4)=b

  deallocate(e)
  deallocate(gamma_fact)
  deallocate(C)
  deallocate(Cbar)
end subroutine cbar_v_e