anchor.diff:          end do
           u(ii2) = s
         end do
  -C
  +!
         c  =   362436./16777216.
         cd =  7654321./16777216.
         cm = 16777213./16777216.
  -C
  +!
         return
         END
  -C.**************************************************************************
  -C
  -C.**************************************************************************
  +!.**************************************************************************
  +!
  +!.**************************************************************************
         subroutine RNDMIN(uin,cin,cdin,cmin,iin,jin)
  -      REAL*4 uin(97),cin,cdin,cmin
  -      REAL*4 u,c,cd,cm
  +      REAL uin(97),cin,cdin,cmin
  +      REAL u,c,cd,cm
         COMMON/RANDOM/ u(97),c,cd,cm,i,j
  -C
  +!
         do kkk = 1,97
           u(kkk) = uin(kkk)
         end do
  -C
  +!
         c  = cin
         cd = cdin
         cm = cmin
         i  = iin
         j  = jin
  -C
  +!
         return
         END
  -C.**************************************************************************
  -C
  -C.**************************************************************************
  +!.**************************************************************************
  +!
  +!.**************************************************************************
         subroutine RNDMOU(UOUT,COUT,CDOUT,CMOUT,IOUT,JOUT)
  -      REAL*4 uout(97),cout,cdout,cmout
  -      REAL*4 u,c,cd,cm
  +      REAL uout(97),cout,cdout,cmout
  +      REAL u,c,cd,cm
         COMMON/RANDOM/ u(97),c,cd,cm,i,j
  -C
  +!
         do kkk = 1,97
           uout(kkk) = u(kkk)
         end do
  -C
  +!
         COUT  = C
         CDOUT = CD
         CMOUT = CM
         IOUT  = I
         JOUT  = J
  -C
  +!
         return
         END
  -C.**************************************************************************
  -C
  -C.**************************************************************************
  +!.**************************************************************************
  +!
  +!.**************************************************************************
         subroutine RNDMTE(IO)
  -C
  -C.    *******************************************************************
  -C.    *  SUBROUTINE RNDMTE(IO)                                          
  -C.    *******************************************************************
  -C
  -      REAL*4 uu(97)
  -      REAL*4 u(6),x(6),d(6)
  -      DATA u / 6533892.0 , 14220222.0 ,  7275067.0 ,
  -     +         6172232.0 ,  8354498.0 , 10633180.0 /
  -C
  +!
  +!.    *******************************************************************
  +!.    *  SUBROUTINE RNDMTE(IO)
  +!.    *******************************************************************
  +!
  +      REAL uu(97)
  +      REAL u(6),x(6),d(6)
  +      DATA u / 6533892.0 , 14220222.0 ,  7275067.0 ,  6172232.0 ,  8354498.0 , 10633180.0 /
  +!
         call RNDMOU(UU,CC,CCD,CCM,II,JJ)
         call RNDMST(12,34,56,78)
  -C
  +!
         do ii1 = 1,20000
           xx = rndm4()
         end do
  -C
  +!
         sd = 0.0
         do ii2 = 1,6
           x(II2)  = 4096.*(4096.*rndm4())
           d(II2)  = x(II2)-u(II2)
           sd = sd+d(II2)
         end do
  -C
  +!
         call RNDMIN(uu,cc,ccd,ccm,ii,jj)
         if (IO.EQ.1 .OR. SD .NE. 0.) write(6,10) (u(I),x(I),d(I),i=1,6)
  -C
  - 10   format('  === TEST OF THE RANDOM-GENERATOR ===',/,
  -     +       '    EXPECTED VALUE    CALCULATED VALUE     DIFFERENCE',/,
  -     +       6(F17.1,F20.1,F15.3,/),
  -     +       '  === END OF TEST ;',
  -     +       '  GENERATOR HAS THE SAME STATUS AS BEFORE CALLING RNDMTE')
  +!
  + 10   format('  === TEST OF THE RANDOM-GENERATOR ===',/, &
  +     &       '    EXPECTED VALUE    CALCULATED VALUE     DIFFERENCE',/, &
  +     &       6(F17.1,F20.1,F15.3,/), &
  +     &       '  === END OF TEST ;', &
  +     &       '  GENERATOR HAS THE SAME STATUS AS BEFORE CALLING RNDMTE')
         return
         END
   
  
positive.code:               x = x + 0.5 * Rlength * XP
                y = y + 0.5 * Rlength * YP
                CALL CALC_ION_LOSS_CRY(IS,PC,Rlength,DESt)
                PC=PC - DESt*Rlength 
                else
                Dxp = (Length-Lrefl)/Rcurv 
                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)*xpcrit ! [mrad]
                CALL CALC_ION_LOSS_CRY(IS,PC,Rlength,DESt)
                PC=PC - 0.5*DESt*Rlength  !Daniele: "added" energy loss once captured    
                endif
              endif                       
            ENDIF                        
  C.  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                          
               PROC='AM'
               x = x + 0.5 * s_length * xp
               y = y + 0.5 * s_length * yp
              if(ZN .gt. 0) then                        
  !           CALL MOVE_AM_(IS,NAM,s_length,DES(IS),DLYi(IS),DLRi(IS), 
  !     +          xp,yp,PC)
               CALL CALC_ION_LOSS_CRY(IS,PC,s_length,DESt)
               CALL MOVE_AM_(IS,NAM,s_length,DESt,DLYi(IS),DLRi(IS), 
       +            xp,yp,PC)
              endif
              x = x + 0.5 * s_length * xp
              y = y + 0.5 * s_length * yp            
            else
  c            Pvr=0.5*erf((xp_rel-(L_chan/Rcurv)-xpcrit)/
  c     +          (2.0*4.0*xpcrit*xpcrit)**0.5)+0.5
  c            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
              PROC='TRVR'
              x = x + xp * Srefl
              y = y + yp * Srefl
  c            Dxp=(2.0*Ang_avr-Ang_rms)/(2.0*L_chan/Rcurv+2.0*xpcrit)*
  c     +         xp_rel+L_chan/Rcurv*(2.0*Ang_avr-Ang_rms)/
  c     +         (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
  c            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)
              CALL CALC_ION_LOSS_CRY(IS,PC,s_length-Srefl,DESt)
              CALL MOVE_AM_(IS,NAM,s_length-Srefl,DESt,DLYi(IS),
       +          DLRi(IS),xp ,yp,PC)
              x = x + 0.5 * xp * (s_length - Srefl)
              y = y + 0.5 * yp * (s_length - Srefl)
              else
              PROC='TRAM'
              x = x + xp * Srefl
              y = y + yp * Srefl
  c            Dxp=-1.0*(13.6/PC)*SQRT(s_length/DLRi(IS))*xp_rel/
  c     +      (2.0*xpcrit)+(13.6/PC)*SQRT(s_length/DLRi(IS))*
  c     +      (1.0+(L_chan/Rcurv)/(2.0*xpcrit))
              Dxp=-1.0*(13.6/PC)*SQRT(s_length/DLRi(IS))*0.001*xp_rel/
       +          (2.0*xpcrit)+(13.6/PC)*SQRT(s_length/DLRi(IS))*0.001*
       +          (1.0+(L_chan/Rcurv)/(2.0*xpcrit))
  c            write(*,*) xp_rel, Dxp
  c            xp=xp+Dxp+(13.6/PC)*SQRT(s_length/DLRi(IS))*RAN_GAUSS(1.)
              xp=xp+Dxp
              x=x+0.5*xp*(s_length-Srefl)
              y=y+0.5*yp*(s_length-Srefl)
              CALL CALC_ION_LOSS_CRY(IS,PC,s_length-Srefl,DESt)
              CALL MOVE_AM_(IS,NAM,s_length-Srefl,DESt,DLYi(IS),
       +          DLRi(IS),xp ,yp,PC)
              x = x + 0.5 * xp * (s_length - Srefl)
              y = y + 0.5 * yp * (s_length - Srefl)
              endif
            endif
          endif                         
         ENDIF                                 
  
  c      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) 
  c      counter=1
  c      endif
  c      
  c      WRITE(*,*)'Crystal process: ',PROC
  c      write(*,*)'xp_final :', xp
  c      WRITE(*,*)'Crystal process: ',PROC,'Chann Angle',Ch_angle/1000,   +
  c     1'Critical angle: ', xpcrit/1000
  c     2 DLRI(IS),DLYI(IS),AI(IS),DES(IS),eUm(IS),IS,ZN,NAM,C_orient     !,W
        END
  
  C.**************************************************************************
  C     subroutine for the calculazion of the energy loss by ionization
  C.**************************************************************************
        SUBROUTINE CALC_ION_LOSS_CRY(IS,PC,DZ,EnLo)
  
        IMPLICIT none
        integer IS
        double precision PC,DZ,EnLo
        double precision rho(4),z(4),Ime(4) !Daniele: material parameters for dE/dX calculation  
        double precision k,re,me,mp !Daniele: parameters for dE/dX calculation (const,electron radius,el. mass, prot.mass)
        double precision enr,mom,betar,gammar,bgr !Daniele: energy,momentum,beta relativistic, gamma relativistic
        double precision Tmax,plen !Daniele: maximum energy tranfer in single collision, plasma energy (see pdg)
  
negative.code: 
  !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.---------
  
  