anchor.diff:        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)
  +      write(lout,"(a)")      "CHEBY> ERROR R2 cannot be larger than domain of Chebyshev polynomials!"
  +      goto 10 
       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
  +   
  +    ! checks on maps
  +    if ( cheby_lMap(jj) ) then
  +      if (cheby_mapXmax(jj)<cheby_mapXmin(jj)) then
  +        ! swap
  +        tmpFlt=cheby_mapXmax(jj)
  +        cheby_mapXmax(jj)=cheby_mapXmin(jj)
  +        cheby_mapXmin(jj)=tmpFlt
  +      end if
  +      if (cheby_mapXmax(jj)==cheby_mapXmin(jj)) then
  +        write(lout,"(a)")                     "CHEBY> ERROR X-extremes for map coincide!"
  +        write(lout,"(a,1pe22.15,a,1pe22.15)") "CHEBY>       xmin [mm]: ",cheby_mapXmin(jj), &
  +                                                        " - xmax [mm]: ",cheby_mapXmax(jj)
  +        goto 10 
  +      end if
  +      if (cheby_mapNx(jj)<1) then
  +        write(lout,"(a)")    "CHEBY> ERROR wrong X-stepping for map!"
  +        write(lout,"(a,i0)") "CHEBY>       must be >0 - got: ",cheby_mapNx(jj)
  +        goto 10 
  +      end if
  +      if (cheby_mapYmax(jj)<cheby_mapYmin(jj)) then
  +        ! swap
  +        tmpFlt=cheby_mapYmax(jj)
  +        cheby_mapYmax(jj)=cheby_mapYmin(jj)
  +        cheby_mapYmin(jj)=tmpFlt
  +      end if
  +      if (cheby_mapYmax(jj)==cheby_mapYmin(jj)) then
  +        write(lout,"(a)")                     "CHEBY> ERROR Y-extremes for map coincide!"
  +        write(lout,"(a,1pe22.15,a,1pe22.15)") "CHEBY>       ymin [mm]: ",cheby_mapYmin(jj), &
  +                                                        " - ymax [mm]: ",cheby_mapYmax(jj)
  +        goto 10 
  +      end if
  +      if (cheby_mapNy(jj)<1) then
  +        write(lout,"(a)")    "CHEBY> ERROR wrong Y-stepping for map!"
  +        write(lout,"(a,i0)") "CHEBY>       must be >0 - got: ",cheby_mapNy(jj)
  +        goto 10 
  +      end if
  +    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)
  +      write(lout,"(a,i0,a)") 'CHEBY> status of chebyshev lens #',jj," - name: '"//trim(bez(kk))//"'"
  +      write(lout,"(a)")          "CHEBY> - filename         : '"//trim(cheby_filename(cheby_itable(jj)))//"'"
  +      write(lout,"(a,1pe22.15)") "CHEBY> - R2           [mm]: ",cheby_r2(jj)
  +      write(lout,"(a,1pe22.15)") "CHEBY> - R1           [mm]: ",cheby_r1(jj)
  +      write(lout,"(a,1pe22.15)") "CHEBY> - tilt angle  [deg]: ",cheby_angle(jj)
  +      write(lout,"(a,1pe22.15)") "CHEBY> - hor offset   [mm]: ",cheby_offset_x(jj)
  +      write(lout,"(a,1pe22.15)") "CHEBY> - ver offset   [mm]: ",cheby_offset_y(jj)
  +      write(lout,"(a,1pe22.15)") "CHEBY> - lens powering [A]: ",cheby_I(jj)
  +      write(lout,"(a,1pe22.15)") "CHEBY> - scaling factor []: ",cheby_scalingFact(jj)
  +      if ( cheby_lMap(jj) ) then
  +        write(lout,"(a)") 'CHEBY> requested dump of potential map'
  +        write(lout,"(a,a)")        "CHEBY> - map filename : ",cheby_mapFileName(jj)
  +        write(lout,"(a,1pe22.15)") "CHEBY> - xmin     [mm]: ",cheby_mapXmin(jj)
  +        write(lout,"(a,1pe22.15)") "CHEBY> - xmax     [mm]: ",cheby_mapXmax(jj)
  +        write(lout,"(a,1pe22.15)") "CHEBY> - Nx         []: ",cheby_mapNx  (jj)
  +        write(lout,"(a,1pe22.15)") "CHEBY> - ymin     [mm]: ",cheby_mapYmin(jj)
  +        write(lout,"(a,1pe22.15)") "CHEBY> - ymin     [mm]: ",cheby_mapYmax(jj)
  +        write(lout,"(a,1pe22.15)") "CHEBY> - Ny         []: ",cheby_mapNy  (jj)
  +        call cheby_potentialMap(jj,kk)
  +      end if
       end if
     end do
  -
  +  return
  +
  +10 continue
  +   write(lout,"(a,i0,a)") "CHEBY>       concerned lens #", jj," - name: '"//trim(bez(kk))//"'"
  +   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 subroutine cheby_postInput
   
   
  
positive.code:   integer, intent(in) :: ix
    integer, intent(in) :: n
    
    real(kind=fPrec) xx, yy, rr, frr, dxp, dyp
    real(kind=fPrec) theta, radio, angle_rad
    integer          j
    logical          lrotate
  
    ! rotation angle
    lrotate = cheby_angle(icheby(ix)).ne.zero
    angle_rad = (cheby_angle(icheby(ix))/c180e0)*pi
  
    do j=1,napx
  
      ! apply offset
      xx=xv1(j)-cheby_offset_x(icheby(ix))
      yy=xv2(j)-cheby_offset_y(icheby(ix))
  
      ! check that particle is within the domain of chebyshev polynomials
      rr=sqrt(xx**2+yy**2)
  
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 )
  