anchor.diff:      if ( (ax < r1lteps .and. ay < r1lteps) .or. (ax > r2gteps .or. ay > r2gteps) ) cycle
   
       ! compute kick from cheby map
  -    call cheby_getKick( xx, yy, dxp, dyp, cheby_itable(icheby(ix)) )
  +    ! - normalised variables
  +    uu=xx/cheby_refR(cheby_itable(icheby(ix)))
  +    vv=yy/cheby_refR(cheby_itable(icheby(ix)))
  +    ! - polynomials:
  +    Tx(0)=one
  +    Ty(0)=one
  +    Tx(1)=uu
  +    Ty(1)=vv
  +    ! - derivatives:
  +    Tpx(0)=zero
  +    Tpy(0)=zero
  +    Tpx(1)=one
  +    Tpy(1)=one
  +    do nn=2,cheby_maxOrder(cheby_itable(icheby(ix)))
  +      Tx(nn)=two*(uu*Tx(nn-1))-Tx(nn-2)
  +      Ty(nn)=two*(vv*Ty(nn-1))-Ty(nn-2)
  +      ! use definition of derivative applied to definition of polynomial:
  +      ! d/du Tn(u)= d/du[ 2u*T_(n-1)(u)-T_(n-2)(u) ] = 2T_(n-1)(u) +2u*T'_(n-1)(u) -T'_(n-2)(u)
  +      !           = 2(T_(n-1)(u) +u*T'_(n-1)(u))-T'_(n-2)(u)
  +      Tpx(nn)=two*(Tx(nn-1)+uu*Tpx(nn-1))-Tpx(nn-2)
  +      Tpy(nn)=two*(Ty(nn-1)+vv*Tpy(nn-1))-Tpy(nn-2)
  +    end do
  +    ! - get actual kicks
  +    dxp=zero
  +    dyp=zero
  +    do nn=0,cheby_maxOrder(cheby_itable(icheby(ix)))
  +      do jj=0,nn
  +        if (cheby_coeffs(jj,nn-jj,cheby_itable(icheby(ix))) /= zero ) then
  +          dxp=dxp+(cheby_coeffs(jj,nn-jj,cheby_itable(icheby(ix)))*Tpx(jj))*Ty (nn-jj)
  +          dyp=dyp+(cheby_coeffs(jj,nn-jj,cheby_itable(icheby(ix)))*Tx (jj))*Tpy(nn-jj)
  +        end if
  +      end do
  +    end do
  +    dxp=-(dxp*c1e3)/(cheby_refR(cheby_itable(icheby(ix)))*c1m3) ! ref radius in [mm], kick in [mrad]
  +    dyp=-(dyp*c1e3)/(cheby_refR(cheby_itable(icheby(ix)))*c1m3) ! ref radius in [mm], kick in [mrad]
  +    ! - check kicks are not NULL
       if (abs(dxp)<epsilon) dxp=zero
       if (abs(dyp)<epsilon) dyp=zero
       if ( dxp == zero .and. dyp == zero ) cycle
  
positive.code: 
    ! A. Mereghetti (CERN, BE-ABP-HSS)
    ! last modified: 02-04-2020
    ! apply kick of Chebyshev lenses
  
    use mod_common, only : beta0, napx, brho
    use mod_common_main
    use mathlib_bouncer
    use numerical_constants, only : zero, c180e0, pi, two, c1m15
    use physical_constants, only: clight
  
    implicit none
  
negative.code:   integer, intent(in) :: ix
    integer, intent(in) :: n
  
    real(kind=fPrec) xx, yy, ax, ay, rr, dxp, dyp, uu, vv
    real(kind=fPrec) theta, angle_rad, epsilon, gteps, lteps, r1lteps, r2gteps
    integer          ii, jj, nn
    logical          lrotate
  
    real(kind=fPrec) :: Tx (0:cheby_maxOrder(cheby_itable(icheby(ix)))), &
                        Ty (0:cheby_maxOrder(cheby_itable(icheby(ix)))), &
                        Tpx(0:cheby_maxOrder(cheby_itable(icheby(ix)))), &
                        Tpy(0:cheby_maxOrder(cheby_itable(icheby(ix))))
    
    epsilon=c1m15
    gteps=one+epsilon
    lteps=one-epsilon
  