! SUBROUTINE THAT CREATES THE SCATTERING TABLE/flag_mech = 1 ==> isotropic scattering process, 2 ==> polar optical phonons
subroutine TempVar_sc_table(nBin) ! CB bin subroutine intervalley_rate(i_count,valley_i,valley_f,nFinalValleys,coupling_const,delta_Efi_eV,EpIntact_eV,out_file_1,out_file_2)
  use particles
  use scattering
  use simulation
  integer nBin,Temp_nTBin,i_count,valley_i,valley_f,err ! initial and final valley
  real ( kind = 8 ) :: Tempature_A,Tempature_O,Polarf_p,ivGLf_p,ivGXf_p,ivLGf_p,ivLLf_p,ivLXf_p,ivXGf_p,ivXLf_p,ivXXf_p
! CREATE TABLE FOR THE GAMMA VALLEY, Scattering mechanism: 1: acoustic phonons, 2-3: polar optical phonons  4-5: G-L iv  6-7: G-X iv
  valley_i = 1
  i_count = 0  
  Temp_nTBin = TBinFromCB(nBin)
  Tempature_A = Bin_TemperAP(Temp_nTBin)
  Tempature_O = Bin_TemperOP(Temp_nTBin)
  Polarf_p =  1/(exp((polar_Ep_G*e_c)/(kB*Tempature_O)) - 1)
  if (LOmodeOccPhonSim.eq.1) then ! only LOmode1.. kp should be smaller 
    ivGLf_p =  1/(exp((IntVal_Ep_G_L*e_c)/(kB*Tempature_A)) - 1)
    ivGXf_p =  1/(exp((IntVal_Ep_G_X*e_c)/(kB*Tempature_A)) - 1)
    ivLGf_p =  1/(exp((IntVal_Ep_L_G*e_c)/(kB*Tempature_A)) - 1)
    ivLLf_p =  1/(exp((IntVal_Ep_L_L*e_c)/(kB*Tempature_A)) - 1)
    ivLXf_p =  1/(exp((IntVal_Ep_L_X*e_c)/(kB*Tempature_A)) - 1)
    ivXGf_p =  1/(exp((IntVal_Ep_X_G*e_c)/(kB*Tempature_A)) - 1)
    ivXLf_p =  1/(exp((IntVal_Ep_X_L*e_c)/(kB*Tempature_A)) - 1)
    ivXXf_p =  1/(exp((IntVal_Ep_X_X*e_c)/(kB*Tempature_A)) - 1)  
  else 
    ivGLf_p =  1/(exp((IntVal_Ep_G_L*e_c)/(kB*Tempature_O)) - 1)
    ivGXf_p =  1/(exp((IntVal_Ep_G_X*e_c)/(kB*Tempature_O)) - 1)
    ivLGf_p =  1/(exp((IntVal_Ep_L_G*e_c)/(kB*Tempature_O)) - 1)
    ivLLf_p =  1/(exp((IntVal_Ep_L_L*e_c)/(kB*Tempature_O)) - 1)
    ivLXf_p =  1/(exp((IntVal_Ep_L_X*e_c)/(kB*Tempature_O)) - 1)
    ivXGf_p =  1/(exp((IntVal_Ep_X_G*e_c)/(kB*Tempature_O)) - 1)
    ivXLf_p =  1/(exp((IntVal_Ep_X_L*e_c)/(kB*Tempature_O)) - 1)
    ivXXf_p =  1/(exp((IntVal_Ep_X_X*e_c)/(kB*Tempature_O)) - 1)   
  endif  
  if (Polarf_p.ge.0) then
  else
    err = 0
  endif     
  if(acoustic_G.eq.1)then !    Acoustic phonons scattering rate
    call TempVar_acou_rate(i_count,valley_i,nBin,Tempature_A)
  endif 
  call TempVar_polar_rate(i_count,valley_i,nBin,Polarf_p)   
  if(intervalley_G_L.eq.1)then !    Intervalley scattering: G to L 
     valley_f = 2   ! final_valley number to L (2)
     call TempVar_iv_rate(i_count,valley_i,valley_f,nBin,ivGLf_p) 
  endif
  if(intervalley_G_X.eq.1)then !    Intervalley scattering: G to X valley
     valley_f = 3   ! final_valley number to L (2)
     call TempVar_iv_rate(i_count,valley_i,valley_f,nBin,ivGXf_p) 
  endif   
  if(IonImpurity_Scat_G.eq.1)then       ! Add 10/2013
     call TempVar_ionimpurity_rate(nBin,i_count,valley_i)
  endif
  if(Alloy_Scat_G.eq.1)then 
     call TempVar_alloyscat_rate(nBin,i_count,valley_i)
  endif  
  call TempVar_renorm_table(nBin,valley_i,i_count) ! renormalize
  valley_i = 2  !     CREATE TABLE FOR THE L VALLEY Scattering mechanism:   1. acoustic phonons 2-3 polar optical phonons, 4-5: L-G 6-7: L-L 8-9: L-X intervalley
  i_count = 0  
  if(acoustic_L.eq.1)then !    Acoustic phonons scattering rate
    call TempVar_acou_rate(i_count,valley_i,nBin,Tempature_A)
  endif
  if(polar_L.eq.1)then  ! polar_rate(i_count,valley,EpIntact_eV,out_file_1,out_file_2) !
     call TempVar_polar_rate(i_count,valley_i,nBin,Polarf_p) 
  endif
  if(intervalley_L_G.eq.1)then !    Intervalley scattering: L to G 
     valley_f = 1   ! final_valley number to L (2)
     call TempVar_iv_rate(i_count,valley_i,valley_f,nBin,ivLGf_p) 
  endif
  if(intervalley_L_L.eq.1)then !    Intervalley scattering: L to L 
     valley_f = 2   ! final_valley number to L (2)
     call TempVar_iv_rate(i_count,valley_i,valley_f,nBin,ivLLf_p) 
  endif
  if(intervalley_L_X.eq.1)then !    Intervalley scattering: L to X valley
     valley_f = 3   ! final_valley number to L (2)
     call TempVar_iv_rate(i_count,valley_i,valley_f,nBin,ivLXf_p) 
  endif    
  if(IonImpurity_Scat_L.eq.1)then       
     call TempVar_ionimpurity_rate(nBin,i_count,valley_i)
  endif
  if(Alloy_Scat_L.eq.1)then  
     call TempVar_alloyscat_rate(nBin,i_count,valley_i)
  endif 
  call TempVar_renorm_table(nBin,valley_i,i_count) ! renormalize
  valley_i = 3 !     CREATE TABLE FOR THE X VALLEY Scattering mechanism:  ! 1- acoustic phonons, 2-3: polar optical phonons, 4-5: X-G, 6-7: X-L, 8-9: X-X 
  i_count = 0
  if(acoustic_X.eq.1)then !    Acoustic phonons scattering rate
    call TempVar_acou_rate(i_count,valley_i,nBin,Tempature_A)
  endif
  if(polar_X.eq.1)then  ! polar_rate(i_count,valley,EpIntact_eV,out_file_1,out_file_2) !
     call TempVar_polar_rate(i_count,valley_i,nBin,Polarf_p) 
  endif
  if(intervalley_X_G.eq.1)then !    Intervalley scattering: X to L 
     valley_f = 1   ! final_valley number to L (2)
     call TempVar_iv_rate(i_count,valley_i,valley_f,nBin,ivXGf_p) 
  endif
  if(intervalley_X_L.eq.1)then !    Intervalley scattering: X to L 
     valley_f = 2   ! final_valley number to L (2)
     call TempVar_iv_rate(i_count,valley_i,valley_f,nBin,ivXLf_p) 
  endif
  if(intervalley_X_X.eq.1)then !    Intervalley scattering: X to X valley
     valley_f = 3   ! final_valley number to L (2)
     call TempVar_iv_rate(i_count,valley_i,valley_f,nBin,ivXXf_p) 
  endif    
  if(IonImpurity_Scat_X.eq.1)then       ! Add 10/2013
     call TempVar_ionimpurity_rate(nBin,i_count,valley_i)
  endif
  if(Alloy_Scat_X.eq.1)then  
     call TempVar_alloyscat_rate(nBin,i_count,valley_i)
  endif 
  call TempVar_renorm_table(nBin,valley_i,i_count) ! renormalize
  return
end

subroutine TempVar_ionimpurity_rate(nBin,i_count,valley_i) !    TempVar_ionimpurity_rate TempVar_alloyscat_rate 10/2013 
  use particles
  use scattering
  use simulation
  integer nBin, Temp_nAlBin, i, i_count, valley_i ! initial and final valley
  i_count = i_count + 1 ! sigma
  Temp_nAlBin = AlBinFromCB(nBin) 
  do i = 1, num_EnergyLevel
    TempVar_scat_table(nBin,i,i_count,valley_i) = impurity_constTemp(Temp_nAlBin,valley_i,i)   
  enddo
  return
end
subroutine TempVar_alloyscat_rate(nBin,i_count,valley_i) !    Subroutine for the calculation of acoustic phonons scattering rate ! rate calculation
  use particles
  use scattering
  use simulation
  integer nBin, Temp_nAlBin, i, i_count, valley_i ! initial and final valley
  i_count = i_count + 1 ! sigma
  Temp_nAlBin = AlBinFromCB(nBin) 
  do i = 1, num_EnergyLevel
    TempVar_scat_table(nBin,i,i_count,valley_i) = alloy_constTemp(Temp_nAlBin,valley_i,i)   
  enddo
  return
end
subroutine TempVar_acou_rate(i_count,valley_i,nBin,Tempature_A) !    Subroutine for the calculation of acoustic phonons scattering rate ! rate calculation
  use particles
  use scattering
  use simulation
  integer nBin, Temp_nAlBin, i, i_count, valley_i ! initial and final valley
  real ( kind = 8 ) :: Tempature_A
  i_count = i_count + 1 ! sigma
  Temp_nAlBin = AlBinFromCB(nBin) 
  do i = 1, num_EnergyLevel
    TempVar_scat_table(nBin,i,i_count,valley_i) = acou_constTemp(Temp_nAlBin,valley_i,i)*Tempature_A   
  enddo
  return
end

subroutine TempVar_polar_rate(i_count,valley_i,nBin,Polarf_p)  ! valley,G-1,L-2,X-3 !     Subroutine for the calculation of the POLAR OPTICAL PHONONS scattering rate (absorption + emission)
  use particles
  use scattering
  use simulation
  integer nBin,Temp_nAlBin, i, i_count, valley_i ! initial and final valley
  real ( kind = 8 ) :: Polarf_p
  i_count = i_count + 1  ! 1 - Absorption 2 - Emission
  Temp_nAlBin = AlBinFromCB(nBin) 
  do i = 1, num_EnergyLevel ! total energy level 
     TempVar_scat_table(nBin,i,i_count,valley_i) = polarAb_constTemp(Temp_nAlBin,valley_i,i)*Polarf_p !    if (TempVar_scat_table(nBin,i,i_count,valley_i).ge.0) then  !    else !      err = 0 !    endif 
  enddo
  i_count = i_count + 1
  do i = 1, num_EnergyLevel
    TempVar_scat_table(nBin,i,i_count,valley_i) = polarEm_constTemp(Temp_nAlBin,valley_i,i)*(1.+Polarf_p) !    if (TempVar_scat_table(nBin,i,i_count,valley_i).ge.0) then !    else !      err = 0 !    endif 
  enddo  
  return
end   
! Generic subroutine for the calculation of INTERVALLEY PHONONS scattering rate  (absorption + emission)
subroutine TempVar_iv_rate(i_count,valley_i,valley_f,nBin,ivf_p) 
  use simulation
  use scattering
  integer i_count,Temp_nAlBin,valley_i,valley_f,nFinalValleys,i ! initial and final valley
  real ( kind = 8 ) :: ivf_p
  i_count = i_count + 1 !  (a) Scattering rate - absorption
  Temp_nAlBin = AlBinFromCB(nBin) 
  do i = 1, num_EnergyLevel
    TempVar_scat_table(nBin,i,i_count,valley_i) = ivAb_constTemp(Temp_nAlBin,valley_i,valley_f,i)*ivf_p
  enddo    
  i_count = i_count + 1 !    (b) Scattering rate - emission
  do i = 1, num_EnergyLevel
    TempVar_scat_table(nBin,i,i_count,valley_i) = ivEm_constTemp(Temp_nAlBin,valley_i,valley_f,i)*(1.+ivf_p)
  enddo
  return
end

subroutine TempVar_renorm_table(nBin,valley_i,i_count) ! call renormalize_table(valley_i,i_count) 
    use scattering
    integer nBin, Temp_nAlBin,valley_i, i_count, i_max, i, k !, err
    real ( kind = 8 ) :: gammamax ! local 
    max_scatt_mech(valley_i) = i_count
    if(max_scatt_mech(valley_i).ge.1)then
        do k = 1, num_EnergyLevel ! scatt_table(i,i_count,valley_i)
            TempVar_renscat_table(nBin,k,1,valley_i) = TempVar_scat_table(nBin,k,1,valley_i) 
        enddo
        if (max_scatt_mech(valley_i).gt.1) then
          do i = 2, i_count
            do k = 1, num_EnergyLevel ! scatt_table(i,i_count,valley_i)
 !             if (TempVar_scat_table(nBin,k,i,valley_i).ge.0) then
              TempVar_renscat_table(nBin,k,i,valley_i) = TempVar_scat_table(nBin,k,i,valley_i) 
              TempVar_renscat_table(nBin,k,i,valley_i) = TempVar_renscat_table(nBin,k,i,valley_i) + TempVar_renscat_table(nBin,k,i-1,valley_i)
            enddo !
          enddo
        endif
        i_max = max_scatt_mech(valley_i)
        gammamax = 0.
        do i = 1,num_EnergyLevel
          if(TempVar_renscat_table(nBin,i,i_max,valley_i).gt.gammamax) gammamax = TempVar_renscat_table(nBin,i,i_max,valley_i)   ! max_tau
        enddo
        do i = 1, i_count
          do k = 1, num_EnergyLevel
            TempVar_renscat_table(nBin,k,i,valley_i) = TempVar_renscat_table(nBin,k,i,valley_i)/gammamax ! relative tau
          enddo
        enddo
        TempVar_tau_max(nBin,valley_i) = 1./gammamax  ! max !        print*,'valley index = ',valley_i,'  tau_max =',tau_max !        print*,'  '
    endif
    return
end

subroutine TempVar_scatter_carrier(nBin,Temp_wvk,Temp_ene,Temp_valley,Temp_ScatMech)   ! 0709
    use simulation
    use scattering
    use particles
    real ( kind = 8 ), dimension(DIM) :: Temp_wvk ! carrier position and wave vector k 
    integer :: Temp_valley, Temp_valley_f, Temp_ScatMech,nBin,nAlBin,loc,i_fix, i_top, i, select_mech
    real ( kind = 8 ) :: Temp_ene,rr, bound_lower, bound_upper, ee, err1    
    Temp_ScatMech = 0 ! Phonon 0709 
    nAlBin = AlBinFromCB(nBin) 
    ee = Temp_ene ! Calculate index to the scattering table
    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.TempVar_renscat_table(nBin,loc,i_top,Temp_valley))then ! 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 221  ! self-scattering
    endif
    if(rr.lt.TempVar_renscat_table(nBin,loc,1,Temp_valley))then ! the first mechanism... 
      i_fix = 1
      ScatMechFreq(i_fix,Temp_valley) = ScatMechFreq(i_fix,Temp_valley)+1 ! counting   
      goto 112
    endif
    if(i_top.gt.1)then
        do i=1,i_top-1
           bound_lower = TempVar_renscat_table(nBin,loc,i,Temp_valley)
           bound_upper = TempVar_renscat_table(nBin,loc,i+1,Temp_valley)
           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 112
           endif
        enddo
    endif
112   continue
    Temp_ScatMech = i_fix ! 0709 Perform scattering (change energy and randomize momentum)   
    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) ! i_fix... mechanism,
    elseif(select_mech.eq.2)then
      call polar_optical_angle(nAlBin,i_fix,Temp_wvk,Temp_ene,Temp_valley) 
    elseif(select_mech.eq.3)then   
      call Coulomb_angle_BH(nAlBin,i_fix,Temp_wvk,Temp_ene,Temp_valley)
    endif
    Temp_valley = Temp_valley_f !   if ((Temp_ene.le.MaxE_eV).and.(Temp_ene.gt.0)) then      else        err1 = 1  ! for debugging    endif
221   continue ! self scattering 
return
end