anchor.diff:              y=y+0.5*yp*(s_length-Srefl)
   !            write(*,*) "s_length-Srefl", s_length-Srefl
               CALL CALC_ION_LOSS_CRY(IS,is2,PC,s_length-Srefl,DESt)
  -            CALL MOVE_AM_(IS,is2,NAM,s_length-Srefl,DESt,DLYi(IS2),DLRi(IS2),xp ,yp,PC)
  +            CALL move_am(IS,is2,NAM,s_length-Srefl,DESt,DLYi(IS2),DLRi(IS2),xp ,yp,PC)
               x = x + 0.5 * xp * (s_length - Srefl)
               y = y + 0.5 * yp * (s_length - Srefl)
               endif
  
positive.code: 
  !            write(*,*) "s_length-Sdech", s_length-Sdech
              CALL CALC_ION_LOSS_CRY(IS,is2,PC,s_length-Sdech,DESt)
              CALL MOVE_AM_(IS,is2,NAM,s_length-Sdech,DESt,DLYi(IS2),DLRi(IS2),xp,yp,PC)
             !next line new from sasha
              x = x + 0.5*(s_length-Sdech)*xp
              y = y + 0.5*(s_length-Sdech)*yp
  
negative.code:       END
  
  
  ! ================================================================================================ !
  !  Subroutine for the movement in the amorphous
  ! ================================================================================================ !
  subroutine move_am(is,is2,nam,dz,dei,dly,dlr,xp,yp,pc)
  
    use mod_ranlux
    use mod_funlux
    use floatPrecision
    use coll_crystal
    use coll_materials, only : anuc, hcut, bnref, csref
  
    implicit none
  
    integer,          intent(in)    :: is
    integer,          intent(in)    :: is2
    integer,          intent(in)    :: nam
    real(kind=fPrec), intent(in)    :: dz
    real(kind=fPrec), intent(in)    :: dei
    real(kind=fPrec), intent(in)    :: dly
    real(kind=fPrec), intent(in)    :: dlr
    real(kind=fPrec), intent(inout) :: xp
    real(kind=fPrec), intent(inout) :: yp
    real(kind=fPrec), intent(inout) :: pc
  
    integer i,length_cry,ichoix
    real(kind=fPrec) t,xran_cry(1),bn(4),cs(0:5,1:4),freep(4),zlm,xp_in,yp_in,xm2,xln15s,tz,tx,tlow,  &
      thigh,teta,pptot,ppsd,ppel,pc_in,kymcs,kxmcs,ecmsq,dya,bsd,bpp,aran
    real(kind=fPrec), external :: ruth_cry
  
    xp_in = xp
    yp_in = yp
    pc_in = pc
  
    ! New treatment of scattering routine based on standard sixtrack routine
    ! useful calculations for cross-section and event topology calculation 
    ecmsq  = two*0.93828_fPrec*pc
    xln15s = log(0.15*ecmsq)
  
    ! New models, see Claudia's thesis
    pptot = 0.041084d0-0.0023302d0*log(ecmsq)+0.00031514d0*log(ecmsq)**2
    ppel  = (11.7d0-1.59d0*log(ecmsq)+0.134d0*log(ecmsq)**2)/1000
    ppsd  = (4.3d0+0.3d0*log(ecmsq))/1000
    bpp   = 7.156d0+1.439d0*log(sqrt(ecmsq))
  
    ! Distribution for Ruth. scatt.
    tlow      = tlcut_cry
    mcurr_cry = is2
    thigh     = hcut(is2)
    call funlxp(ruth_cry,cgen_cry(1,is),tlow,thigh)
  
    ! Cross-section calculation
    ! freep: number of nucleons involved in single scattering
    freep(is) = freeco_cry * anuc(is2)**(one/three)
  
    ! Compute pp and pn el+single diff contributions to cross-section (both added : quasi-elastic or qel later)
    cs(3,is) = freep(is) * ppel
    cs(4,is) = freep(is) * ppsd
  
    ! Correct TOT-CSec for energy dependence of qel
    ! TOT CS is here without a Coulomb contribution
    cs(0,is) = csref(0,is2) + freep(is) * (pptot - pptref_cry)
    bn(is)   = bnref(is2) * cs(0,is) / csref(0,is2)
  
    ! Also correct inel-CS
    cs(1,is) = csref(1,is2) * cs(0,is) / csref(0,is2)
  
    ! Nuclear Elastic is TOT-inel-qel ( see definition in RPP)
    cs(2,is) = cs(0,is) - cs(1,is) - cs(3,is) - cs(4,is)
    cs(5,is) = csref(5,is2)
    ! Now add Coulomb
    cs(0,is) = cs(0,is) + cs(5,is)
  
    ! Calculate cumulative probability
    do i=1,4
      cprob_cry(i,is) = cprob_cry(i-1,is)+cs(i,is)/cs(0,is)
    end do
  
    ! Multiple Coulomb Scattering
  
    xp  = xp*c1e3
    yp  = yp*c1e3
    pc  = pc-dei*dz  ! Energy lost because of ionization process[GeV]
  
    dya   = (13.6/pc)*sqrt(dz/dlr) ! rms of coloumb scattering MCS (mrad)
    kxmcs = dya*ran_gauss(one)
    kymcs = dya*ran_gauss(one)
  
    xp = xp+kxmcs
    yp = yp+kymcs
  
    if(nam == 0) return ! Turn on/off nuclear interactions
  
    ! Can nuclear interaction happen?
    zlm = -collnt_cry(is)*log(rndm4())
  
    if(zlm < dz) then
      ! Choose nuclear interaction
      aran = rndm4()
      i=1
  10  if(aran > cprob_cry(i,is)) then
        i = i+1
        goto 10
      end if
      ichoix = i
  
      ! Do the interaction
      t = 0 ! default value to cover ichoix=1
      select case(ichoix)
      case(1) ! Deep inelastic, impinging p disappeared
        iProc = proc_absorbed
  
      case(2) ! p-n elastic
        iProc = proc_pne
        t     = -log(rndm4())/bn(is)
  
      case(3) ! p-p elastic
        iProc = proc_ppe
        t     = -log(rndm4())/bpp
  
      case(4) ! Single diffractive
        iProc = proc_diff
        xm2   = exp(rndm4()*xln15s)
        pc    = pc*(one - xm2/ecmsq)
        if(xm2 < two) then
          bsd = two*bpp
        else if(xm2 >= two .and. xm2 <= five) then
          bsd = (106.d0 - 17.d0*xm2)* bpp/36.d0
        else if(xm2 > five) then
          bsd = 7.d0*bpp/12.d0
        end if
        t = -log(rndm4())/bsd
  
      case(5)
        iProc      = proc_ruth
        length_cry = 1
        call funlux(cgen_cry(1,is),xran_cry,length_cry)
        t = xran_cry(1)
  
      end select
  
      ! Calculate the related kick
      if(ichoix == 4) then
        teta = sqrt(t)/pc_in ! DIFF has changed PC
      else
        teta = sqrt(t)/pc
      end if
  
      tx = teta*ran_gauss(one)*c1e3
      tz = teta*ran_gauss(one)*c1e3
  
      ! Change p angle
      xp = xp + tx
      yp = yp + tz
    end if
  
    xp = xp/c1e3
    yp = yp/c1e3
  
  end subroutine move_am
  
  ! ================================================================================================ !
  !  Subroutine for check if a nuclear interaction happen while in channeling
  