subroutine electron_positron (lat_1, lat_pos, lat_ele, co_positron, co_electron, con, global )

! program to compute beta, alpha for electron and positrons

  use nonlin_mod
  use bmadz_mod
  use constraints_mod

  implicit none

  type (lat_struct) lat_pos, lat_ele, lat_1
  type (lat_struct), save :: lat
  type (coord_struct), allocatable, save :: co(:)
  type (coord_struct), allocatable :: co_positron(:), co_electron(:)
  
  type (global_struct) global
  type (constraint_struct) con


  integer i, status
  integer ipretz
  integer, save :: ix_positron = 0, ix_electron = 0
  
  logical :: debug = .false.
  
  call reallocate_coord( co, lat_1%n_ele_max )


  lat = lat_1

  do ipretz = 1,2
    if (ipretz == 1) then
      lat%param%particle = positron$
      co(0)%vec(:) = co_positron(0)%vec(:)
      do i = 1,lat%n_ele_track
        if (lat%ele(i)%key /= beambeam$) cycle
        if (lat%ele(i)%type == 'ELECTRON') lat%ele(i)%value(charge$) = 0
        if (lat%ele(i)%type == 'POSITRON') lat%ele(i)%value(charge$) = -1
      end do
      !    print *
      !    print *, ' positrons '
    endif
    
    if (ipretz == 2) then 
      lat%param%particle = electron$
      co(0)%vec(:) = co_electron(0)%vec(:)
      do i = 1,lat%n_ele_track
        if (lat%ele(i)%key /= beambeam$) cycle
        if (lat%ele(i)%type == 'ELECTRON') lat%ele(i)%value(charge$) = -1
        if (lat%ele(i)%type == 'POSITRON') lat%ele(i)%value(charge$) = 0 
      end do
      !   print *
      !   print *,' electrons '
    endif
    
    !  co(0)%vec = 0.
    co(0)%species = lat%param%particle

    call clear_lat_1turn_mats( lat )
    if (con%circular_machine) call closed_orbit_calc( lat, co, 4 )
    
    call track_all (lat, co)
    
    call lat_make_mat6(lat,-1,co)
    lat%ele(0) = lat%ele_init !for not closed
    if (con%circular_machine) call twiss_at_start(lat, status = status)
    if (status /= ok$) then
      print *,' ELECTRON_POSITRON calling TWISS_AT_START:    ' 
      if (lat%param%particle == electron$) print *,'   ELECTRON optics unstable'
      if (lat%param%particle == positron$) print *,'   POSITRON optics unstable'
      do i=0,lat%n_ele_track
        lat%ele(i)%a%beta = 0.
        lat%ele(i)%b%beta = 0.
      end do
    else
      call twiss_propagate_all (lat)
    endif
    
    if (lat%param%particle == positron$) then
      lat_pos = lat
      forall (i=0:lat%n_ele_track) 
        co_positron(i) = co(i)
      end forall
      if ((any(con%c(1:con%n_constraint)%variable == d_xemit$) .or. &
           any(con%c(1:con%n_constraint)%variable == d_curly_d$)) .and. &
           lat%ele(0)%a%beta /= 0.) &
           call radiation_integrals (lat, co, global%mode_positron, ix_positron)
      if (debug) then
        print *,' positron:', lat%ele(0)%a%beta
        print *,' lat%ele(0)%a%beta ', lat%ele(0)%a%beta
      endif
    
    endif

    if (lat%param%particle == electron$) then
      lat_ele = lat
      forall (i=0:lat%n_ele_track) 
        co_electron(i) = co(i)
      end forall
      if ((any(con%c(1:con%n_constraint)%variable == d_xemit$) .or. &
           any(con%c(1:con%n_constraint)%variable == d_curly_d$)) .and. &
           lat%ele(0)%a%beta /= 0.) &
           call radiation_integrals (lat, co, global%mode_electron, ix_electron) 
      if (debug) then
        print *,' electron:'
        print *,' lat%ele(0)%a%beta ', lat%ele(0)%a%beta
      endif
    endif
        
  end do
  
  return
end subroutine electron_positron