!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Select Scattering Mechanism & Perform the Scattering Part Modifying Particle Atributes (kx, ky, kz, iv, energy)    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine scatter_carrier(nAlBin,Temp_wvk,Temp_ene,Temp_valley,Temp_ScatMech)   ! 0709 ! Add 08/2013 for Al Content Var
    use simulation
    use scattering
    use particles
    real ( kind = 8 ), dimension(DIM) :: Temp_wvk ! carrier position and wave vector k 
    integer :: nAlBin,Temp_valley, Temp_valley_f, Temp_ScatMech ! Add 08/2013 for Al Content Var
    real ( kind = 8 ) :: Temp_ene !, Temp_tau
    integer loc, i_fix, i_top, i, select_mech
    real ( kind = 8 ):: rr, bound_lower, bound_upper, ee, err1    
    Temp_ScatMech = 0 ! Phonon 0709 
! Calculate index to the scattering table
    ee = Temp_ene
    loc = int(ee/StepE_eV)        
    if(loc.eq.0)loc=1 ! if zero energy .. it has minimum energy 
    if(loc.gt.num_EnergyLevel) loc = num_EnergyLevel ! maximum 
!     Select scattering mechanism
    i_top = max_scatt_mech(Temp_valley) ! i_count ... total # of scattering algorithms..
    call random_number(harvest = RandomNum)
    rr = RandomNum
    if(rr.ge.ren_scatt_table(nAlBin,loc,i_top,Temp_valley))then !  ! Add 08/2013 for Al Content Var  if larger than total scattering rate = sum over all scattering algorithms
     ScatMechFreq(i_top+1,Temp_valley)=ScatMechFreq(i_top+1,Temp_valley)+1 ! counting for self scattering  
     goto 222  ! self-scattering
    endif
    if(rr.lt.ren_scatt_table(nAlBin,loc,1,Temp_valley))then ! the first mechanism... 
     i_fix = 1
     ScatMechFreq(i_fix,Temp_valley) = ScatMechFreq(i_fix,Temp_valley)+1 ! counting   
     goto 111
    endif
    if(i_top.gt.1)then
        do i=1,i_top-1
           bound_lower = ren_scatt_table(nAlBin,loc,i,Temp_valley) ! Add 08/2013 for Al Content Var
           bound_upper = ren_scatt_table(nAlBin,loc,i+1,Temp_valley)! Add 08/2013 for Al Content Var
           if ((rr.ge.bound_lower).and.(rr.lt.bound_upper))then
              i_fix = i + 1  ! index
              ScatMechFreq(i_fix,Temp_valley)=ScatMechFreq(i_fix,Temp_valley)+1 ! count 
              goto 111
           endif
        enddo
    endif
111   continue
!     Perform scattering (change energy and randomize momentum)    
    Temp_ScatMech = i_fix ! 0709
    select_mech = flag_mech(i_fix,Temp_valley)
    Temp_valley_f = i_valley(i_fix,Temp_valley) !  valley destination.. 
    if (Temp_ene.ge.1) then
      err1 = 1.0
    endif
    if(select_mech.eq.1)then  
      call isotropic(nAlBin,i_fix,Temp_wvk,Temp_ene,Temp_valley,Temp_valley_f) ! Add 08/2013 for Al Content Var i_fix... mechanism,
      if ((Temp_ene.ge.0).and.(Temp_ene.le.1)) then
        err1 = 0
      else
        err1 = 1
      endif
    elseif(select_mech.eq.2)then
      call polar_optical_angle(nAlBin,i_fix,Temp_wvk,Temp_ene,Temp_valley) ! Add 08/2013 for Al Content Var
      if ((Temp_ene.ge.0).and.(Temp_ene.le.1)) then
        err1 = 0
      else
        err1 = 1
      endif
    elseif(select_mech.eq.3)then   
      call Coulomb_angle_BH(nAlBin,i_fix,Temp_wvk,Temp_ene,Temp_valley,Temp_valley_f)
      if ((Temp_ene.ge.0).and.(Temp_ene.le.1)) then
        err1 = 0
      else
        err1 = 1
      endif 
    endif
    Temp_valley = Temp_valley_f ! if ((Temp_ene.le.MaxE_eV).and.(Temp_ene.gt.0)) then      else        err1 = 1  ! for debugging    endif
222   continue ! self scattering 
return
end      

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!       SCATTERING SUBROUTINES (CHANGE ENERGY AND WAVEVECTORS OF PARTICLES)In the definition of the scattering rates, a variable called 'flag_mech' has been defined.  
!       The values assigned to this variable correspond to: 1: Isotropic scattering (acoustic, intervalley)    2: Polar_optical ==> anisotropic scattering (small angle)  3: Coulomb scattering ==> small-angle scattering
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
!  ISOTROPIC SCATTERING PROCESS - uniform probability density for scattering in all directions, Acoustic, intervalley, Alloy
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine isotropic(nAlBin,i_fix,Temp_wvk,Temp_ene,Temp_valley,Temp_valley_f) ! Add 08/2013 for Al Content Var
    use particles
    use scattering
    use simulation    
    real ( kind = 8 ), dimension(DIM) :: Temp_wvk ! carrier position and wave vector k 
    integer :: nAlBin,Temp_valley,Temp_valley_f ! Add 08/2013 for Al Content Var
    real ( kind = 8 ) :: Temp_ene !, Temp_tau
    integer i_fix, j
    real ( kind = 8 ) :: fi, ct, st, rknew    
    real ( kind = 8 ) :: Test_ene, err1 !, Temp_tau   ! Test
    real ( kind = 8 )::  sq_wvk_sum, abs_wvk, factor   
!   Update carrier energy
    if ((Temp_ene + E_change(i_fix,Temp_valley)).le.0) goto 333 
    Temp_ene = Temp_ene + E_change(i_fix,Temp_valley)
    if ((Temp_ene.ge.0).and.(Temp_ene.le.1)) then
      err1 = 0
    else
      err1 = 1
    endif 
!   Update carrier wavevector
    rknew = semh(nAlBin,Temp_valley_f)*sqrt(Temp_ene*(1.+nonParabolPara(nAlBin,Temp_valley_f)*Temp_ene)) ! Add 08/2013 for Al Content Var
    call random_number(harvest = RandomNum)
    fi = 2.*pi*RandomNum
    call random_number(harvest = RandomNum)
    ct = 1. - 2.*RandomNum
    st = sqrt(1.-ct*ct)
   ! isotropic random... 
    Temp_wvk(1) = rknew*st*cos(fi)
    Temp_wvk(2) = rknew*st*sin(fi)
    Temp_wvk(3) = rknew*ct       
    !Test
    sq_wvk_sum = 0.
    do j = 1,DIM !col_len = RCutoff_SiIn   
      sq_wvk_sum = sq_wvk_sum + Temp_wvk(j)*Temp_wvk(j) 
    enddo
    abs_wvk = sqrt(sq_wvk_sum)
    factor = hhem(nAlBin,Temp_valley_f)*sq_wvk_sum ! Add 08/2013 for Al Content Var
    Test_ene = 2*factor/(1.+sqrt(1.+nonParabolPara4(nAlBin,Temp_valley_f)*factor)) ! Add 08/2013 for Al Content Var
    if (Test_ene.ne.Temp_ene) then
      err1 = 1.0
    endif
    if (Test_ene.ge.1) then
      err1 = 1.0
    endif
333 return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!      POLAR OPTICAL PHONONS SCATTERING ANGLE 0- Randomize the polar angle according to the notes
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine polar_optical_angle(nAlBin,i_fix,Temp_wvk,Temp_ene,Temp_valley) ! Add 08/2013 for Al Content Var
    use particles
    use scattering
    use simulation    
    integer nAlBin,i_fix  ! Add 08/2013 for Al Content Var
    real ( kind = 8 ), dimension(DIM) :: Temp_wvk ! carrier position and wave vector k 
    integer :: Temp_valley
    real ( kind = 8 ) :: Temp_ene !, Temp_tau
    real ( kind = 8 ) :: k,kxy,kxp,kyp,kzp,kp
    real ( kind = 8 ) :: enew, cth0, sth0, cfi0, sfi0, ge, gnew, zeta, cth, sth, cfi, sfi
    real ( kind = 8 ) :: Test_ene, err1 !, Temp_tau    ! Test
    real ( kind = 8 )::  sq_wvk_sum, abs_wvk, factor
    integer j 
!   Update carrier energy
    enew = Temp_ene + E_change(i_fix,Temp_valley)
    if (enew.le.0) goto 332 
!   Calculate the rotation angles
    kxy = sqrt(Temp_wvk(1)*Temp_wvk(1)+Temp_wvk(2)*Temp_wvk(2))
    k = sqrt(Temp_wvk(1)*Temp_wvk(1)+Temp_wvk(2)*Temp_wvk(2)+Temp_wvk(3)*Temp_wvk(3)) ! total k 
    cth0 = Temp_wvk(3)/k  ! cos theta  = kz/k
    sth0 = kxy/k ! sin theta  = kxy/k
    cfi0 = Temp_wvk(1)/kxy ! cos phi.. = kx/kxy
    sfi0 = Temp_wvk(2)/kxy ! sin phi.. = ky/kxy
!   Randomize momentum in the rotated coordinate system
    kp = semh(nAlBin,Temp_valley)*sqrt(enew*(1.+nonParabolPara(nAlBin,Temp_valley)*enew))
    ge = Temp_ene*(1.+nonParabolPara(nAlBin,Temp_valley)*Temp_ene)
    gnew = enew*(1.+nonParabolPara(nAlBin,Temp_valley)*enew)
    zeta = 2.*sqrt(ge*gnew)/(ge+gnew-2.*sqrt(ge*gnew))
    call random_number(harvest = RandomNum)
    cth = ((zeta+1.)-(2.*zeta+1)**RandomNum)/zeta
    sth = sqrt(1.-cth*cth)
    call random_number(harvest = RandomNum)
    fi = 2.*pi*RandomNum
    cfi = cos(fi)
    sfi = sin(fi)
    kxp = kp*sth*cfi
    kyp = kp*sth*sfi
    kzp = kp*cth
    Temp_wvk(1) = kxp*cfi0*cth0-kyp*sfi0+kzp*cfi0*sth0 !   Return back to the original coordinate system
    Temp_wvk(2) = kxp*sfi0*cth0+kyp*cfi0+kzp*sfi0*sth0
    Temp_wvk(3) = -kxp*sth0+kzp*cth0
    Temp_ene = enew   
    sq_wvk_sum = 0.      !Test
    do j = 1,DIM !col_len = RCutoff_SiIn   
       sq_wvk_sum = sq_wvk_sum + Temp_wvk(j)*Temp_wvk(j) 
    enddo
    abs_wvk = sqrt(sq_wvk_sum)
    factor = hhem(nAlBin,Temp_valley)*sq_wvk_sum
    Test_ene = 2*factor/(1.+sqrt(1.+nonParabolPara4(nAlBin,Temp_valley)*factor))
    if (Test_ene.ne.Temp_ene) then
       err1 = 1.0
    endif
    if (Test_ene.ge.1) then
       err1 = 1.0
    endif
332 return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!      COULOMB SCATTERING ANGLE - Randomize the polar angle according to the expressions derived in the notes
!      The assumption is that the scattering process is elastic
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 subroutine  Coulomb_angle_BH(nAlBin,i_fix,Temp_wvk,Temp_ene,Temp_valley)
    use particles
    use scattering
    use simulation   
    integer nAlBin,i_fix  ! Add 08/2013 for Al Content Var
    real ( kind = 8 ), dimension(DIM) :: Temp_wvk ! carrier position and wave vector k 
    integer :: Temp_valley
    real ( kind = 8 ) :: Temp_ene !, Temp_tau
    real ( kind = 8 ) :: k,kxy,kxp,kyp,kzp,kp
    real ( kind = 8 ) :: enew, cth0, sth0, cfi0, sfi0, ge, gnew, zeta, cth, sth, cfi, sfi
    real ( kind = 8 ) :: Test_ene, err1 !, Temp_tau    ! Test
    real ( kind = 8 )::  sq_wvk_sum, abs_wvk, factor
    integer j 
 !!     Update carrier energy
    enew = Temp_ene + E_change(i_fix,Temp_valley)
!     Calculate the rotation angles
    kxy = sqrt(Temp_wvk(1)*Temp_wvk(1)+Temp_wvk(2)*Temp_wvk(2))
    k = sqrt(Temp_wvk(1)*Temp_wvk(1)+Temp_wvk(2)*Temp_wvk(2)+Temp_wvk(3)*Temp_wvk(3)) ! total k 
    cth0 = Temp_wvk(3)/k  ! cos theta  = kz/k
    sth0 = kxy/k ! sin theta  = kxy/k
    cfi0 = Temp_wvk(1)/kxy ! cos phi.. = kx/kxy
    sfi0 = Temp_wvk(2)/kxy ! sin phi.. = ky/kxy
!!     Randomize momentum in the rotated coordinate system
    kp = semh(nAlBin,Temp_valley)*sqrt(enew*(1.+nonParabolPara(nAlBin,Temp_valley)*enew))
    ge = Temp_ene*(1.+nonParabolPara(nAlBin,Temp_valley)*Temp_ene) 
    call random_number(harvest = RandomNum)
    cth = 1. - 2*RandomNum/(1.+4.*(1-RandomNum)*ge/Energy_debye(nAlBin,Temp_valley))
    sth = sqrt(1.-cth*cth)
    call random_number(harvest = RandomNum)
    fi = 2.*pi*RandomNum
    cfi = cos(fi)
    sfi = sin(fi)
    kxp = kp*sth*cfi
    kyp = kp*sth*sfi
    kzp = kp*cth
    Temp_wvk(1) = kxp*cfi0*cth0-kyp*sfi0+kzp*cfi0*sth0 !   Return back to the original coordinate system
    Temp_wvk(2) = kxp*sfi0*cth0+kyp*cfi0+kzp*sfi0*sth0
    Temp_wvk(3) = -kxp*sth0+kzp*cth0
    return
end  
!      ge = e*(1.+af(iv)*e)
!      rr = ran(iso)
!      cth = 1. - 2*rr/(1.+4.*(1-rr)*ge/Energy_debye)
!      sth = sqrt(1.-cth*cth)

!      fi = two_pi*ran(iso)
!      cfi = cos(fi)
!      sfi = sin(fi)
!      kxp = k*sth*cfi
!      kyp = k*sth*sfi
!      kzp = k*cth
!!     Return back to the original coordinate system
!      kx = kxp*cfi0*cth0-kyp*sfi0+kzp*cfi0*sth0
!      ky = kxp*sfi0*cth0+kyp*cfi0+kzp*sfi0*sth0
!     kz = -kxp*sth0+kzp*cth0
!    e = enew