!+           
! Subroutine Implement_Ecloud(ring,sigx,sigy,sigz,density, bunch, ecloud_type, beamfield_file, bbi_length, d_bbi) 
!
! Insert a beambeam element at each bend as a electron cloud 
!
! Input:
!   ring -- lat_struct: New elements will be inserted
!   sig_x -- real: cloud width
!   sig_y -- real: cloud height
!   sig_z -- real: bunch length
!   density -- real: density of the cloud
!   bbi_length -- real: parameter that defines longitudinal dependence
! Output:
!   ring -- lat_struct: ring with new elements inserted
!

  subroutine implement_ecloud(ring,sigx,sigy,sigz,density,bunch, ecloud_type, beamfield_file, bbi_length, d_bbi) 
  use bmad
  use bookkeeper_mod

 implicit none

 type (lat_struct) ring
 type (ele_struct), allocatable :: ecloud_ele(:)

 integer i,j
 integer j_max
 integer ix, ii
 integer nslice/1/
 integer bunch

 real(rp) sigx,sigy,sigz
  real(rp), save :: sigx_s,sigy_s,sigz_s
 real(rp) delta_s
 real(rp) density
 real(rp), allocatable, save :: length(:)
 real(rp) bbi_length, d_bbi

 character*29 beamfield_file
 character*3 ecloud_type
 character*2 number

 logical ecloud/.false./


 ring%param%n_part = sigx*sigy*sigz * density
 
 if( .not. ecloud)then
  j=0
  do i=1,ring%n_ele_track
   if(ring%ele(i)%key == sbend$ .or. ring%ele(i)%key == rbend$) then
    j=j+1
   endif
  end do

  j_max = j
  allocate (ecloud_ele(1:j), length(1:j))

  j=0
  number = '00'
  do i=1,ring%n_ele_track
   if(ring%ele(i)%key == sbend$ .or. ring%ele(i)%key == rbend$) then
    j=j+1
    length(j) = ring%ele(i)%value(l$)
    if(j<10)write(number(2:2),'(i1)')j
    if(j>=10)write(number(1:2),'(i2)')j
    if(j>99) print *,' Too many bends for implement beambeam'
    call init_ele(ecloud_ele(j))
!     print '(a,i10,1x,a)', ' j, beambeam_ele(j)%name ', j, beambeam_ele(j)%name
      ecloud_ele(j)%key = custom$
    ecloud_ele(j)%s = ring%ele(i-1)%s
    ecloud_ele(j)%tracking_method = custom$
    ecloud_ele(j)%mat6_calc_method = tracking$
    ecloud_ele(j)%value(p0c_start$)=ring%ele(0)%value(p0c_start$)
    ecloud_ele(j)%value(e_tot_start$)=ring%ele(0)%value(e_tot_start$)
    ecloud_ele(j)%value(p0c$)=ring%ele(0)%value(p0c$)
    ecloud_ele(j)%value(e_tot$)=ring%ele(0)%value(e_tot$)

    if(ecloud_type == 'BBI')then
      ecloud_ele(j)%name = 'ECLOUD_BBI'//number//ring%ele(i)%name(1:len_trim(ring%ele(i)%name))
      ecloud_ele(j)%value(sig_x$) = sigx
      ecloud_ele(j)%value(sig_y$) = sigy
      ecloud_ele(j)%value(sig_z$) = sigz
      ecloud_ele(j)%value(n_slice$) = nslice
      ecloud_ele(j)%value(charge$) = -1
      !!! ecloud_ele(j)%value(scratch1$) = bbi_length
      !!! ecloud_ele(j)%value(scratch2$) = d_bbi
      ecloud_ele(j)%value(bbi_const$) =  -ring%param%n_part*(length(j)/sigz) * ecloud_ele(j)%value(charge$) * classical_radius_factor /  &
                                          (2 * pi * ecloud_ele(j)%value(p0c$) * (sigx + sigy))
    endif
    if(ecloud_type == 'MAP')then
      ecloud_ele(j)%name = 'ECLOUD_MAP'//number//ring%ele(i)%name(1:len_trim(ring%ele(i)%name))
      !!! ecloud_ele(j)%value(scratch1$) = bunch
      !!! ecloud_ele(j)%value(scratch2$) = length(j)
      ecloud_ele(j)%alias = beamfield_file
    endif      
   endif
  end do
 
  do j = 1,j_max
    ix = element_at_s(ring,ecloud_ele(j)%s, .true.)
!   print *,j,rf_ele(j),ix, ring%ele(ix)%name
    call insert_element(ring, ecloud_ele(j),ix)
    if (.not. bmad_com%auto_bookkeeper) call lattice_bookkeeper(ring)
    call lat_make_mat6(ring,-1)
    print *,' ring%n_ele_max, ix = ', ring%n_ele_max,ix
    print '(i10,1x,2a16)',ix,ring%ele(ix)%name, ring%ele(ix+1)%name

  end do

 ecloud =.true.

 print *,' ECLOUD IMPLEMENTED'

 print *,' this is the real deal'
!deallocate (ecloud_ele)
endif ! if not ecloud



  return
  end subroutine