anchor.diff:      end if
       ichoix = i
   
  -    ! Do the interaction 
  +    ! Do the interaction
       select case(ichoix)
       case(1) ! deep inelastic, impinging p disappeared
         iProc = proc_ch_absorbed
  -    
  +
       case(2) ! p-n elastic
         iProc  = proc_ch_pne
         bn(IS) = bnref(is2) * cs(0,IS) / csref_tot_rsc
  
positive.code:       end
  
  
  !.**************************************************************************
  !     SUBROUTINE FOR THE MOVEMENTS OF THE PARTICLES IN THE CRYSTAL
  !.**************************************************************************
        SUBROUTINE CRYST(IS,is2,x,xp,y,yp,PC,Length,j)
  
        use mod_ranlux
        use mod_funlux
        use mod_common_main
        use floatPrecision
        use coll_crystal
        use coll_materials, only : zatom, exenergy, rho, anuc
  
  !     Simple tranport protons in crystal 2
  !-----------------------------------------------------------C
  !      J -   number of element                              C
  !      S - longitudinal coordinate
  !      IS -   number of substance 1-4: Si,W,C,Ge(110)       C
  !      x,xp,y,yp - coordinates at input of crystal          C
  !      PC -   momentum of particle*c [GeV]                  C
  !      W  -   weigth of particle                            C
  !-----------------------------------------------------------C
  !
  !
        IMPLICIT none
  !
  !
        real(kind=fPrec) Length !crystal geometrical parameters
                                                    ! [m],[m],[m],[m],[rad]
        real(kind=fPrec) ymax,ymin       !crystal geometrical parameters
        real(kind=fPrec) s_length             !element length along s
        integer IS,j,is2                            !index of the material
  !      integer counter
        real(kind=fPrec) DESt                  ! Daniele: changed energy loss by ionization now calculated and not tabulated
        integer NAM,ZN                        !switch on/off the nuclear
                                              !interaction (NAM) and the MCS (ZN)
        real(kind=fPrec) x,xp,y,yp,PC         !coordinates of the particle
                                              ![m],[rad],[m],[rad],[GeV]
        real(kind=fPrec) x0,y0                !coordinates of the particle [m]
        real(kind=fPrec) s                    !long coordinates of the particle [m]
        real(kind=fPrec) a_eq,b_eq,c_eq,Delta !second order equation param.
        real(kind=fPrec) Ang_rms, Ang_avr     !Volume reflection mean angle [rad]
        real(kind=fPrec) c_v1                 !fitting coefficient
        real(kind=fPrec) c_v2                 !fitting coefficient
        real(kind=fPrec) Dechan               !probability for dechanneling
        real(kind=fPrec) Lrefl, Srefl         !distance of the reflection point [m]
        real(kind=fPrec) Vcapt                !volume capture probability
        real(kind=fPrec) Chann                !channeling probability
        real(kind=fPrec) N_atom               !probability for entering
                                              !channeling near atomic planes
        real(kind=fPrec) Dxp                  !variation in angle
        real(kind=fPrec) xpcrit               !critical angle for curved crystal[rad]
        real(kind=fPrec) xpcrit0              !critical angle for str. crystal [rad]
        real(kind=fPrec) Rcrit                !critical curvature radius [m]
        real(kind=fPrec) ratio                !x=Rcurv/Rcrit
        real(kind=fPrec) Cry_length
        real(kind=fPrec) TLdech2              !tipical dechanneling length(1) [m]
        real(kind=fPrec) TLdech1              !tipical dechanneling length(2) [m]
        real(kind=fPrec) tdech, Ldech,Sdech   !angle, lenght, and S coordinate
                                              !of dechanneling point
        real(kind=fPrec) Rlength, Red_S       !reduced length/s coordinate
                                              !(in case of dechanneling)
        real(kind=fPrec) Am_length            !Amorphous length
        real(kind=fPrec) Length_xs, Length_ys !Amorphous length
        real(kind=fPrec) L_chan, tchan
        real(kind=fPrec) xp_rel               !xp-miscut angle in mrad
  
        real(kind=fPrec) const_dech,xpin,ypin
        real(kind=fPrec) alpha !Daniele: par for new chann prob
        real(kind=fPrec) Pvr !Daniele: prob for VR->AM transition
  
  !
  !
        NAM=1 !switch on/off the nuclear interaction (NAM) and the MCS (ZN)
        ZN=1
  
  !-------------Daniele: dE/dX and dechanneling length calculation--------------------
  
  
  !      write(*,*) "p entering routine", PC
        mom=PC*1.0d3 ! [GeV/c] -> [MeV/c]
        enr=(mom*mom+mp*mp)**0.5 ! [MeV]
        gammar=enr/mp
        betar=mom/enr
        bgr=betar*gammar
  
        Tmax=(2.0d0*me*bgr**2)/(1.0d0+2*gammar*me/mp+(me/mp)**2) ![MeV]
  
        plen=((rho(is2)*zatom(is2)/anuc(is2))**0.5)*28.816d-6 ![MeV]
  
        const_dech=(256.0/(9.0*(4.D0*DATAN(1.D0))**2))* &
       & (1.0/(log(2.0*me*gammar/(exenergy(is2)*c1e3))-1.0))*((aTF*dP)/(re*me))  ![m/MeV]
  
        const_dech=const_dech*1.0d3    ![m/GeV]
  
  !      write(*,*)DESt, const_dech
  
  !----------------------------------------------------------
  
  
  !      miscut=0.001000
  !
  !      write(*,*)"last miscut angle =",miscut
  !      write(*,*) 'enter crystal subroutine'
  !      write(*,*) 'particle energy Gev :', PC
  !      write(*,*) 'x_initial :', x
  !      write(*,*) 'Length [m]:', Length
  !      write(*,*) 'Random:', rndm4()
  !      write(*,*)'xp',xp,'x',x , 's', s
  !      DESt = 0
        s=0
  !      write(*,*) "s_length", Rcurv, length
        s_length=Rcurv*(sin(length/Rcurv)) !
        L_chan=length
  
        if ( miscut .lt. 0 &
       &     .and. x .gt. 0 & !should be useless
       &     .and. x .lt. -length*tan(miscut)) then
              L_chan=-x/sin(miscut)
        endif
        tchan=L_chan/Rcurv
        xp_rel=xp-miscut
  !  FIRST CASE: p don't interact with crystal
        ymin = - C_ymax/2
        ymax =  C_ymax/2
        IF (y.LT.ymin .or. y.GT.ymax .or. x.gt.C_xmax) THEN
          x = x+xp*s_length
          y = y+yp*s_length
          iProc = proc_out
          GOTO 111
  ! SECOND CASE: p hits the amorphous layer
        ELSEIF ( (x.LT.Alayer) .or.  ((y-ymin).LT.Alayer) .or. ((ymax-y).lt.Alayer)  ) THEN
          x0=x
          y0=y
          a_eq=(1+(xp)**2)
          b_eq=(2*(x)*(xp)-2*(xp)*Rcurv)
          c_eq=(x)**2-2*(x)*Rcurv
          Delta=b_eq**2-4*a_eq*c_eq
          s=((-b_eq+sqrt(Delta))/(2*a_eq))
          if (s .ge. s_length) s=s_length
          x=(xp)*s+x0
          Length_xs=sqrt((x-x0)**2+s**2)
          if ( (yp .ge.0 .and. (y+yp*s).le.ymax)) then
            Length_ys = yp*Length_xs
          elseif (yp.lt.0 .and. (y+yp*s).ge. ymin) then
            Length_ys = yp*Length_xs
          else
            s=(ymax-y)/yp
            Length_ys = sqrt((ymax-y)**2+s**2)
            x=x0+xp*s
            Length_xs=sqrt((x-x0)**2+s**2)
          endif
          Am_length   = sqrt(Length_xs**2+Length_ys**2)
          s=s/2
          x=x0+xp*s
          y=y0+yp*s
          iProc = proc_AM
          CALL CALC_ION_LOSS_CRY(IS,is2,PC,AM_Length,DESt)
          CALL move_am(IS,is2,NAM,Am_Length,DESt,DLYi(IS2),dlri(is2),xp,yp,PC)
          x=x+xp*(s_length-s)
          y=y+yp*(s_length-s)
          GOTO 111
        ELSEIF ((x.GT.(C_xmax-Alayer)) .and. x.LT.(C_xmax)  ) THEN
          iProc = proc_AM
          CALL CALC_ION_LOSS_CRY(IS,is2,PC,s_length,DESt)
          CALL move_am(IS,is2,NAM,s_length,DESt,DLYi(IS2),DLRi(IS2), xp,yp,PC)
          WRITE(*,*)'Fix here!'
          GOTO 111
        END IF
  !
  ! THIRD CASE: the p interacts with the crystal.
  !. Define typical angles/probabilities for orientation 110
  !
        xpcrit0 = (2.e-9*eUm(IS2)/PC)**0.5       ! critical angle (rad) for
                                                ! straight crystals
        Rcrit  = PC/(2.e-6*eUm(IS2))*AI(IS2)      ! critical curvature radius [m]
                                                ! if R>Rcritical=>no channeling is
                                                ! possible (ratio<1)
        ratio = Rcurv/Rcrit                     ! parameter Rcry/Rcritical
  !      write(*,*) "Critical Radius: ",Rcrit
        xpcrit = xpcrit0*(Rcurv-Rcrit)/Rcurv    ! critical angle for curved crystal
  !----------------valentina approx-----------
  !      xpcrit = xpcrit0*(1-(Rcrit/Rcurv))**0.5
  !      write(*,*)(Rcurv-Rcrit)/Rcurv,(1-(Rcrit/Rcurv))**0.5
  
                                                ! NB: if ratio<1 => xpcrit<0
        c_v1 = 1.7                              ! fitting coefficient ??!
        c_v2 = -1.5                             ! fitting coefficient ???
        if (ratio .le. 1.) then                 ! case 1:no possibile channeling
          Ang_rms = c_v1*0.42*xpcrit0*sin(1.4*ratio)  ! rms scattering
          Ang_avr = c_v2*xpcrit0*0.05*ratio           ! average angle reflection
          Vcapt = 0.0                                 ! probability of VC
          elseif (ratio .le. 3) then              ! case 2: strongly bent xstal
            Ang_rms = c_v1*0.42*xpcrit0*sin(1.571*0.3*ratio+0.85)! rms scattering
            Ang_avr = c_v2*xpcrit0*(0.1972*ratio-0.1472)  ! avg angle reflection
  !          Vcapt   = 0.01*(ratio-0.7)/(PC**2)
            Vcapt   = 0.0007*(ratio-0.7)/PC**0.2 !correction by sasha drozdin/armen
            !K=0.00070 is taken based on simulations using CATCH.f (V.Biryukov)
          else                                       ! case 3: Rcry >> Rcrit
            Ang_rms = c_v1*xpcrit0*(1./ratio)        !
            Ang_avr = c_v2*xpcrit0*(1.-1.6667/ratio) ! average angle for VR
  !          Vcapt = 0.01*(ratio-0.7)/(PC**2)        ! probability for VC
            Vcapt = 0.0007*(ratio-0.7)/PC**0.2  !correction by sasha drozdin/armen
            ! K=0.0007 is taken based on simulations using CATCH.f (V.Biryukov)
        endif
  !c----------------valentina approx-----------
  !      Ang_avr=-(xpcrit+xpcrit0)
  !      Ang_rms=(xpcrit0-xpcrit)/2
  !c-----------end valentina approx--------------
  
  !      write(*,*) "Rcrit" , Rcrit,"Rcurv",Rcurv,
  !     c "Ratio: ",ratio,"average VR angle:", ang_avr*1e6,"+-",
  !     c ang_rms*1e6, "ang crit:", xpcrit0*1e6,xpcrit*1e6
  !
        if(C_orient .eq. 2) then
          Ang_avr = Ang_avr * 0.93                     ! for (111)
          Ang_rms = Ang_rms * 1.05
          xpcrit  = xpcrit * 0.98
        endif
  !
  !. case 3-1: channeling
  !      IF (abs(xp_rel) .lt. xpcrit) THEN              ! if R' < R'c (ok CH) (1)
  !        Chann  = (xpcrit**2-xp_rel**2)**0.5/xpcrit0  ! probability of CH/VC  OCCHIO DANIELE LINEA COMMENTATA PER PROVARE LA SUCCESSIVA PROB.
  !        N_atom = 0.1                                ! probability of entering OCCHIO DANIELE, prob cambiata in accordo con nuova chann
  !--------------DAN CHAN prob------
  !         alpha = xp_rel/xpcrit
  !         Chann = ((0.64-(1/ratio)*(1/ratio)*0.64)**0.5)*     !DANIELE saturation at 80%
  !     &           (1-alpha*alpha)**0.5
  
  !         Chann = ((0.8-(1/ratio)*(1/ratio)*0.8)**0.5)*     !DANIELE saturation at 90%
  !     &           (1-alpha*alpha)**0.5
  
  !         N_atom=0.14                                      !DANIELE for sat. at 80%
  !         N_atom=0.1                                       !DANIELE for sat. at 90%
  
  !      IF (abs(xp_rel) .lt. xpcrit) THEN
  !         Chann = ((0.9-(1/ratio)*(1/ratio)*0.9)**0.5)*     !DANIELE saturation at 95%
  !     &           (1-alpha*alpha)**0.5
  !        N_atom = 0.1
  
        IF (abs(xp_rel) .lt. xpcrit) THEN
           alpha = xp_rel/xpcrit
           Chann = ((0.9-alpha*alpha*0.9)**0.5)*(1-(1/ratio))**0.5         !DANIELE saturation at 95%
           N_atom = 0.1
  
  !      IF (abs(xp_rel) .lt. xpcrit) THEN
  !         alpha = xp_rel/xpcrit
  !         Chann = ((0.8-alpha*alpha*0.8)**0.5)*(1-(1/ratio))**0.5         !DANIELE saturation at 90%
  !         N_atom = 0.1
  
  !--------------end DAN CHAN prob------
  
  
  
  
                                                       ! close to atomic planes
          IF (rndm4() .le. Chann) then      ! if they can channel: 2 options
                                            ! option 1:channeling
  !          TLdech1= 0.00054*PC*(1.-1./ratio)**2 ! calculate dechanneling length
  !          TLdech1= 0.0005*PC*(1.-1./ratio)**2 !calculate tipical dech. length(m)
            TLdech1= const_dech*PC*(1.-1./ratio)**2 !Daniele: updated calculate tipical dech. length(m)
            IF (rndm4() .le. N_atom) then
  !            TLdech1= 0.000004*PC*(1.-1./ratio)**2! calculate tipical dechanneling
                                                    !length near atomic planes(m)
             !next line new from sasha
  !            TLdech1= 2.0e-6*PC*(1.-1./ratio)**2  ! dechanneling length (m)
              TLdech1= (const_dech/200.d0)*PC*(1.-1./ratio)**2  ! Daniele: updated dechanneling length (m)
  
                                 !for short crystal for high amplitude particles
            ENDIF
  
  !          TLdech1=TLdech1/100 !!!!CHECK
  
            Dechan = -log(rndm4())                 ! probability of dechanneling
            Ldech  = TLdech1*Dechan                ! actual dechan. length
                       ! careful: the dechanneling lentgh is along the trajectory
                       ! of the particle -not along the longitudinal coordinate...
            if(Ldech .LT. L_chan) THEN
              iProc = proc_DC
              Dxp= Ldech/Rcurv             ! change angle from channeling [mrad]
              Sdech=Ldech*cos(miscut+0.5*Dxp)
  
              x  = x+ Ldech*(sin(0.5*Dxp+miscut))   ! trajectory at channeling exit
              xp = xp + Dxp + 2.0*(rndm4()-0.5)*xpcrit
              y= y + yp * Sdech
  !            write(*,*) "Ldech", Ldech
  !            write(*,*) "DESt", DESt
              CALL CALC_ION_LOSS_CRY(IS,is2,PC,Ldech,DESt)
              PC = PC - 0.5*DESt*Ldech          ! energy loss to ionization while in CH [GeV]
  
              x = x + 0.5*(s_length-Sdech)*xp
              y = y + 0.5*(s_length-Sdech)*yp
  
  !            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
            else
              iProc = proc_CH
              xpin=XP
              ypin=YP
  
  !            write(*,*) "angles before entering CH subroutine (x,y):", xp, yp
              CALL MOVE_CH(IS,is2,NAM,L_chan,X,XP,YP,PC,Rcurv,Rcrit)  !daniele:check if a nuclear interaction happen while in CH
  !            write(*,*) "angles after exiting CH subroutine (x,y):", xp, yp
  !            stop
  
              if(iProc /= proc_CH) then             !daniele: if an nuclear interaction happened, move until the middle with initial xp,yp
              x = x + 0.5 * L_chan * xpin           !then propagate until the "crystal exit" with the new xp,yp
              y = y + 0.5 * L_chan * ypin           !accordingly with the rest of the code in "thin lens approx"
              x = x + 0.5 * L_chan * XP
              y = y + 0.5 * L_chan * YP
  !            write(*,*) "Length", Length
              CALL CALC_ION_LOSS_CRY(IS,is2,PC,Length,DESt)
              PC = PC - DESt*Length       ! energy loss to ionization [GeV]
              else
              Dxp= L_chan/Rcurv + 0.5*RAN_GAUSS(1.0d0)*xpcrit ! change angle[rad]
              xp = Dxp
              !next line new from sasha
              x  = x+ L_chan*(sin(0.5*Dxp+miscut)) ! trajectory at channeling exit
  !            xp = xp + Dxp + 2.0*(rndm4()-0.5)*xpcrit
              y = y + s_length * yp
              CALL CALC_ION_LOSS_CRY(IS,is2,PC,Length,DESt)
              PC = PC - 0.5*DESt*Length       ! energy loss to ionization [GeV]
              endif
            endif
          ELSE                                   !option 2: VR
                                                 ! good for channeling
                                                 ! but don't channel         (1-2)
            iProc = proc_VR
  !          Dxp=0.5*(xp_rel/xpcrit+1)*Ang_avr
  !          xp=xp+Dxp+Ang_rms*RAN_GAUSS(1.)
              !next line new from sasha
            xp=xp+0.45*(xp/xpcrit+1)*Ang_avr
            x = x + 0.5*s_length * xp
            y = y + 0.5*s_length * yp
            CALL CALC_ION_LOSS_CRY(IS,is2,PC,s_length,DESt)
            CALL move_am(IS,is2,NAM,s_length,DESt,DLYi(IS2),DLRi(IS2),xp ,yp,PC)
            x = x + 0.5*s_length * xp
            y = y + 0.5*s_length * yp
          ENDIF                                    !
  ! case 3-2: no good for channeling. check if the  can VR
        ELSE
          Lrefl =  (xp_rel)*Rcurv                  ! distance of refl. point [m]
  !        Srefl = sin(xp) * Lrefl
          Srefl = sin(xp_rel/2+miscut) * Lrefl
          if(Lrefl .gt. 0. .and. Lrefl .lt. Length) then
                  ! VR point inside
                  !2 options: volume capture and volume reflection
            IF (rndm4() .gt. Vcapt .or. ZN .eq. 0.) THEN   !opt. 1: VR
              iProc = proc_VR
              x = x + xp * Srefl
              y = y + yp * Srefl
              Dxp= Ang_avr
              xp = xp + Dxp + Ang_rms*RAN_GAUSS(1.0d0)
              x = x + 0.5* xp * (s_length - Srefl)
              y = y + 0.5* yp * (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)
              x = x + 0.5 * xp * (s_length - Srefl)
              y = y + 0.5 * yp * (s_length - Srefl)
            ELSE                                      !opt 2: VC
              x = x + xp * Srefl
              y = y + yp * Srefl
  !            TLdech2= 0.00011*PC**0.25*(1.-1./ratio)**2 ! dechanneling length(m)
  !            Dechan = log(1.-rndm4())
  !            Ldech  = -TLdech2*Dechan
             !next 2 lines new from sasha - different dechanneling
             !probability
   !           TLdech2= 0.01*PC*(1.-1./ratio)**2   ! typical dechanneling length(m)
   !           Ldech  = 0.005*TLdech2*(sqrt(0.01-log(rndm4())) -0.1)**2 ! DC length
              TLdech2= (const_dech/10.0d0)*PC*(1.-1./ratio)**2   ! Daniele: updated typical dechanneling length(m)
              Ldech  = TLdech2*(sqrt(0.01-log(rndm4())) -0.1)**2 ! daniele: updated DC length
              tdech=Ldech/Rcurv
              Sdech=Ldech*cos(xp+0.5*tdech)
              IF(Ldech .LT. (Length-Lrefl)) then
                iProc = proc_DC
                Dxp= Ldech/Rcurv + 0.5*ran_gauss(1.0d0)*xpcrit
                x  = x+ Ldech*(sin(0.5*Dxp+xp))   ! trajectory at channeling exit
                y = y + Sdech * yp
                xp =  Dxp
                Red_S = s_length-Srefl -Sdech
                x = x + 0.5 * xp * Red_S
                y = y + 0.5 * yp * Red_S
  !              write(*,*) "Srefl", Srefl
                CALL CALC_ION_LOSS_CRY(IS,is2,PC,Srefl,DESt)
                PC=PC - DESt * Srefl !Daniele: "added" energy loss before capture
  !              write(*,*) "Sdech", Sdech
                CALL CALC_ION_LOSS_CRY(IS,is2,PC,Sdech,DESt)
                PC=PC - 0.5 * DESt * Sdech !Daniele: "added" energy loss while captured
                CALL CALC_ION_LOSS_CRY(IS,is2,PC,Red_S,DESt)
                CALL move_am(IS,is2,NAM,Red_S,DESt,DLYi(IS2),DLRi(IS2),xp,yp,PC)
                x = x + 0.5 * xp * Red_S
                y = y + 0.5 * yp * Red_S
              else
                iProc = proc_VC
                Rlength = Length-Lrefl
                tchan = Rlength / Rcurv
                Red_S=Rlength*cos(xp+0.5*tchan)
  !              write(*,*) "Lrefl", Lrefl
                CALL CALC_ION_LOSS_CRY(IS,is2,PC,Lrefl,DESt)
                PC=PC - DESt*Lrefl  !Daniele: "added" energy loss before capture
                xpin=XP
                ypin=YP
                CALL MOVE_CH(IS,is2,NAM,Rlength,X,XP,YP,PC,Rcurv,Rcrit)  !daniele:check if a nuclear interaction happen while in CH
  
                if(iProc /= proc_VC) then             !daniele: if an nuclear interaction happened, move until the middle with initial xp,yp
                x = x + 0.5 * Rlength * xpin           !then propagate until the "crystal exit" with the new xp,yp
                y = y + 0.5 * Rlength * ypin           !accordingly with the rest of the code in "thin lens approx"
                x = x + 0.5 * Rlength * XP
                y = y + 0.5 * Rlength * YP
  !              write(*,*) "Rlength", Rlength
                CALL CALC_ION_LOSS_CRY(IS,is2,PC,Rlength,DESt)
                PC=PC - DESt*Rlength
                else
                Dxp = (Length-Lrefl)/Rcurv
  !              write(*,*) "Dxp", Dxp
                x  = x+ sin(0.5*Dxp+xp)*Rlength     ! trajectory at channeling exit
                y = y + red_S * yp
                xp =  Length/Rcurv + 0.5*ran_gauss(1.0d0)*xpcrit ! [mrad]
  !              write(*,*) "xp at channeling exit", xp
  !              write(*,*) "Rlength", Rlength
                CALL CALC_ION_LOSS_CRY(IS,is2,PC,Rlength,DESt)
                PC=PC - 0.5*DESt*Rlength  !Daniele: "added" energy loss once captured
                endif
              endif
            ENDIF
  !.  case 3-3: move in amorphous substance (big input angles)---------------  MODIFIED FOR TRANSITION VRAM DANIELE
          else
             if(xp_rel .GT. L_chan/Rcurv+2.0*xpcrit .or. xp_rel .lt. -xpcrit) then
               iProc = proc_AM
               x = x + 0.5 * s_length * xp
               y = y + 0.5 * s_length * yp
              if(ZN .gt. 0) then
               CALL CALC_ION_LOSS_CRY(IS,is2,PC,s_length,DESt)
               CALL move_am(IS,is2,NAM,s_length,DESt,DLYi(IS2),DLRi(IS2), xp,yp,PC)
              endif
              x = x + 0.5 * s_length * xp
              y = y + 0.5 * s_length * yp
            else
  !            Pvr=0.5*erf((xp_rel-(L_chan/Rcurv)-xpcrit)/(2.0*4.0*xpcrit*xpcrit)**0.5)+0.5
  !            Pvr=0.2+(0.6/(2.0*xpcrit))*(xp_rel-(L_chan/Rcurv))
              Pvr=((xp_rel-(L_chan/Rcurv))/(2.0*xpcrit))
              if(rndm4() .gt. Pvr) then
              iProc = proc_TRVR
              x = x + xp * Srefl
              y = y + yp * Srefl
  !            Dxp=(2.0*Ang_avr-Ang_rms)/(2.0*L_chan/Rcurv+2.0*xpcrit)*
  !     +         xp_rel+L_chan/Rcurv*(2.0*Ang_avr-Ang_rms)/
  !     +         (2.0*L_chan/Rcurv+2.0*xpcrit)-Ang_avr
              Dxp=-3.0*Ang_rms*xp_rel/(2.0*xpcrit)+Ang_avr+(3.0*Ang_rms*(L_chan/Rcurv)/(2.0*xpcrit))
  !            write(*,*) xp_rel, Dxp, Ang_avr, Ang_rms
  !            xp=xp+Dxp+Ang_rms*RAN_GAUSS(1.)
              xp=xp+Dxp
              x=x+0.5*xp*(s_length-Srefl)
              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)
              x = x + 0.5 * xp * (s_length - Srefl)
              y = y + 0.5 * yp * (s_length - Srefl)
              else
              iProc = proc_TRAM
              x = x + xp * Srefl
              y = y + yp * Srefl
              Dxp=-1.0*(13.6/PC)*SQRT(s_length/DLRi(IS2))*0.001*xp_rel/    &
       &          (2.0*xpcrit)+(13.6/PC)*SQRT(s_length/DLRi(IS2))*0.001*   &
       &          (1.0+(L_chan/Rcurv)/(2.0*xpcrit))
              xp=xp+Dxp
              x=x+0.5*xp*(s_length-Srefl)
              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)
              x = x + 0.5 * xp * (s_length - Srefl)
              y = y + 0.5 * yp * (s_length - Srefl)
              endif
            endif
          endif
         ENDIF
  
  !       write(*,*) "p at the end of routine", PC
  
  !      if (counter .eq. 0) then
  111   write(833,*)'crystal parameters:\n Length:',Length, '\n Rcurv:'   &
       & , Rcurv ,'\n Critical Radius:', Rcrit, 'ratio',ratio             &
       &, '\n Critical angle for straight:',                              &
       & xpcrit0,'\n critical angle for curved crystal:', xpcrit,' \n Leng&
       &th:', Length, '\n xmax:', C_xmax, ' ymax:', ymax, '  C_orient: '  &
       &, C_orient                                                        &
       &, '\n Avg angle reflection:', Ang_avr, '\n full channeling angle: &
       &',(Length/Rcurv)
        END
  
  ! ================================================================================================ !
  !  Subroutine for the calculazion of the energy loss by ionisation
  ! ================================================================================================ !
  subroutine calc_ion_loss_cry(is,is2,pc,dz,enlo)
  
    use mod_ranlux
    use mod_funlux
  
negative.code:   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)
  
  