subroutine do_synrad (walls, u, ring, gen_params, window, auto)

  use synrad_mod
  use cesrv_struct
  use cesrv_interface
  use synrad_window_mod

  implicit none

  type (walls_struct) walls
  type (universe_struct), target :: u
  type (lat_struct) ring, p_ring, e_ring
  type (coord_struct), allocatable :: e_orb(:), p_orb(:)
  type (synrad_param_struct) gen_params
  type (ele_power_struct), allocatable :: p_power(:), e_power(:)
  type (crotch_window_struct) window(:)
  type (normal_modes_struct) modes

  logical, optional :: auto
  logical :: run_auto

  integer ix, ios
  character*80 line

  run_auto = .false.
  if (present(auto)) then
    run_auto = auto
  end if

  allocate(e_orb(0:ring%n_ele_max))
  allocate( p_orb(0:ring%n_ele_max))
  allocate( e_power(ring%n_ele_max))
  allocate( p_power(ring%n_ele_max))

  gen_params%i_beam = .2

  if (.not. run_auto) then
    call get_input_string ('Current per beam <CR = 0.2 Amps>:', line)
    call string_trim (line, line, ix)
    if (ix /= 0) then
      read (line, *, iostat = ios) gen_params%i_beam
      if (ios /= 0) then
        print *, 'Cannot read number.'
        return
      endif
    endif
  end if

  call radiation_integrals(ring, u%orb, modes)

  gen_params%epsilon_y = max((modes%a%emittance * .02), modes%b%emittance)

  if (.not. run_auto) then
    print *, "a%emit: ", modes%a%emittance, " b%emit: ", modes%b%emittance
    print *, 'Default emittance is: ',gen_params%epsilon_y
    call get_input_string ('Vertical Emittance <CR = DEFAULT>:', line)
    call string_trim (line, line, ix)
    if (ix /= 0) then
      read (line, *, iostat = ios) gen_params%epsilon_y
      if (ios /= 0) then
        print *, 'Cannot read number.'
        return
      endif
      print *, ' Vertical Emittance set to: ',gen_params%epsilon_y
    endif
  end if

! initialize the window ray counters

  call window_ray_reset (window)

! calculate twiss, closed orbits

  if (ring%param%particle == electron$) then
    e_ring = ring
    e_ring%param%particle = electron$
    call closed_orbit_calc (e_ring, e_orb, 4)
    call lat_make_mat6 (e_ring, -1, e_orb)
    call twiss_at_start (e_ring)
    call twiss_propagate_all (e_ring)
    p_ring = ring
    p_ring%param%particle = positron$
    call twiss_at_start (p_ring)
    call closed_orbit_calc (p_ring, p_orb, 4)
    call lat_make_mat6 (p_ring, -1, p_orb)
    call twiss_at_start (p_ring)
    call twiss_propagate_all (p_ring)
  else
    p_ring = ring
    p_ring%param%particle = positron$
    call twiss_at_start (p_ring)
    call closed_orbit_calc (p_ring, p_orb, 4)
    call lat_make_mat6 (p_ring, -1, p_orb)
    call twiss_at_start (p_ring)
    call twiss_propagate_all (p_ring)
    e_ring = ring
    e_ring%param%particle = electron$
    call closed_orbit_calc (e_ring, e_orb, 4)
    call lat_make_mat6 (e_ring, -1, e_orb)
    call twiss_at_start (e_ring)
    call twiss_propagate_all (e_ring)
  endif

  call init_wall (walls%positive_x_wall)
  call init_wall (walls%negative_x_wall)

  call calculate_window_power (p_ring, p_orb, +1, p_power, walls, gen_params, window)
  call calculate_window_power (e_ring, e_orb, -1, e_power, walls, gen_params, window)

  call sigma_at_windows (window, gen_params)

  deallocate(e_orb)
  deallocate(p_orb)
  deallocate(e_power)
  deallocate(p_power)

  print *,' Synchrotron radiation modeling and propagation is complete.'

  !-----------------------------------------------
  contains

  subroutine init_wall (wall)

  use synrad_struct
  use synrad_interface

  implicit none

  type (wall_struct) wall

  !

  wall%seg(:)%power%power_tot = 0
  wall%seg(:)%power%power_per_len = 0
  wall%seg(:)%power%power_per_area = 0
  wall%seg(:)%power%n_source = 0
  wall%seg(:)%power%main_source%ix_ele = 0
  wall%seg(:)%power%main_source%power_per_len = 0
  wall%seg(:)%power%main_source%s = 0

  end subroutine init_wall

end subroutine