anchor.diff:  
     ! variables to save parameters for tracking etc.
     integer, save          :: cheby_itable(ncheby)=0         ! index of chebyshev table
  -  real(kind=fPrec), save :: cheby_offset_x(ncheby)=zero    ! hor. offset [mm]
  -  real(kind=fPrec), save :: cheby_offset_y(ncheby)=zero    ! ver. offset [mm]
  -  real(kind=fPrec), save :: cheby_angle(ncheby)=zero       ! rotation angle about the longitudinal axis [rad]
  -  real(kind=fPrec), save :: cheby_scalingFact(ncheby)=zero ! scaling factor []
  +  real(kind=fPrec), save :: cheby_r2(ncheby)     = zero    ! outer radius R2 [mm] (optional)
  +  real(kind=fPrec), save :: cheby_r1(ncheby)     = zero    ! inner radius R1 [mm] (optional)
  +  real(kind=fPrec), save :: cheby_angle(ncheby)=zero       ! rotation angle about the longitudinal axis [rad] (optional)
  +  real(kind=fPrec), save :: cheby_offset_x(ncheby)=zero    ! hor. offset [mm] (optional)
  +  real(kind=fPrec), save :: cheby_offset_y(ncheby)=zero    ! ver. offset [mm] (optional)
  +  real(kind=fPrec), save :: cheby_L(ncheby)=-one           ! actual length of (thick) lens [m] (optional)
  +  real(kind=fPrec), save :: cheby_I(ncheby)=-one           ! actual powering of lens [A] (optional)
  +  real(kind=fPrec), save :: cheby_scalingFact(ncheby)=one  ! scaling factor [] (computed internally)
   
     ! tables with chebyshev coefficients
     integer, parameter     :: cheby_max_order=30             ! max order of chebyshev polynomials currently supported
     character(len=mFNameLen), save:: cheby_filename(ncheby_tables) = " "! file names
     real(kind=fPrec), save :: cheby_coeffs(0:cheby_max_order,0:cheby_max_order,ncheby_tables) = zero ! coefficients
     integer, save          :: cheby_maxOrder(ncheby_tables)  = 0    ! max order of the current map
  -  real(kind=fPrec), save :: cheby_refCurr(ncheby_tables)   = zero ! reference current [A]
  -  real(kind=fPrec), save :: cheby_refRadius(ncheby_tables) = zero ! reference radius [mm]
  -  real(kind=fPrec), save :: cheby_refBeta(ncheby_tables)   = zero ! reference e-beta []
  +  real(kind=fPrec), save :: cheby_refI(ncheby_tables) = one  ! reference current [A] (optional)
  +  real(kind=fPrec), save :: cheby_refR(ncheby_tables) = zero ! reference radius [mm] (mandatory)
  +  real(kind=fPrec), save :: cheby_refL(ncheby_tables) = one  ! reference length [m]  (optional)
   
   contains
   
  
positive.code: 
  
  subroutine cheby_postInput
    use mathlib_bouncer
    use utils
  
    integer j
    logical exist
    
    ! Parse files with coefficients for Chebyshev polynomials:
     do j=1,mcheby
      inquire(file=cheby_filename(j), exist=exist)
      if(.not. exist) then
        write(lout,"(a)") "CHEBY> ERROR Problems with file with coefficients for Chebyshev polynominals: ", &
              trim(cheby_filename(j))
        call prror(-1)
      end if
      call parseChebyFile(j)
    end do
  
  end subroutine cheby_postInput
  
  
  subroutine cheby_kick(jcheby)
  
    ! A. Mereghetti (CERN, BE-ABP-HSS)
    ! last modified: 20-02-2018
    ! apply kick of electron lens
  
    use crcoall
    use mod_common
    use mod_common_main
    use mathlib_bouncer
    use numerical_constants
    use physical_constants
  
    real(kind=fPrec) :: xx, yy, rr, frr, dxp, dyp
    real(kind=fPrec) :: theta, radio, angle_rad, brho_b, beta_b
    integer          :: j, jcheby
    logical          :: lrotate
  
    angle_rad = zero ! -Wmaybe-uninitialized
  
    ! rotation angle
    lrotate=cheby_angle(jcheby).ne.zero
  
    ! Brho of beam
    beta_b=e0f*gammar/pma
    brho_b=e0f/(clight*c1m6)
  
    do j=1,napx
  
      ! apply offset
      xx=xv1(j)-cheby_offset_x(jcheby)
      yy=xv2(j)-cheby_offset_y(jcheby)
  
      ! check that particle is within the domain of chebyshev polynomials
      rr=sqrt(xx**2+yy**2)
      if (rr.gt.cheby_refRadius(cheby_itable(jcheby))) then
        write(lout,"(a,3(e12.6,1x),a,e12.6)") "CHEBY> ERROR in cheby_kick: particle at position (x,y,r): ",&
          xv1(j), xv2(j), rr,' is outside radial domain of Chebyshev polinomials: ',cheby_refRadius(cheby_itable(jcheby))
        call prror
      end if
  
      ! in case of non-zero tilt angle, rotate coordinates
      if (lrotate) then
        theta = atan2_mb(yy, xx)-angle_rad
        xx = rr * cos_mb(theta)
        yy = rr * sin_mb(theta)
      end if
  
      ! apply kick
      call cheby_getKick( xx, yy, dxp, dyp, cheby_itable(jcheby), brho_b, beta_b )
      ! take into account scaling factor
      dxp=dxp *cheby_scalingFact(jcheby)
      dyp=dyp *cheby_scalingFact(jcheby)
  
      ! in case cheby has a non-zero angle, rotate kicks
      if (lrotate) then
        ! NB: cheby_angle(jcheby) is the rotation angle of the cheby
        theta = atan2_mb(dyp, dxp)+angle_rad
        radio = sqrt(dxp**2 + dyp**2)
        dxp = radio * cos_mb(theta)
        dyp = radio * sin_mb(theta)
      end if
  
      ! apply kicks, taking into account magnetic rigidity of particle being tracked;
      yv1(j)=yv1(j)+dxp *oidpsv(j)
      yv2(j)=yv2(j)+dyp *oidpsv(j)
    end do
  
  end subroutine cheby_kick
  
  
    subroutine cheby_getKick( xx, yy, dxp, dyp, iTable, brho_b, beta_b )
  
      ! A. Mereghetti (CERN, BE-ABP-HSS)
      ! last modified: 28-02-2018
      ! compute kicks from Chebyshev polinomials - see FermiLAB-FN-0972-APC
  
      use mathlib_bouncer
      use physical_constants
      use numerical_constants
  
      ! interface vars
      real(kind=fPrec) :: xx, yy, dxp, dyp, brho_b, beta_b
      integer          :: iTable
  
      ! temp vars
      real(kind=fPrec) :: uu, vv, Tx (0:cheby_maxOrder(iTable)), Ty (0:cheby_maxOrder(iTable)), &
                          kx, ky, Tpx(0:cheby_maxOrder(iTable)), Tpy(0:cheby_maxOrder(iTable)), &
                          fu, fv
      integer          :: ii, jj
  
      ! normalised variables
      uu=xx/cheby_refRadius(iTable)
      vv=yy/cheby_refRadius(iTable)
      ! normalisation factors of derivatives
      fu=(one-uu)*(one+uu)
      fv=(one-vv)*(one+vv)
  
      ! polynomials:
      Tx(0)=one
      Ty(0)=one
      Tx(1)=uu
      Ty(1)=vv
      do ii=2,cheby_max_order
         Tx(ii)=two*uu*Tx(ii-1)-Tx(ii-2)
         Ty(ii)=two*vv*Ty(ii-1)-Ty(ii-2)
      end do
      ! derivatives:
      Tpx(0)=zero
      Tpy(0)=zero
      do ii=1,cheby_max_order
         Tpx(ii)=real(ii,fPrec)*(Tx(ii-1)-uu*Tx(ii))/fu
         Tpy(ii)=real(ii,fPrec)*(Ty(ii-1)-vv*Ty(ii))/fv
      end do
  
      ! get kicks
      dxp=zero
      dyp=zero
      do ii=0,cheby_maxOrder(iTable)
         do jj=0,ii
            dxp=dxp+cheby_coeffs(jj,ii,iTable)*Tpx(jj)*Ty (ii-jj)
            dyp=dyp+cheby_coeffs(jj,ii,iTable)*Tx (jj)*Tpy(ii-jj)
         end do
      end do
      dxp=-dxp/cheby_refRadius(iTable)
      dyp=-dyp/cheby_refRadius(iTable)
  
      ! take into account Brho and beta
      dxp=dxp/(brho_b*clight*beta_b)
      dyp=dyp/(brho_b*clight*beta_b)
  
     end subroutine cheby_getKick
  
  ! ================================================================================================ !
  !  Last modified: 2019-02-25
  
negative.code:   end if
  
    if(inLine(1:1) == "I") then
      ! Read reference current of lens
      if(nSplit < 3) then
        write(lout,"(a)") "CHEBY> ERROR Not enough arguments for expressing ref lens current [A]."
        write(lout,"(a)") "CHEBY>       Correct format:"
        write(lout,"(a)") "I : value"
        goto 30
      end if
      call chr_cast(lnSplit(3),cheby_refI(ifile),spErr)
      lDefI=.false.
  
    else if(inLine(1:1) == "L") then
      ! Read reference length of lens
      if(nSplit < 3) then
        write(lout,"(a)") "CHEBY> ERROR Not enough arguments for expressing ref lens length [m]."
        write(lout,"(a)") "CHEBY>       Correct format:"
        write(lout,"(a)") "L : value"
        goto 30
      end if
      call chr_cast(lnSplit(3),cheby_refL(ifile),spErr)
      lDefL=.false.
  
    else if(inLine(1:1) == "R") then
      ! Read reference radius e-beam in e-lens
      if(nSplit < 3) then
        write(lout,"(a)") "CHEBY> ERROR Not enough arguments for expressing ref lens radius [mm]."
        write(lout,"(a)") "CHEBY>       Correct format:"
        write(lout,"(a)") "R : value"
        goto 30
      end if
      call chr_cast(lnSplit(3),cheby_refR(ifile),spErr)
  
    else
      ! Read chebyshev coefficients
      if(nSplit /= 4) then
        goto 30
      end if
      call chr_cast(lnSplit(1),ii,spErr)
      if(ii > cheby_max_order) then
        goto 30
      end if
      call chr_cast(lnSplit(2),jj,spErr)
      if(jj > cheby_max_order) then
        goto 30
      end if
      call chr_cast(lnSplit(4),cheby_coeffs(ii,jj,ifile),spErr)
  
    end if ! close if for keyword identification
    goto 10
  