!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!    Subroutine for the calculation of acoustic phonons scattering rate
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine acoustic_rate(nAlBin,i_count,valley,sigma,out_file) ! Al Content Variation Add 08/2013 
   use simulation
    use scattering
    integer i_count,valley,i,nAlBin
    character*30 out_file !  ,  
    real ( kind = 8 ) :: sigma,rhouu,acou_const, acou_rate, tempEe_eV, tempNP_Ee_eV
    !   Calculate constant
    rhouu = Al_rho_GaAs(nAlBin)*Al_u_sound(nAlBin)*Al_u_sound(nAlBin)  ! rho u^2 ! Al Content Variation Add 08/2013 
    acou_const = sqrt(2.)*kB*Tsys/(pi*hbar*rhouu)*(eff_mass(nAlBin,valley)/hbar)*(sqrt(eff_mass(nAlBin,valley))/hbar*sqrt(e_c))*(echbar*e_c)*sigma*sigma !Tp_eV
    i_count = i_count + 1 ! sigma   !    Create scattering table  
    if (PrintoutScatRate.eq.1) then ! 0709
      open(unit = 10, file=out_file, status='unknown')
      write(10,*)'energy  ',out_file
    endif
    do i = 1, num_EnergyLevel
     tempEe_eV = StepE_eV*real(i,kind = 8)
     tempNP_Ee_eV = tempEe_eV*(1.+nonParabolPara(nAlBin,valley)*tempEe_eV) ! for nonparabolic... ! Al Content Variation Add 08/2013 
     acou_rate = acou_const*sqrt(tempNP_Ee_eV)*(1.+nonParabolPara2(nAlBin,valley)*tempEe_eV) ! Al Content Variation Add 08/2013 
     scatt_table(nAlBin,i,i_count,valley) = acou_rate  ! Al Content Variation Add 08/2013 
     if (PrintoutScatRate.eq.1) then ! 0709
       write(10,'(2X,F8.4,2X,E14.6)'),tempEe_eV,acou_rate
     endif
    enddo
    if (PrintoutScatRate.eq.1) then ! 0709
      close(10)
    endif
    flag_mech(i_count,valley) = 1 ! isotropic scattering
    E_change(i_count,valley) = 0. ! energy change, elastic 
    i_valley(i_count,valley) = valley ! intravalley 
    return
end
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     Subroutine for the calculation of the POLAR OPTICAL PHONONS scattering rate (absorption + emission)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine polar_rate(nAlBin,i_count,valley,EpIntact_eV,out_file_1,out_file_2) ! valley,G-1,L-2,X-3
 use simulation
 use scattering
 integer nAlBin,i_count,valley,i
 character*30 out_file_1, out_file_2
 real ( kind = 8 ) :: EpIntact_eV  ! parameters
 real ( kind = 8 ) :: tempEe_eV,tempEef_eV,tempNP_Ee_eV,tempNP_Eef_eV  ! local variables - energy
 real ( kind = 8 ) :: PolarConst,polar_ab,polar_ab_rate,polar_em,polar_em_rate  ! local variables - rate 
 real ( kind = 8 ) :: f_ph,rnum,denom,factor  ! local variables - etc
!   Calculate constant    
 f_ph = 1./(exp(EpIntact_eV/Tp_eV)-1.)   ! phonon population
 PolarConst = echbar*echbar*e_c*EpIntact_eV*sqrt(eff_mass(nAlBin,valley)/2./e_c)/4./pi*(1./nBin_eps_high(nAlBin) - 1./nBin_eps_low(nAlBin)) ! for polar optical eps_high,eps_low (4.*pi) ?
!     (a) Scattering rate - absorption
 i_count = i_count + 1  ! 1 - Absorption 2 - Emission
 if (PrintoutScatRate.eq.1) then ! 0709
   open(unit=10, file=out_file_1, status='unknown')
   write(10,*)'energy  ',out_file_1
 endif
 polar_ab = f_ph*PolarConst 
 do i = 1, num_EnergyLevel ! total energy level 
     tempEe_eV = StepE_eV*real(i,kind = 8) ! energy level 
     tempEef_eV = tempEe_eV + EpIntact_eV     ! Final carrier energy
     tempNP_Ee_eV = tempEe_eV*(1.+nonParabolPara(nAlBin,valley)*tempEe_eV) ! for Init momentum w/ nonparabolicity
     tempNP_Eef_eV = tempEef_eV*(1.+nonParabolPara(nAlBin,valley)*tempEef_eV)! for Final momentum w/ nonparabolicity
     if (PolarOpticalTest.eq.1) then
       rnum = sqrt(tempEe_eV) + sqrt(tempEef_eV) !
       denom = sqrt(EpIntact_eV) ! 
       factor = (1.+nonParabolPara2(nAlBin,valley)*tempEef_eV)/sqrt(tempNP_Ee_eV)*log(abs(rnum/denom))
     elseif (PolarOpticalTest.eq.2) then
       rnum = sqrt(tempNP_Ee_eV) + sqrt(tempNP_Eef_eV) !
       denom = sqrt(abs(tempNP_Ee_eV - tempNP_Eef_eV)) ! 
       factor = (1.+nonParabolPara2(nAlBin,valley)*tempEef_eV)/sqrt(tempNP_Ee_eV)*log(abs(rnum/denom))
     else
       rnum = sqrt(tempNP_Ee_eV) + sqrt(tempNP_Eef_eV) !
       denom = sqrt(tempNP_Ee_eV) - sqrt(tempNP_Eef_eV) ! 
       factor = (1.+nonParabolPara2(nAlBin,valley)*tempEef_eV)/sqrt(tempNP_Ee_eV)*log(abs(rnum/denom))
     endif
 !    A = (2*(1.+nonParabolPara2(valley)*tempEe_eV)*(1.+nonParabolPara(valley)*tempEef_eV) + nonParabolPara(valley)*(tempNP_Ee_eV+tempNP_Eef_eV))**2.
 !    B = -nonParabolPara2(valley)*sqrt(tempNP_Ee_eV*tempNP_Eef_eV)*(4.*(1.+nonParabolPara(valley)*tempEe_eV)*(1.+nonParabolPara(valley)*tempEef_eV) + nonParabolPara(valley)*(tempNP_Ee_eV+tempNP_Eef_eV))
 !    C = 4.*(1.+nonParabolPara(valley)*tempEe_eV)*(1.+nonParabolPara(valley)*tempEef_eV)*(1.+nonParabolPara2(valley)*tempEe_eV)*(1.+nonParabolPara2(valley)*tempEef_eV)
 !    A = 4.      C = 4.      B = 0.
 !    factor = (1.+nonParabolPara2(valley)*tempEef_eV)/sqrt(ge)*(A*log(abs(rnum/denom))+B)/C
     polar_ab_rate = polar_ab*factor
     scatt_table(nAlBin,i,i_count,valley) = polar_ab_rate
     if (PrintoutScatRate.eq.1) then ! 0709
       write(10,'(2X,F8.4,2X,E14.6)'),tempEe_eV,polar_ab_rate
     endif
  enddo
  if (PrintoutScatRate.eq.1) then ! 0709
    close(10)
  endif  
  flag_mech(i_count,valley) = 2 ! anisotropic polar scattering
  E_change(i_count,valley) = EpIntact_eV ! energy change
  i_valley(i_count,valley) = valley ! final valley - intravalley 
!    (b) Scattering rate - emission
  i_count = i_count + 1
  if (PrintoutScatRate.eq.1) then ! 0709
    open(unit=11, file=out_file_2, status='unknown')
    write(11,*)'energy  ',out_file_2
  endif
  polar_em = (1.+f_ph)*PolarConst
  do i = 1, num_EnergyLevel
     tempEe_eV = StepE_eV*real(i,kind = 8)
     tempEef_eV = tempEe_eV - EpIntact_eV
     if(tempEef_eV.le.0)then  ! PolarConst,polar_ab,polar_ab_rate,polar_em,polar_em_rate
       polar_em_rate = 0
     else !tempEe_eV,tempEef_eV,tempNP_Ee_eV,tempNP_Eef_eV 
       tempNP_Ee_eV = tempEe_eV*(1.+nonParabolPara(nAlBin,valley)*tempEe_eV)
       tempNP_Eef_eV = tempEef_eV*(1.+nonParabolPara(nAlBin,valley)*tempEef_eV)
       if (PolarOpticalTest.eq.1) then
         rnum = sqrt(tempEe_eV) + sqrt(tempEef_eV) !
         denom = sqrt(EpIntact_eV) ! 
         factor = (1.+nonParabolPara2(nAlBin,valley)*tempEef_eV)/sqrt(tempNP_Ee_eV)*log(abs(rnum/denom))
       elseif (PolarOpticalTest.eq.2) then
         rnum = sqrt(tempNP_Ee_eV) + sqrt(tempNP_Eef_eV) !
         denom = sqrt(abs(tempNP_Ee_eV - tempNP_Eef_eV)) ! 
         factor = (1.+nonParabolPara2(nAlBin,valley)*tempEef_eV)/sqrt(tempNP_Ee_eV)*log(abs(rnum/denom))
       else
         rnum = sqrt(tempNP_Ee_eV) + sqrt(tempNP_Eef_eV) !
         denom = sqrt(tempNP_Ee_eV) - sqrt(tempNP_Eef_eV) ! 
         factor = (1.+nonParabolPara2(nAlBin,valley)*tempEef_eV)/sqrt(tempNP_Ee_eV)*log(abs(rnum/denom))
       endif  
 !       A = (2*(1.+nonParabolPara2(valley)*tempEe_eV)*(1.+nonParabolPara(valley)*tempEef_eV)+ nonParabolPara(valley)*(tempNP_Ee_eV+tempNP_Eef_eV))**2.
 !       B = -nonParabolPara2(valley)*sqrt(tempNP_Ee_eV*tempNP_Eef_eV)* (4.*(1.+nonParabolPara(valley)*tempEe_eV)*(1.+nonParabolPara(valley)*tempEef_eV)+ nonParabolPara(valley)*(tempNP_Ee_eV+tempNP_Eef_eV))
 !       C = 4.*(1.+nonParabolPara(valley)*tempEe_eV)*(1.+nonParabolPara(valley)*tempEef_eV)*(1.+nonParabolPara2(valley)*tempEe_eV)*(1.+nonParabolPara2(valley)*tempEef_eV)
 !       A = 4.         C = 4.         B = 0.  !       factor = (1.+af2(valley)*tempEef_eV)/sqrt(ge)*(A*log(abs(rnum/denom))+B)/C 
       polar_em_rate = polar_em*factor
     endif
     scatt_table(nAlBin,i,i_count,valley) = polar_em_rate
     if (PrintoutScatRate.eq.1) then ! 0709
       write(11,'(2X,F8.4,2X,E14.6)'),tempEe_eV,polar_em_rate
     endif
  enddo
  if (PrintoutScatRate.eq.1) then ! 0709
    close(11)                
  endif
  flag_mech(i_count,valley) = 2 ! anisotropic polar scattering
  E_change(i_count,valley) = -EpIntact_eV ! energy change     
  i_valley(i_count,valley) = valley ! final valley - intravalley 
  return
end   
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Generic subroutine for the calculation of INTERVALLEY PHONONS scattering rate  (absorption + emission)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine intervalley_rate(nAlBin,i_count,valley_i,valley_f,nFinalValleys,coupling_const,delta_Efi_eV,EpIntact_eV,out_file_1,out_file_2)
    use simulation
    use scattering
    integer nAlBin,i_count,valley_i,valley_f,nFinalValleys,i ! initial and final valley
    real ( kind = 8 ) :: coupling_const,delta_Efi_eV,EpIntact_eV
    character*30 out_file_1, out_file_2
    real ( kind = 8 ) :: f_ph,final_mass,ivconst,iv_ab,iv_em,tempEe_eV,tempEef_eV,tempNP_Eef_eV,factor,iv_ab_rate,iv_em_rate
    !    Calculate constants
    f_ph = 1./(exp(EpIntact_eV/Tp_eV)-1.)
    final_mass = eff_mass(nAlBin,valley_f)
    ivconst = nFinalValleys*(coupling_const**2.)*e_c*sqrt(e_c)/(sqrt(2.)*pi*Al_rho_GaAs(nAlBin)*EpIntact_eV)*(final_mass/hbar)*sqrt(final_mass)/hbar    
    !  (a) Scattering rate - absorption
    i_count = i_count + 1
    if (PrintoutScatRate.eq.1) then ! 0709
      open(unit=10, file=out_file_1, status='unknown')
      write(10,*)'energy  ',out_file_1
    endif
    iv_ab = f_ph*ivconst
    do i = 1, num_EnergyLevel
      tempEe_eV = StepE_eV*real(i,kind = 8)
      tempEef_eV = tempEe_eV + EpIntact_eV - delta_Efi_eV
      tempNP_Eef_eV = tempEef_eV*(1.+ nonParabolPara(nAlBin,valley_f)*tempEef_eV)
      if(tempEef_eV.le.0)then  ! tempEef_eV (final energy) < 0 absorption 0
        iv_ab_rate = 0.
      else
        factor = sqrt(tempNP_Eef_eV)*(1.+nonParabolPara2(nAlBin,valley_f)*tempEef_eV)
        iv_ab_rate = iv_ab*factor
      endif
      scatt_table(nAlBin,i,i_count,valley_i) = iv_ab_rate
      if (PrintoutScatRate.eq.1) then ! 0709
        write(10,'(2X,F8.4,2X,E14.6)'),tempEe_eV,iv_ab_rate
      endif
    enddo
    if (PrintoutScatRate.eq.1) then ! 0709
      close(10)      
    endif
    flag_mech(i_count,valley_i) = 1
    E_change(i_count,valley_i) = EpIntact_eV - delta_Efi_eV
    i_valley(i_count,valley_i) = valley_f
    !    (b) Scattering rate - emission
    i_count = i_count + 1
    if (PrintoutScatRate.eq.1) then ! 0709
      open(unit=11, file=out_file_2, status='unknown')
      write(11,*)'energy  ',out_file_2
    endif
    iv_em = (1.+f_ph)*ivconst
    do i = 1, num_EnergyLevel
      tempEe_eV = StepE_eV*real(i,kind = 8)
      tempEef_eV = tempEe_eV - EpIntact_eV - delta_Efi_eV
      tempNP_Eef_eV = tempEef_eV*(1.+nonParabolPara(nAlBin,valley_f)*tempEef_eV)
      if(tempEef_eV.le.0)then !  tempEef_eV < 0 0 
        iv_em_rate = 0.
      else
        factor = sqrt(tempNP_Eef_eV)*(1.+nonParabolPara2(nAlBin,valley_f)*tempEef_eV)
        iv_em_rate = iv_em*factor
      endif
      scatt_table(nAlBin,i,i_count,valley_i) = iv_em_rate
      if (PrintoutScatRate.eq.1) then ! 0709
        write(11,'(2X,F8.4,2X,E14.6)'),tempEe_eV,iv_em_rate
      endif
    enddo
    if (PrintoutScatRate.eq.1) then ! 0709
      close(11)
    endif
    flag_mech(i_count,valley_i) = 1 ! isotropic scattering
    E_change(i_count,valley_i) = - EpIntact_eV - delta_Efi_eV ! energy change     
    i_valley(i_count,valley_i) = valley_f ! final valley
    return
end

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     Subroutine for the calculation of the Ionized Impurity / Alloy_Scattering 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine ionimpurity_rate(nAlBin,i_count,valley,out_file_1)  ! valley,G-1,L-2,X-3
    use simulation
    use scattering
    integer i_count,valley,i,nAlBin
    character*30 out_file !  ,  
    real ( kind = 8 ) ::  impurity_const, impurity_rate, tempEe_eV, tempNP_Ee_eV, gammaSquare, Enefactor   !   Calculate constant
    ! impurity_const = n_d*eff_mass(nAlBin,valley)*echbar*sqrt(2.*eff_mass(nAlBin,valley))/pi*echbar*DebyeEps(nAlBin,valley)*echbar*DebyeEps(nAlBin,valley)*echbar*sqrt(e_c)
    i_count = i_count + 1 ! ! acou_const = sqrt(2.)*kB*Tsys/(pi*hbar*rhouu)*(eff_mass(nAlBin,valley)/hbar)*(sqrt(eff_mass(nAlBin,valley))/hbar*sqrt(e_c))*(echbar*e_c)*sigma*sigma !Tp_eV sigma   !    Create scattering table  
    if (PrintoutScatRate.eq.1) then ! 0709
      open(unit = 10, file=out_file, status='unknown')
      write(10,*)'energy  ',out_file
    endif
    do i = 1, num_EnergyLevel
     tempEe_eV = StepE_eV*real(i,kind = 8)
     tempNP_Ee_eV = tempEe_eV*(1.+nonParabolPara(nAlBin,valley)*tempEe_eV) ! for nonparabolic... ! Al Content Variation Add 08/2013 
     gammaSquare = 8*eff_mass(nAlBin,valley)*e_c*tempEe_eV*Debye_length(nAlBin)*Debye_length(nAlBin)/hbar/hbar
     impurity_const = (log(1+gammaSquare) - gammaSquare/(gammaSquare+1))*(sqrt(tempEe_eV*e_c)*tempEe_eV*e_c)
     impurity_rate = Energy_debye(nAlBin,valley)/impurity_const
     if (impurity_rate.gt.1.5D15) then ! limit the rate
        impurity_rate =  1.5D15
     endif
    ! impurity_rate = impurity_const*sqrt(tempNP_Ee_eV)*(1.+nonParabolPara2(nAlBin,valley)*tempEe_eV)/(1.+4.*tempNP_Ee_eV/Energy_debye(nAlBin,valley)) ! Add 10/2013
     scatt_table(nAlBin,i,i_count,valley) = impurity_rate  ! Al Content Variation Add 08/2013 
     if (PrintoutScatRate.eq.1) then ! 0709
       write(10,'(2X,F8.4,2X,E14.6)'),tempEe_eV,impurity_rate
     endif
    enddo
    if (PrintoutScatRate.eq.1) then ! 0709
      close(10)
    endif
    flag_mech(i_count,valley) = 3 ! isotropic scattering
    E_change(i_count,valley) = 0. ! energy change, elastic 
    i_valley(i_count,valley) = valley ! intravalley 
    return
end   

subroutine alloyscat_rate(nAlBin,i_count,valley,out_file_1) ! valley,G-1,L-2,X-3
    use simulation
    use scattering
    integer i_count,valley,i,nAlBin
    character*30 out_file !    
    real ( kind = 8 ) :: alloy_const,alloy_rate,tempEe_eV, tempNP_Ee_eV !   Calculate constant
    if (Al_Dist_Content(nAlBin).eq.0.0) then
        alloy_const = 0  
    else
        alloy_const = AlloyPotDiff*e_c*AlloyPotDiff*e_c*sqrt(2.*e_c)*Al_PrimV(nAlBin)*Al_Dist_Content(nAlBin)*(1-Al_Dist_Content(nAlBin))/(pi*hbar*hbar*hbar*hbar)*eff_mass(nAlBin,valley)*sqrt(eff_mass(nAlBin,valley)) 
    endif   !    acou_const = sqrt(2.)*kB*Tsys/(pi*hbar*rhouu)*(eff_mass(nAlBin,valley)/hbar)*(sqrt(eff_mass(nAlBin,valley))/hbar*sqrt(e_c))*(echbar*e_c)*sigma*sigma !Tp_eV
    i_count = i_count + 1 ! sigma   !    Create scattering table  
    if (PrintoutScatRate.eq.1) then ! 0709
      open(unit = 10, file=out_file, status='unknown')
      write(10,*)'energy  ',out_file
    endif
    do i = 1, num_EnergyLevel
     tempEe_eV = StepE_eV*real(i,kind = 8)
     tempNP_Ee_eV = tempEe_eV*(1.+nonParabolPara(nAlBin,valley)*tempEe_eV) ! for nonparabolic... ! Al Content Variation Add 08/2013 
     alloy_rate  =  alloy_const*sqrt(tempNP_Ee_eV)*(1.+nonParabolPara2(nAlBin,valley)*tempEe_eV) 
     scatt_table(nAlBin,i,i_count,valley) = alloy_rate  ! Al Content Variation Add 08/2013 
     if (PrintoutScatRate.eq.1) then ! 0709
       write(10,'(2X,F8.4,2X,E14.6)'),tempEe_eV,alloy_rate
     endif
    enddo
    if (PrintoutScatRate.eq.1) then ! 0709
      close(10)
    endif
    flag_mech(i_count,valley) = 1 ! isotropic scattering
    E_change(i_count,valley) = 0. ! energy change, elastic 
    i_valley(i_count,valley) = valley ! intravalley 
    return
end