subroutine rotate_cbar(ele, tilt, cbar, data_in)

    use cesrv_struct, only: rp, ele_struct
    implicit none

    type(ele_struct) ele
    real(rp) :: cbar(2,2), cbar_bmad(2,2)
    real(rp) :: beta_rel, Arel, Brel, phi_a_rel, phi_b_rel
    real(rp) :: A(2,2), B(2,2) ! A,B((x,y),(in-phase,out-of-phase)), in lab frame
    real(rp) :: A_obs(2,2), B_obs(2,2) ! Aobs,Bobs((x,y),(in-phase,out-of-phase)), in "observed" (rotated) frame
    real(rp) :: A_twid(2,2), B_twid(2,2) ! same notation as A,B above; projected s.th. A(1,2) = 0 and B(2,1) = 0.
    real(rp) :: phi_a_shift, phi_b_shift ! phase shifts necessary to transform A,B --> Atwid, Btwid

    real(rp) :: tilt ! in radians
    real(rp) :: R_tilt(2,2) ! rotation matrix corresponding to the tilt
    real(rp) :: R_phi_a(2,2), R_phi_b(2,2) ! rotation matrices corresponding to phase shifts phi_a_shift, phi_b_shift

    integer :: ix_lat

    logical, optional, intent(in) :: data_in
    logical data

    if (present(data_in)) then
       data = data_in
    else
       data = .false.
    endif

    ix_lat = ele%ix_ele

     beta_rel = ele%a%beta / ele%b%beta

     if (data) then
        cbar_bmad = ele%c_mat
     else
        call c_to_cbar(ele, cbar_bmad)
     endif

     ! define amplitudes of a- and b-mode motion:
     A(1,1) = ele%gamma_c * sqrt(ele%a%beta)
     A(1,2) = 0.
     A(2,1) = -sqrt(ele%b%beta) * cbar_bmad(2,2)
     A(2,2) = -sqrt(ele%b%beta) * cbar_bmad(1,2)

     B(1,1) =  sqrt(ele%a%beta) * cbar_bmad(1,1)
     B(1,2) = -sqrt(ele%a%beta) * cbar_bmad(1,2)
     B(2,1) = ele%gamma_c * sqrt(ele%b%beta)
     B(2,2) = 0.

     ! construct the rotation matrix corresponding to the BPM tilt:
     R_tilt(1,1) =  cos(tilt)
     R_tilt(1,2) =  sin(tilt)
     R_tilt(2,1) = -sin(tilt)
     R_tilt(2,2) =  cos(tilt)

     ! rotate the A and B matrices to generate "observed" (A_obs, B_obs) amplitudes:
     A_obs = matmul(R_tilt, A)
     B_obs = matmul(R_tilt, B)

     ! phase-shift A_obs terms such that A_obs(1,2) = 0
     phi_a_shift = acos(A_obs(1,2) / sqrt(A_obs(1,1)**2 + A_obs(1,2)**2))
     phi_a_shift = sign(phi_a_shift, A_obs(1,1)) - 3.14159265358979323/2. ! if theta = 0, phi_a_shift = -pi/2 and horiz. motion is then cos(psi_a)

     R_phi_a(1,1) =  cos(-phi_a_shift)
     R_phi_a(1,2) =  sin(-phi_a_shift)
     R_phi_a(2,1) = -sin(-phi_a_shift)
     R_phi_a(2,2) =  cos(-phi_a_shift)

     A_twid(1,1) = sqrt(A_obs(1,1)**2 + A_obs(1,2)**2)
     A_twid(1,2) = 0.
     A_twid(2,:) = matmul(R_phi_a, A_obs(2,:)) ! mix the A_obs(2,:) amplitudes to compensate for phase-shift

     ! now do the same for B_obs terms, such that B_obs(2,2) = 0:
     phi_b_shift = acos(B_obs(2,2) / sqrt(B_obs(2,1)**2 + B_obs(2,2)**2))
     phi_b_shift = sign(phi_b_shift, B_obs(2,1)) - 3.14159265358979323/2. ! if theta = 0, phi_b_shift = -pi/2 and vert. motion is then cos(psi_b)

     R_phi_b(1,1) =  cos(-phi_b_shift)
     R_phi_b(1,2) =  sin(-phi_b_shift)
     R_phi_b(2,1) = -sin(-phi_b_shift)
     R_phi_b(2,2) =  cos(-phi_b_shift)

     B_twid(1,:) = matmul(R_phi_b, B_obs(1,:)) ! mix the B_obs(2,:) amplitudes to compensate for phase-shift
     B_twid(2,1) = sqrt(B_obs(2,1)**2 + B_obs(2,2)**2)
     B_twid(2,2) = 0.

     Arel = sqrt(A_twid(2,1)**2 + A_twid(2,2)**2) / sqrt(A_twid(1,1)**2 + A_twid(1,2)**2)
     Brel = sqrt(B_twid(1,1)**2 + B_twid(1,2)**2) / sqrt(B_twid(2,1)**2 + B_twid(2,2)**2)

     phi_a_rel = atan2(A_twid(2,2), A_twid(2,1))
     phi_b_rel = atan2(B_twid(1,2), B_twid(1,1))

     cbar(2,2) = -ele%gamma_c * sqrt(beta_rel)    * Arel * cos(phi_a_rel)
     cbar(1,2) = -ele%gamma_c * sqrt(beta_rel)    * Arel * sin(phi_a_rel)
     cbar(1,1) = ele%gamma_c * sqrt(1./beta_rel) * Brel * cos(phi_b_rel)

  end subroutine rotate_cbar