anchor.diff:    Tpx(1)=one
     Tpy(1)=one
     do nn=2,cheby_maxOrder(iTable)
  -     Tx(nn)=two*(uu*Tx(nn-1))-Tx(nn-2)
  -     Ty(nn)=two*(vv*Ty(nn-1))-Ty(nn-2)
  -     Tpx(nn)=(real(nn,fPrec)*(Tx(nn-1)-uu*Tx(nn)))/fu
  -     Tpy(nn)=(real(nn,fPrec)*(Ty(nn-1)-vv*Ty(nn)))/fv
  +    Tx(nn)=two*(uu*Tx(nn-1))-Tx(nn-2)
  +    Ty(nn)=two*(vv*Ty(nn-1))-Ty(nn-2)
  +    Tpx(nn)=(real(nn,fPrec)*(Tx(nn-1)-uu*Tx(nn)))/fu
  +    Tpy(nn)=(real(nn,fPrec)*(Ty(nn-1)-vv*Ty(nn)))/fv
     end do
   
     ! get kicks
     dxp=zero
     dyp=zero
     do nn=0,cheby_maxOrder(iTable)
  -     do jj=0,nn
  -        dxp=dxp+(cheby_coeffs(jj,nn-jj,iTable)*Tpx(jj))*Ty (nn-jj)
  -        dyp=dyp+(cheby_coeffs(jj,nn-jj,iTable)*Tx (jj))*Tpy(nn-jj)
  -     end do
  +    do jj=0,nn
  +      dxp=dxp+(cheby_coeffs(jj,nn-jj,iTable)*Tpx(jj))*Ty (nn-jj)
  +      dyp=dyp+(cheby_coeffs(jj,nn-jj,iTable)*Tx (jj))*Tpy(nn-jj)
  +    end do
     end do
     dxp=-(dxp*c1e3)/(cheby_refR(iTable)*c1m3) ! ref radius in [mm], kick in [mrad]
     dyp=-(dyp*c1e3)/(cheby_refR(iTable)*c1m3) ! ref radius in [mm], kick in [mrad]
   
    end subroutine cheby_getKick
   
  +
  +subroutine cheby_getPotential( xx, yy, zz, iTable )
  +
  +  ! A. Mereghetti (CERN, BE-ABP-HSS)
  +  ! last modified: 01-03-2019
  +  ! compute potential from Chebyshev polinomials - see FermiLAB-FN-0972-APC
  +
  +  use numerical_constants, only : zero, one, two
  +
  +  ! interface vars
  +  real(kind=fPrec), intent(in ) :: xx
  +  real(kind=fPrec), intent(in ) :: yy
  +  real(kind=fPrec), intent(out) :: zz
  +  integer,          intent(in ) :: iTable
  +
  +  ! temp vars
  +  real(kind=fPrec) :: uu, vv, Tx(0:cheby_maxOrder(iTable)), Ty(0:cheby_maxOrder(iTable))
  +  integer          :: nn, jj
  +
  +  ! normalised variables
  +  uu=xx/cheby_refR(iTable)
  +  vv=yy/cheby_refR(iTable)
  +
  +  ! - polynomials:
  +  Tx(0)=one
  +  Ty(0)=one
  +  Tx(1)=uu
  +  Ty(1)=vv
  +  do nn=2,cheby_maxOrder(iTable)
  +    Tx(nn)=two*(uu*Tx(nn-1))-Tx(nn-2)
  +    Ty(nn)=two*(vv*Ty(nn-1))-Ty(nn-2)
  +  end do
  +
  +  ! get kicks
  +  zz=zero
  +  do nn=0,cheby_maxOrder(iTable)
  +    do jj=0,nn
  +      zz=zz+(cheby_coeffs(jj,nn-jj,iTable)*Tx(jj))*Ty(nn-jj)
  +    end do
  +  end do
  +
  + end subroutine cheby_getPotential
  +
   end module cheby
  
positive.code:       cheby_r1(jj)=tmpFlt
      end if
      if (cheby_r2(jj)>cheby_refR(cheby_itable(jj))) then
        write(lout,"(a)") "CHEBY> ERROR R2 cannot be larger than domain of Chebyshev polynomials!"
        write(lout,"(a,1pe22.15,a,1pe22.15)") "CHEBY>       R2 [mm]: ",cheby_r2(jj), &
             " - reference radius [mm]:",cheby_refR(cheby_itable(jj))
        call prror(-1)
      end if
      if (cheby_I (jj)<=zero) then
        cheby_I (jj)=cheby_refI(cheby_itable(jj))
      else
        cheby_scalingFact(jj)=cheby_scalingFact(jj)*((cheby_I(jj))/(cheby_refI(cheby_itable(jj))))
      end if
      if(st_quiet < 2) then
        ! find name
        do kk=1,nele
          if(kz(kk) == cheby_kz .and. icheby(kk) == jj ) then
            exit
          end if
        end do
        write(lout,"(a)") ''
        write(lout,"(a,i0,a)") ' status of chebyshev lens #',jj," - name: '"//trim(bez(kk))//"'"
        write(lout,"(a)")          " - filename         : '"//trim(cheby_filename(cheby_itable(jj)))//"'"
        write(lout,"(a,1pe22.15)") " - R2           [mm]: ",cheby_r2(jj)
        write(lout,"(a,1pe22.15)") " - R1           [mm]: ",cheby_r1(jj)
        write(lout,"(a,1pe22.15)") " - tilt angle  [deg]: ",cheby_angle(jj)
        write(lout,"(a,1pe22.15)") " - hor offset   [mm]: ",cheby_offset_x(jj)
        write(lout,"(a,1pe22.15)") " - ver offset   [mm]: ",cheby_offset_y(jj)
        write(lout,"(a,1pe22.15)") " - lens powering [A]: ",cheby_I(jj)
        write(lout,"(a,1pe22.15)") " - scaling factor []: ",cheby_scalingFact(jj)
      end if
    end do
  
  end subroutine cheby_postInput
  
  
  
negative.code:      
        ! take into account scaling factor, Brho of beam and its relativistic beta,
        !    and magnetic rigidity and relativistic beta of particle being tracked
        dxp=((dxp*cheby_scalingFact(icheby(ix)))/(brho*clight*betrel)*moidpsv(jj))*rvv(jj)
        dyp=((dyp*cheby_scalingFact(icheby(ix)))/(brho*clight*betrel)*moidpsv(jj))*rvv(jj)
        
        ! apply kicks
        yv1(jj)=yv1(jj)+dxp
        yv2(jj)=yv2(jj)+dyp
      end if
    end do
  
  end subroutine cheby_kick
  
  
  subroutine cheby_potentialMap(iTable,ix)
  
    ! A. Mereghetti (CERN, BE-ABP-HSS)
    ! last modified: 01-03-2019
    ! dump map of potential
  
    use mod_common, only : bez
    use mod_common_main
    use mathlib_bouncer
    use numerical_constants, only : zero, c180e0, pi
    use mod_units
  
    integer, intent(in) :: iTable
    integer, intent(in) :: ix
    
    real(kind=fPrec) xx, yy, rr, zz, dx, dy
    real(kind=fPrec) theta, radio, angle_rad
    integer          ii, jj, inside, fUnit
    logical          lrotate, err
  
    write(lout,"(a)") "CHEBY> Dumping potential map...'"
    call f_requestUnit(cheby_mapFileName(iTable),fUnit)
    call f_open(unit=fUnit,file=cheby_mapFileName(iTable),mode='w',err=err,formatted=.true.,status="new")
    if(err) then
      write(lout,"(a)") "CHEBY> ERROR Failed to open file."
      call prror(-1)
    end if
   
    ! rotation angle
    lrotate = cheby_angle(iTable).ne.zero
    angle_rad = (cheby_angle(iTable)/c180e0)*pi
    dx=(cheby_mapXmax(iTable)-cheby_mapXmin(iTable))/real(cheby_mapNx(iTable),fPrec)
    dy=(cheby_mapYmax(iTable)-cheby_mapYmin(iTable))/real(cheby_mapNy(iTable),fPrec)
  
    ! write header
    write(fUnit,'("# ",a,i0)')       'name of chebyshev lens: "'//trim(bez(ix))//'" - id: ',iTable
    write(fUnit,"('# ',a)")          "filename         : '"//trim(cheby_filename(iTable))//"'"
    write(fUnit,"('# ',a,1pe22.15)") "R2           [mm]: ",cheby_r2(iTable)
    write(fUnit,"('# ',a,1pe22.15)") "R1           [mm]: ",cheby_r1(iTable)
    write(fUnit,"('# ',a,1pe22.15)") "tilt angle  [deg]: ",cheby_angle(iTable)
    write(fUnit,"('# ',a,1pe22.15)") "hor offset   [mm]: ",cheby_offset_x(iTable)
    write(fUnit,"('# ',a,1pe22.15)") "ver offset   [mm]: ",cheby_offset_y(iTable)
    write(fUnit,"('# ',a,1pe22.15)") "lens powering [A]: ",cheby_I(iTable)
    write(fUnit,"('# ',a,1pe22.15)") "scaling factor []: ",cheby_scalingFact(iTable)
    write(fUnit,"('# ',a,a)")        "map filename : ",cheby_mapFileName(iTable)
    write(fUnit,"('# ',a,1pe22.15)") "xmin [mm]: ",cheby_mapXmin(iTable)
    write(fUnit,"('# ',a,1pe22.15)") "xmax [mm]: ",cheby_mapXmax(iTable)
    write(fUnit,"('# ',a,1pe22.15)") "Nx     []: ",cheby_mapNx  (iTable)
    write(fUnit,"('# ',a,1pe22.15)") "ymin [mm]: ",cheby_mapYmin(iTable)
    write(fUnit,"('# ',a,1pe22.15)") "ymin [mm]: ",cheby_mapYmax(iTable)
    write(fUnit,"('# ',a,1pe22.15)") "Ny     []: ",cheby_mapNy  (iTable)
    write(fUnit,"('# ',a)") "x [mm], y [mm], V [V], inside [0:False,1:True]"
    
    ! get map
    do jj=0,cheby_mapNy(iTable)
      yy=cheby_mapYmin(iTable)+jj*dy-cheby_offset_y(iTable)
       
      do ii=0,cheby_mapNx(iTable)
        xx=cheby_mapXmin(iTable)+ii*dx-cheby_offset_x(iTable)
  
        ! check that particle is within the domain of chebyshev polynomials
        rr=sqrt(xx**2+yy**2)
        if (rr.ge.cheby_r1(iTable).and.rr.lt.cheby_r2(iTable)) then
          inside=1
        else
          ! rr<r1 || rr>=r2 -> no kick from lens
          inside=0
        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
        ! compute kick from cheby map
        call cheby_getPotential( xx, yy, zz, iTable )
        write(fUnit,'(3(1X,1pe22.15),1X,i0)') xx, yy, zz, inside
          
      end do
    end do
  
    call f_close(fUnit)
    
  end subroutine cheby_potentialMap
  
  
  subroutine cheby_getKick( xx, yy, dxp, dyp, iTable )
  