!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!    Subroutine for the calculation of acoustic phonons scattering rate
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine acoustic_rate(i_count,valley,sigma,out_file)
    use simulation
    use scattering
    integer i_count,valley,i
    character*30 out_file ! rho_GaAs, u_sound
    real ( kind = 8 ) :: sigma,rhouu,acou_const, acou_rate, tempEe_eV, tempNP_Ee_eV
    !   Calculate constant
    rhouu = rho_GaAs*u_sound*u_sound  ! rho u^2
    acou_const = sqrt(2.)*kB*Tsys/(pi*hbar*rhouu)*(eff_mass(valley)/hbar)*(sqrt(eff_mass(valley))/hbar*sqrt(e_c))*(echbar*e_c)*sigma*sigma !Tp_eV
    !    Create scattering table  
    i_count = i_count + 1 ! sigma
    open(unit = 10, file=out_file, status='unknown')
    write(10,*)'energy  ',out_file
    do i = 1, num_EnergyLevel
     tempEe_eV = StepE_eV*real(i,kind = 8)
     tempNP_Ee_eV = tempEe_eV*(1.+nonParabolPara(valley)*tempEe_eV) ! for nonparabolic... 
     acou_rate = acou_const*sqrt(tempNP_Ee_eV)*(1.+nonParabolPara2(valley)*tempEe_eV)
     scatt_table(i,i_count,valley) = acou_rate
     write(10,'(2X,F8.4,2X,E14.6)'),tempEe_eV,acou_rate
    enddo
    close(10)
    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(i_count,valley,EpIntact_eV,out_file_1,out_file_2) ! valley,G-1,L-2,X-3
 use simulation
 use scattering
 integer 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(valley)/2./e_c)/4./pi*(1./eps_high - 1./eps_low) ! for polar optical eps_high,eps_low (4.*pi) ?
!     (a) Scattering rate - absorption
 i_count = i_count + 1  ! 1 - Absorption 2 - Emission
 open(unit=10, file=out_file_1, status='unknown')
 write(10,*)'energy  ',out_file_1
 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(valley)*tempEe_eV) ! for Init momentum w/ nonparabolicity
     tempNP_Eef_eV = tempEef_eV*(1.+nonParabolPara(valley)*tempEef_eV)! for Final momentum w/ nonparabolicity
     rnum = sqrt(tempNP_Ee_eV) + sqrt(tempNP_Eef_eV) !
     denom = sqrt(tempNP_Ee_eV) - sqrt(tempNP_Eef_eV) ! 
 !    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
     factor = (1.+nonParabolPara2(valley)*tempEef_eV)/sqrt(tempNP_Ee_eV)*log(abs(rnum/denom))
     polar_ab_rate = polar_ab*factor
     scatt_table(i,i_count,valley) = polar_ab_rate
     write(10,'(2X,F8.4,2X,E14.6)'),tempEe_eV,polar_ab_rate
  enddo
  close(10)             
  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
  open(unit=11, file=out_file_2, status='unknown')
  write(11,*)'energy  ',out_file_2
  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(valley)*tempEe_eV)
        tempNP_Eef_eV = tempEef_eV*(1.+nonParabolPara(valley)*tempEef_eV)
        rnum = sqrt(tempNP_Ee_eV) + sqrt(tempNP_Eef_eV)
        denom = sqrt(tempNP_Ee_eV) - sqrt(tempNP_Eef_eV)
 !       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
        factor = (1.+nonParabolPara2(valley)*tempEef_eV)/sqrt(tempNP_Ee_eV)*log(abs(rnum/denom))
        polar_em_rate = polar_em*factor
     endif
     scatt_table(i,i_count,valley) = polar_em_rate
     write(11,'(2X,F8.4,2X,E14.6)'),tempEe_eV,polar_em_rate
  enddo
  close(11)                
  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(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 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(valley_f)
    ivconst = nFinalValleys*(coupling_const**2.)*e_c*sqrt(e_c)/(sqrt(2.)*pi*rho_GaAs*EpIntact_eV)*(final_mass/hbar)*sqrt(final_mass)/hbar    
    !  (a) Scattering rate - absorption
    i_count = i_count + 1
    open(unit=10, file=out_file_1, status='unknown')
    write(10,*)'energy  ',out_file_1
    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(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(valley_f)*tempEef_eV)
        iv_ab_rate = iv_ab*factor
     endif
     scatt_table(i,i_count,valley_i) = iv_ab_rate
     write(10,'(2X,F8.4,2X,E14.6)'),tempEe_eV,iv_ab_rate
    enddo
    close(10)      
    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
    open(unit=11, file=out_file_2, status='unknown')
    write(11,*)'energy  ',out_file_2
    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(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(valley_f)*tempEef_eV)
        iv_em_rate = iv_em*factor
     endif
     scatt_table(i,i_count,valley_i) = iv_em_rate
     write(11,'(2X,F8.4,2X,E14.6)'),tempEe_eV,iv_em_rate
    enddo
    close(11)
    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