anchor.diff:        tmpflt=elens_r2(ielens(iElem))
         elens_r2(ielens(iElem)) = elens_r1(ielens(iElem))
         elens_r1(ielens(iElem)) = tmpflt
  -    else if(elens_r2(ielens(iElem)) == elens_r1(ielens(iElem))) then
  -      write(lerr,"(a)") "ELENS> ERROR ELEN R2=R1. Elens does not exist."
  -      iErr = .true.
  -      return
  -    end if
  -    if(elens_r2(ielens(iElem)) <= zero) then
  -      write(lerr,"(a)") "ELENS> ERROR R2<=0!"
  -      iErr = .true.
  -      return
  +    end if
  +    if(elens_r2(ielens(iElem)) < zero) then
  +      write(lerr,"(a)") "ELENS> ERROR R2<0!"
  +      iErr = .true.
  +      return
  +    else if(elens_r2(ielens(iElem)) < pieni) then
  +      elens_r2(ielens(iElem))=zero
       end if
       if(elens_r1(ielens(iElem)) < zero) then
         write(lerr,"(a)") "ELENS> ERROR R1<0!"
         iErr = .true.
         return
  +    else if(elens_r1(ielens(iElem)) < pieni) then
  +      elens_r1(ielens(iElem))=zero
       end if
       if(elens_lThetaR2(ielens(iElem))) then
         if(elens_len(ielens(iElem)) <= zero) then
  
positive.code:          iErr = .true.
           return
        end if
      end if
      
      if(st_debug) then
  
negative.code:     ! 2) calculate radius
      !    radial position of main beam relative to center of elens beam
      !    rr = sqrt(xx**2+yy**2)
      rr_sq=(xx+yy)*(xx+yy)-two*(xx*yy)
      rr=sqrt(rr_sq)
      ! 3) calculate kick
      !    shape function: spatial charge density depends on type:
      !    0    if r <= R1
      !    frr  if R1 < r < R2
      !    1    if r >= R2
      if (rr.gt.elens_r1(ielens(ix)).or.(elens_lZeroThick(ielens(ix)).and.rr.eq.elens_r1(ielens(ix)))) then ! rr<=R1 -> no kick from elens
        if (rr.lt.elens_r2(ielens(ix))) then ! R1<rr<R2
          select case (elens_type(ielens(ix)))
          case (1)
            ! UNIFORM: eLens with uniform profile
            ! formula: (r^2-r1^2)/(r2^2-r1^2)
            frr=rr_sq-r1_sq
          case (2)
            ! GAUSSIAN: eLens with Gaussian profile
            ! formula: (exp(-r1^2/2sig^2)-exp(-r^2/2sig^2))/(exp(-r1^2/2sig^2)-exp(-r2^2/2sig^2))
            frr=exp_mb(-((r1_sq/elSig_sq)/two))-exp_mb(-((rr_sq/elSig_sq)/two))
          case (3)
            ! RADIAL PROFILE: eLens with radial profile as from file
            ! formula: (cumul_J(r)-cumul_J(r1))/(cumul_J(r2)-cumul_J(r1))
            frr=polinterp( rr, &
                  elens_radial_profile_R(0:elens_radial_profile_nPoints(elens_iRadial(ielens(ix))),elens_iRadial(ielens(ix))), &
                  elens_radial_profile_J(0:elens_radial_profile_nPoints(elens_iRadial(ielens(ix))),elens_iRadial(ielens(ix))), &
                  elens_radial_profile_nPoints(elens_iRadial(ielens(ix)))+1, &
                  elens_radial_mpoints(ielens(ix)), elens_radial_jguess(ielens(ix)) )-elens_radial_fr1(ielens(ix))
          case default
            write(lerr,"(a,i0,a)") "ELENS> ERROR elens_kick: elens_type=",elens_type(ielens(ix))," not recognized. "
            write(lerr,"(a)")      "ELENS>       Possible values for type are: 1, 2 and 3"
            call prror
          end select
          ! take into account normalisation factor (geometrical)
          frr=frr/elens_geo_norm(ielens(ix))
        else ! rr>=R2
          frr=one
        endif
        ! 'radial kick'
        frr = (elens_theta_r2(ielens(ix))*((elens_r2(ielens(ix))/rr)*frr))*moidpsv(jj)
        if(elens_lThetaR2(ielens(ix))) then
          if(elens_I(ielens(ix)) < zero) then
            frr = frr*((rvv(jj)+elens_beta_e(ielens(ix))*beta0)/(one+elens_beta_e(ielens(ix))*beta0))
  