anchor.diff:  
   !new models, see Claudia's thesis
   
  -      pptot=0.041084d0-0.0023302d0*log(ecmsq)+0.00031514d0*log(ecmsq)**2     
  +      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))  
  +      ppsd=(4.3d0+0.3d0*log(ecmsq))/1000
  +      bpp=7.156d0+1.439d0*log(sqrt(ecmsq))
   
   !------ distribution for Ruth. scatt.---------
   
  
positive.code: 
  ! For calculation of dE/dX due to ionization
  
        double precision rho(4),z(4),Ime(4),anuc_cry2(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)
  
        real emr_cry(4) !nuclear radius for Ruth scatt.
  
        common/ion/rho,z,Ime,k,re,me,mp,anuc_cry2,emr_cry
  
        data (rho(i),i=1,4)/2.33d0,0.0d0,0.0d0,0.0d0/  !material density [g/cm^3] 
        data (z(i),i=1,4)/14.0d0,0.0d0,0.0d0,0.0d0/    !atomic number
        data (Ime(i),i=1,4)/173.0d-6,0.0d0,0.0d0,0.0d0/!mean exitation energy [MeV]
        data (anuc_cry2(i),i=1,4)/28.08d0,0.0d0,0.0d0,0.0d0/     !implemented only Si
  
negative.code: !----------------------------------------------------------
  
  
  !      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
        s=0
        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
          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)
  