!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------

subroutine show_lightsources (ring, lines, nl)

  use cesrv_struct
  use cesr_element_mod
  use twiss_and_track_mod

  implicit none

  type (lat_struct)  ring, workr
  character(*) lines(*)
  integer nl

  type (ele_struct)  ele
  type (normal_modes_struct) modes
  type (coord_struct),allocatable :: orb(:)
  type (coord_struct) :: ls_orb

  integer i, loop

  character(10) elab
 
  type (ring_master_struct)  ls_struct(20)
  integer n_ls
  real(rp) sigma_energy_sq, sigma_h, eta_sq, sigma_v, etap_sq
  real(rp) diverg_v, diverg_h !, phase_angle_h, phase_angle_v
  real(rp) s, cbar(2,2)
!

  call reallocate_coord( orb, ring%n_ele_max )

  write (lines(nl+1), *) " Model coupling is used unless it looks non-physically small,"
  write (lines(nl+2), *) " then coupling is assumed to be 1%"
  write (lines(nl+3), *) " "

  workr = ring

  workr%param%particle = positron$
  write (lines(nl+4), *) " Positrons:"
  nl = nl + 4

  do loop = 1,2

    if (loop == 2) then
      workr%param%particle = electron$
      write (lines(nl+1), *) " Electrons:"
      nl = nl + 1
    endif
    call twiss_and_track(workr, orb)

    call radiation_integrals (workr, orb, modes)

    if (loop == 2) then
      call cesr_elements_get('ELS*',n_ls,ls_struct)
    else
      call cesr_elements_get('PLS*',n_ls,ls_struct)
    endif

    write (lines(nl+1), *) " Emittance: ", modes%a%emittance
    nl = nl + 1

    do i=1,n_ls
      s = ls_struct(i)%z
      call twiss_and_track_at_s(workr, s, ele, orb, ls_orb)
      elab = ls_struct(i)%name 
      sigma_energy_sq = (modes%sigE_E)**2 !espread sq for H size
      eta_sq = ele%a%eta**2         !eta sq
      etap_sq= ele%a%etap**2        !etap sq
      sigma_h = sqrt(modes%a%emittance*ele%a%beta + sigma_energy_sq*eta_sq)   

      diverg_h = sqrt(modes%a%emittance*ele%a%gamma + sigma_energy_sq*etap_sq)

      call c_to_cbar( ele, cbar )
! assumes 1% coupling if cbar unreasonably small
      if (abs(cbar(1,2)) < 1e-3) then
         write (lines(nl+1), '(1x,a, a,e11.4,a)') elab, " cbar12 is unreasonable: ", &
              cbar(1,2), "  Using 1% coupling."
         nl = nl + 1
         sigma_v = sqrt(0.01*modes%a%emittance*ele%b%beta)
         diverg_v = sqrt(0.01*modes%a%emittance*ele%b%gamma)
       else
         sigma_v = sqrt(abs(cbar(1,2))*modes%a%emittance*ele%b%beta)
         diverg_v = sqrt(abs(cbar(1,2))*modes%a%emittance*ele%b%gamma)
      end if

! I don't think most people are interested in phase_angle, so take out
!      phase_angle_h = atan( 2 * ele%a%alpha / (ele%a%gamma - ele%a%beta) ) / 2
!      phase_angle_v = atan( 2 * ele%b%alpha / (ele%b%gamma - ele%b%beta) ) / 2


      write (lines(nl+1), *) " Label      sigma_x   ang_div_x    sigma_y    ang_div_y   "
      write (lines(nl+2), '(1x, a10, es10.3, 2x, es10.3, 2x, es10.3, 2x, es10.3)') &
          elab, sigma_h, diverg_h, sigma_v, diverg_v
      write (lines(nl+3), *) "            beta_a     beta_b      gamma_x     gamma_y      cbar12"
      write (lines(nl+4), '(11x, es10.3, 2x, es10.3, 2x, es10.3, 2x, es10.3, 2x, es10.3)') &
           ele%a%beta, ele%b%beta, ele%a%gamma, ele%b%gamma, cbar(1,2)
      write (lines(nl+5), *) "               x           x'          y           y'          s"
      write (lines(nl+6), '(11x,es10.3,2x,es10.3,2x,es10.3,2x,es10.3,2x,f10.4)') &
           ls_orb%vec(1), ls_orb%vec(2), &
           ls_orb%vec(3), ls_orb%vec(4), s
      write (lines(nl+7), *) " "
      nl = nl + 7
    enddo

  enddo

  deallocate( orb )

end subroutine show_lightsources