anchor.diff:  #endif
     end if
   
  -  ! Initialise efficiency array
  -  do i=1,iu
  -    sum_ax(i)   = zero
  -    sqsum_ax(i) = zero
  -    sum_ay(i)   = zero
  -    sqsum_ay(i) = zero
  -    nampl(i)    = zero
  -  end do
  -
  -  nspx = zero
  -  nspy = zero
     iturn = 1
     ie    = 1
     n_tot_absorbed = 0
  
positive.code: #endif
  
    if(firstrun) then
  
      if(iturn == 1) then
        sum_ax(ie) = zero
        sum_ay(ie) = zero
      end if
  
      do j = 1, napx
        xj  = (xv1(j)-torbx(ie)) /c1e3
        xpj = (yv1(j)-torbxp(ie))/c1e3
        yj  = (xv2(j)-torby(ie)) /c1e3
        ypj = (yv2(j)-torbyp(ie))/c1e3
        pj  = ejv(j)/c1e3
  
        if(tbetax(ie) > zero) then
          gammax = (one + talphax(ie)**2)/tbetax(ie)
          gammay = (one + talphay(ie)**2)/tbetay(ie)
        else
          gammax = (one + talphax(ie-1)**2)/tbetax(ie-1)
          gammay = (one + talphay(ie-1)**2)/tbetay(ie-1)
        endif
  
        if(part_abs_pos(j) == 0 .and. part_abs_turn(j) == 0) then
          if(tbetax(ie) > 0.) then
            nspx = sqrt(abs( gammax*(xj)**2 + two*talphax(ie)*xj*xpj +   tbetax(ie)*xpj**2 )/c_emitx0_collgap)
            nspy = sqrt(abs( gammay*(yj)**2 + two*talphay(ie)*yj*ypj +   tbetay(ie)*ypj**2 )/c_emity0_collgap)
          else
            nspx = sqrt(abs( gammax*(xj)**2 + two*talphax(ie-1)*xj*xpj + tbetax(ie-1)*xpj**2 )/c_emitx0_collgap)
            nspy = sqrt(abs( gammay*(yj)**2 + two*talphay(ie-1)*yj*ypj + tbetay(ie-1)*ypj**2 )/c_emity0_collgap)
          end if
  
          sum_ax(ie)   = sum_ax(ie) + nspx
          sqsum_ax(ie) = sqsum_ax(ie) + nspx**2
          sum_ay(ie)   = sum_ay(ie) + nspy
          sqsum_ay(ie) = sqsum_ay(ie) + nspy**2
          nampl(ie)    = nampl(ie) + 1
        else
          nspx = zero
          nspy = zero
        end if
      end do
    end if
  
    ! Note: Not in firstrun
    if(dowritetracks) then
      do j=1,napx
        xj  = (xv1(j)-torbx(ie)) /c1e3
        xpj = (yv1(j)-torbxp(ie))/c1e3
  
negative.code:   real(kind=fPrec), private, save :: scale_by0 = zero
    real(kind=fPrec), private, save :: bx_dist   = zero
    real(kind=fPrec), private, save :: by_dist   = zero
  
  #ifdef HDF5
    ! Variables to save hdf5 dataset indices
  