anchor.diff:    if ( abs(RRA)>epsilon**2 ) then
   !FOX  RR=SQRT(RR_SQ) ;
     else
  +!FOX  RR_SQ=ZERO ;
   !FOX  RR=ZERO ;
     end if
     RRA=zero
  +  call dapek(RR,hh,RRA)
  +  if (st_debug) then
  +    write(lout,'(a,1pe23.16)')   'ELENS> ELENS_KICK_FOX computing at RRA=',RRA
  +    write(lout,'(2(a,1pe23.16))')'ELENS>                when R1=',elens_r1(iLens),' and R2=',elens_r2(iLens)
  +    flush(lout)
  +  end if
   
     ! 3) calculate kick
     !    shape function: spatial charge density depends on type:
     !    0    if r <= R1
     !    frr  if R1 < r < R2
     !    1    if r >= R2
  -  call dapek(RR,hh,RRA)
  -  if (st_debug) then
  -    write(lout,'(a,1pe23.16)')   'ELENS> ELENS_KICK_FOX computing at RRA=',RRA
  -    write(lout,'(2(a,1pe23.16))')'ELENS>                when R1=',elens_r1(iLens),' and R2=',elens_r2(iLens)
  -  end if
  -    
  -  if ( RRA>elens_r1(iLens)+epsilon .or. &
  -       (elens_lZeroThick(iLens).and.abs(RRA-elens_r1(iLens))<epsilon) ) then ! rr<=R1 -> no kick from elens
  -    if (RRA<elens_r2(iLens)-epsilon) then ! R1<rr<R2
  -      flush(lout)
  +  if ( RRA > elens_r1(iLens)*lteps ) then ! rr<R1: no kick from elens
  +    if ( RRA < elens_r2(ielens(ix))*lteps ) then ! R1<=rr<R2: type-dependent kick
  +       
         select case (elens_type(iLens))
  -      case (1)
  -        ! UNIFORM: eLens with uniform profile
  -        ! formula: (r^2-r1^2)/(r2^2-r1^2)
  -!FOX    FRR=RR_SQ-ER1_SQ ;
  +        
  +      case (1) ! UNIFORM: eLens with uniform profile
  +        if ( elens_lFull(ielens(ix)) ) then
  +          ! formula: r/r2
  +!FOX      FRR=RR/ELE_R2 ;
  +        else
  +          ! formula: (r^2-r1^2)/(r2^2-r1^2)*r2/r
  +!FOX      FRR=(((RR+ELE_R1)*(RR-ELE_R1))/ELENOR)*(ELE_R2/RR) ;
  +        end if
            
         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))
  -!FOX    FRR=EXP(-HALF*(ER1_SQ/ESG_SQ))-EXP(-HALF*(RR_SQ/ESG_SQ)) ;
  -        
  +        if ( elens_lFull(ielens(ix)) ) then
  +          if ( RRA <= c1m7 ) then
  +            ! formula: r*r2/2sig^2/(exp(-r1^2/2sig^2)-exp(-r2^2/2sig^2))
  +!FOX        FRR=(((RR/ELESIG)*(ELE_R2/ELESIG))/TWO)/ELENOR ;
  +          else
  +            ! formula: (1-exp(-r^2/2sig^2))/(exp(-r1^2/2sig^2)-exp(-r2^2/2sig^2))*r2/r
  +!FOX        FRR=((ONE-EXP(-HALF*((RR/ELESIG)*(RR/ELESIG))))/ELENOR)*(ELE_R2/RR) ;
  +          end if     
  +        else
  +          ! formula: (exp(-r1^2/2sig^2)-exp(-r^2/2sig^2))/(exp(-r1^2/2sig^2)-exp(-r2^2/2sig^2))*r2/r
  +!FOX      FRR=((EXP(-HALF*((ELE_R1/ELESIG)*(ELE_R1/ELESIG)))-EXP(-HALF*((RR/ELESIG)*(RR/ELESIG))))/ELENOR)*(ELE_R2/RR) ;
  +        end if
  +         
         case (3)
           ! RADIAL PROFILE: eLens with radial profile as from file
           ! formula: (cumul_J(r)-cumul_J(r1))/(cumul_J(r2)-cumul_J(r1))
  
positive.code: !FOX      TMPRR=TMPRR*RR ;
          end do
          call dealloc(cof,'cof')
        case default
          write(lerr,"(a,i0,a)") "ELENS> ERROR in elens_kick_fox: 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)
  !FOX    FRR=FRR/ELENOR ;
      else ! rr>=R2
  !FOX    FRR=ONE ;
      endif
      ! 'radial kick'
  !FOX    FRR=(((ELETR2/RR)*ELE_R2)*FRR)*(MTCDA/(ONE+DPDA)) ;
      if(elens_lThetaR2(ielens(ix))) then
        if(elens_I(ielens(ix))/elens_beam_chrg(ielens(ix)) < zero) then
  !FOX    FRR=FRR*((RV+ELEBET*BETA0)/(ONE+ELEBET*BETA0)) ;
  
negative.code:     !    radial position of main beam relative to center of elens beam
      !    rr = sqrt(xx**2+yy**2)
      rr=sqrt((xx+yy)*(xx+yy)-two*(xx*yy))
      if (abs(rr)<epsilon) rr=zero
      
      ! 3) calculate kick
      !    shape function: spatial charge density depends on type:
  