subroutine initial_PhononSimulation()
  use particles
  use scattering
  use simulation
  integer :: nAlBin
  do nAlBin = 1, (num_BinsForAl+1)
    call TempScParaSet(nAlBin)
  enddo
  call InitializeAPTable()
  if (LOmodeOccPhonSim.eq.1) then ! Included in initialize If we choose LOmode1.. kp should be smaller 
    ConstForEneOp = RepresentEpLO*e_c*BinSize/Default_primV ! J/m^2 
  else 
    ConstForEneOp = RepresentEpLO*e_c*BinSize*3.0/Default_primV ! J/m^2    
  endif  
  qsupplyO = HotPhononFactor*qsupply*k_pO/(k_pA+k_pO)  ! A/k_GaAs of heat flows through
  qsupplyA = qsupply - qsupplyO
  return	
end

subroutine initial_PhononTdis()
  use particles
  use scattering
  use simulation
  integer :: i,nBin,Temp_Bin_Lower,Temp_Bin_Upper
  real ( kind = 8 ) :: UnitVolAPEne
  UnitVolAPEne = 0.0d0
  if (QuasiBoundConsid.eq.1)then   ! total energy  Hot Phonon factor 
    Temp_Bin_Lower = num_boundBin+1
    Temp_Bin_Upper = num_Bins-num_boundBin
  else
    Temp_Bin_Lower = 1
    Temp_Bin_Upper = num_Bins
  endif 
  if (PhononSim.eq.1) then ! left constant fixed 
    do nBin = Temp_Bin_Lower,Temp_Bin_Upper ! BinSize
      if (k_pO.gt.0) then 
        Bin_TemperOP(nBin) = TpO_Hot - BinSize*(nBin-1)*qsupplyO/k_pO
      else
        Bin_TemperOP(nBin) = TpO_Hot
      endif 
      Bin_TemperAP(nBin) = TpA_Hot - BinSize*(nBin-1)*qsupplyA/k_pA
    enddo   
  elseif (PhononSim.eq.2) then ! left constant fixed 
    do nBin = Temp_Bin_Lower,Temp_Bin_Upper !col_len = RCutoff_SiIn
      if (k_pO.gt.0) then 
        Bin_TemperOP(num_Bins+1-nBin) = TpO_Cold + BinSize*(nBin-1)*qsupplyO/k_pO
      else
        Bin_TemperOP(num_Bins+1-nBin) = TpO_Cold
      endif 
      Bin_TemperAP(num_Bins+1-nBin) = TpA_Cold + BinSize*(nBin-1)*qsupplyA/k_pA      
    enddo   
  endif 
  RightBoundTempOP = Bin_TemperOP(num_Bins) + (Bin_TemperOP(num_Bins)-Bin_TemperOP(num_Bins-1))
  RightBoundTempAP = Bin_TemperAP(num_Bins) + (Bin_TemperAP(num_Bins)-Bin_TemperAP(num_Bins-1))
!  ExternalResAP = (Bin_TemperAP(num_Bins) - T_room)/qsupplyA
!  if (k_pO.ne.0) then
!    ExternalResOP = ExternalResAP*(k_pA/k_pO)
!  else
!    ExternalResOP = 0
!  endif
    ! PhononSim = 1 Hot Temp constant/ 2 Cold Temp Constant
  do nBin = Temp_Bin_Lower,Temp_Bin_Upper
    Bin_OccupOP(nBin) = 1/(exp((RepresentEpLO*e_c)/(kB*Bin_TemperOP(nBin))) - 1)
    BinQsupO(nBin) = qsupplyO 
    BinQsupA(nBin) = qsupplyA
    BinQsupTotal(nBin) = qsupply
    call APTempToEneCal(Bin_TemperAP(nBin),UnitVolAPEne)
    if (nBin.eq.1) then ! first
      BinAPEne(nBin) = UnitVolAPEne*(0.5*BinSize)
      BinOPEne(nBin) = RepresentEpLO*e_c*Bin_OccupOP(nBin)*(0.5*BinSize)/Default_primV
    elseif (nBin.eq.num_Bins) then   ! last bin
      BinAPEne(nBin) = UnitVolAPEne*(0.5*BinSize)
      BinOPEne(nBin) = RepresentEpLO*e_c*Bin_OccupOP(nBin)*(0.5*BinSize)/Default_primV
    else ! middle of first last bin
      BinAPEne(nBin) = UnitVolAPEne*BinSize
      BinOPEne(nBin) = RepresentEpLO*e_c*Bin_OccupOP(nBin)*BinSize/Default_primV
    endif
  enddo
  if (QuasiBoundConsid.eq.1)then 
    do nBin = 1,num_boundBin !col_len = RCutoff_SiIn
      BinQsupA(nBin) = 0
      BinQsupO(nBin) = 0
      BinQsupTotal(nBin) = 0
      BinAPEne(nBin) = BinAPEne(Temp_Bin_Lower)
      BinOPEne(nBin) = BinOPEne(Temp_Bin_Lower)
      Bin_OccupOP(nBin) = Bin_OccupOP(Temp_Bin_Lower)
      Bin_TemperAP(nBin) = Bin_TemperAP(Temp_Bin_Lower)
      Bin_TemperOP(nBin) = Bin_TemperOP(Temp_Bin_Lower)
    enddo
    do nBin = Temp_Bin_Upper+1,num_Bins !col_len = RCutoff_SiIn
      BinQsupA(nBin) = 0
      BinQsupO(nBin) = 0
      BinQsupTotal(nBin) = 0
      BinAPEne(nBin) = BinAPEne(Temp_Bin_Upper)
      BinOPEne(nBin) = BinOPEne(Temp_Bin_Upper)
      Bin_OccupOP(nBin) = Bin_OccupOP(Temp_Bin_Upper)
      Bin_TemperAP(nBin) = Bin_TemperAP(Temp_Bin_Upper)
      Bin_TemperOP(nBin) = Bin_TemperOP(Temp_Bin_Upper)
    enddo  
  endif  
  do nBin = 1,num_BinsComb   ! PhononSim = 1 Hot Temp constant/ 2 Cold Temp Constant
    call TempVar_sc_table(nBin)  ! constant temperature outside them 
  enddo 
  return	
end
subroutine TempScParaSet(nAlBin)
  use particles
  use scattering
  use simulation
  integer :: i,val,valf, nFinalValleys
  real ( kind = 8 ) :: rhouu, initEe_eV, initNP_Ee_eV, sigma, acou_const, POIntact_eV, PolarConst, impurity_const, alloy_const, gammaSquare
  real ( kind = 8 ) :: PolarAb_Ee_eV, PolarAbNP_Ee_eV, POAb_rnum, POAb_denom, POAb_factor, PolarEm_Ee_eV, PolarEmNP_Ee_eV, POEm_rnum, POEm_denom, POEm_factor 
  real ( kind = 8 ) :: ivIntact_eV, coupling_const, delta_Efi_eV, ivconst,ivAb_Ee_eV, ivAbNP_Ee_eV, ivAb_factor, ivEm_Ee_eV, ivEmNP_Ee_eV, ivEm_factor
  rhouu = Al_rho_GaAs(nAlBin)*Al_u_sound(nAlBin)*Al_u_sound(nAlBin)  ! rho u^2 
  do val = 1, num_Valley ! intravalley interaction
    if (val.eq.1) then 
        POIntact_eV = polar_Ep_G
        sigma = acou_sigma_G        
    elseif (val.eq.2) then 
        POIntact_eV = polar_Ep_L
        sigma = acou_sigma_L        
    else 
        POIntact_eV = polar_Ep_X
        sigma = acou_sigma_X        
    endif
    acou_const = sqrt(2.)*kB/(pi*hbar*rhouu)*(eff_mass(nAlBin,val)/hbar)*(sqrt(eff_mass(nAlBin,val))/hbar*sqrt(e_c))*(echbar*e_c)*sigma*sigma !Tp_eV
    PolarConst = echbar*echbar*e_c*POIntact_eV*sqrt(eff_mass(nAlBin,val)/2./e_c)/4./pi*(1./nBin_eps_high(nAlBin) - 1./nBin_eps_low(nAlBin)) ! for polar optical eps_high,eps_low (4.*pi) ?
  !  impurity_const = n_d*eff_mass(nAlBin,val)*echbar*sqrt(2.*eff_mass(nAlBin,val))/pi*echbar*DebyeEps(nAlBin,val)*echbar*DebyeEps(nAlBin,val)*echbar*sqrt(e_c)
    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,val)*sqrt(eff_mass(nAlBin,val)) 
    endif
    do i = 1, num_EnergyLevel
      initEe_eV = StepE_eV*real(i,kind = 8) ! initial energy 
      initNP_Ee_eV = initEe_eV*(1.+nonParabolPara(nAlBin,val)*initEe_eV) ! for nonparabolic 
      PolarAb_Ee_eV = initEe_eV + POIntact_eV 
      PolarEm_Ee_eV = initEe_eV - POIntact_eV 
      acou_constTemp(nAlBin,val,i) = acou_const*sqrt(initNP_Ee_eV)*(1.+nonParabolPara2(nAlBin,val)*initEe_eV)
      
      gammaSquare = 8*eff_mass(nAlBin,val)*e_c*initEe_eV*Debye_length(nAlBin)*Debye_length(nAlBin)/hbar/hbar ! change 1
      impurity_const = (log(1+gammaSquare) - gammaSquare/(gammaSquare+1))*(sqrt(initEe_eV*e_c)*initEe_eV*e_c) ! change 1
      impurity_constTemp(nAlBin,val,i) = Energy_debye(nAlBin,val)/impurity_const ! impurity_const*sqrt(initNP_Ee_eV)*(1.+nonParabolPara2(nAlBin,val)*initEe_eV)/(1.+4.*initNP_Ee_eV/Energy_debye(nAlBin,val)) ! Add 10/2013
      if (impurity_constTemp(nAlBin,val,i).gt.1.5D15) then ! limit the rate
        impurity_constTemp(nAlBin,val,i) =  1.5D15
      endif
      alloy_constTemp(nAlBin,val,i)  =  alloy_const*sqrt(initNP_Ee_eV)*(1.+nonParabolPara2(nAlBin,val)*initEe_eV) 
          
      PolarAbNP_Ee_eV = PolarAb_Ee_eV*(1.+nonParabolPara(nAlBin,val)*PolarAb_Ee_eV) ! for nonparabolic 
      
      if (PolarOpticalTest.eq.1) then
        POAb_rnum = sqrt(initEe_eV) + sqrt(PolarAb_Ee_eV) !
        POAb_denom = sqrt(POIntact_eV) ! 
        POAb_factor = (1.+nonParabolPara2(nAlBin,val)*PolarAb_Ee_eV)/sqrt(initNP_Ee_eV)*log(abs(POAb_rnum/POAb_denom))
      elseif (PolarOpticalTest.eq.2) then
        POAb_rnum = sqrt(initNP_Ee_eV) + sqrt(PolarAbNP_Ee_eV) !
        POAb_denom = sqrt(abs(initNP_Ee_eV - PolarAbNP_Ee_eV)) ! 
        POAb_factor = (1.+nonParabolPara2(nAlBin,val)*PolarAb_Ee_eV)/sqrt(initNP_Ee_eV)*log(abs(POAb_rnum/POAb_denom))
      else
        POAb_rnum = sqrt(initNP_Ee_eV) + sqrt(PolarAbNP_Ee_eV) !
        POAb_denom = sqrt(initNP_Ee_eV) - sqrt(PolarAbNP_Ee_eV) ! 
        POAb_factor = (1.+nonParabolPara2(nAlBin,val)*PolarAb_Ee_eV)/sqrt(initNP_Ee_eV)*log(abs(POAb_rnum/POAb_denom))
      endif
     
      polarAb_constTemp(nAlBin,val,i) = POAb_factor*PolarConst 
      
      if(PolarEm_Ee_eV.le.0)then  ! PolarConst,polar_ab,polar_ab_rate,polar_em,polar_em_rate
        polarEm_constTemp(nAlBin,val,i) = 0
      else 
        PolarEmNP_Ee_eV = PolarEm_Ee_eV*(1.+nonParabolPara(nAlBin,val)*PolarEm_Ee_eV) ! for nonparabolic  
        if (PolarOpticalTest.eq.1) then
          POEm_rnum = sqrt(initEe_eV) + sqrt(PolarEm_Ee_eV) !
          POEm_denom = sqrt(POIntact_eV) ! 
          POEm_factor = (1.+nonParabolPara2(nAlBin,val)*PolarEm_Ee_eV)/sqrt(initNP_Ee_eV)*log(abs(POEm_rnum/POEm_denom))
        elseif (PolarOpticalTest.eq.2) then
          POEm_rnum = sqrt(initNP_Ee_eV) + sqrt(PolarEmNP_Ee_eV) !
          POEm_denom = sqrt(abs(initNP_Ee_eV - PolarEmNP_Ee_eV)) ! 
          POEm_factor = (1.+nonParabolPara2(nAlBin,val)*PolarEm_Ee_eV)/sqrt(initNP_Ee_eV)*log(abs(POEm_rnum/POEm_denom))
        else
          POEm_rnum = sqrt(initNP_Ee_eV) + sqrt(PolarEmNP_Ee_eV) !
          POEm_denom = sqrt(initNP_Ee_eV) - sqrt(PolarEmNP_Ee_eV) ! 
          POEm_factor = (1.+nonParabolPara2(nAlBin,val)*PolarEm_Ee_eV)/sqrt(initNP_Ee_eV)*log(abs(POEm_rnum/POEm_denom))
        endif
        polarEm_constTemp(nAlBin,val,i) = POEm_factor*PolarConst       
      endif
    enddo 
  enddo
  do val = 1, num_Valley ! intravalley interaction
    do valf = 1, num_Valley ! final valley 
      if (val.eq.1) then
        if (valf.eq.1) then
          ivIntact_eV = 0
          coupling_const = 0
          delta_Efi_eV = 0
          nFinalValleys = 0
        elseif (valf.eq.2) then
          ivIntact_eV = IntVal_Ep_G_L 
          coupling_const = DefPot_G_L 
          delta_Efi_eV = Al_split_L_G(nAlBin) ! energy difference eV
          nFinalValleys = eq_valleys_L ! # of valleys
        else 
          ivIntact_eV = IntVal_Ep_G_X 
          coupling_const = DefPot_G_X 
          delta_Efi_eV = Al_split_X_G(nAlBin) ! energy difference eV
          nFinalValleys = eq_valleys_X ! # of valleys         
        endif
      elseif (val.eq.2) then  ! L -> 
        if (valf.eq.1) then
          ivIntact_eV = IntVal_Ep_L_G 
          coupling_const = DefPot_L_G 
          delta_Efi_eV = -Al_split_L_G(nAlBin) ! energy difference eV
          nFinalValleys = eq_valleys_G ! # of valleys
        elseif (valf.eq.2) then
          ivIntact_eV = IntVal_Ep_L_L 
          coupling_const = DefPot_L_L 
          delta_Efi_eV = 0. ! energy difference eV
          nFinalValleys = eq_valleys_L - 1. ! # of valleys
        else 
          ivIntact_eV = IntVal_Ep_L_X 
          coupling_const = DefPot_L_X 
          delta_Efi_eV = Al_split_X_G(nAlBin) - Al_split_L_G(nAlBin) ! energy difference eV
          nFinalValleys = eq_valleys_X ! # of valleys         
        endif
      else ! X -> other 
        if (valf.eq.1) then
          ivIntact_eV = IntVal_Ep_X_G 
          coupling_const = DefPot_X_G 
          delta_Efi_eV = -Al_split_X_G(nAlBin) ! energy difference eV
          nFinalValleys = eq_valleys_G ! # of valleys
        elseif (valf.eq.2) then
          ivIntact_eV = IntVal_Ep_X_L 
          coupling_const = DefPot_X_L 
          delta_Efi_eV = Al_split_L_G(nAlBin) - Al_split_X_G(nAlBin) ! energy difference eV
          nFinalValleys = eq_valleys_L ! # of valleys
        else 
          ivIntact_eV = IntVal_Ep_X_X 
          coupling_const = DefPot_X_X 
          delta_Efi_eV = 0. ! energy difference eV
          nFinalValleys = eq_valleys_X - 1. ! # of valleys         
        endif
      endif
      if (nFinalValleys.eq.0) then
        ivAb_constTemp(nAlBin,val,valf,:) = 0
        ivEm_constTemp(nAlBin,val,valf,:) = 0
      else
        ivconst = nFinalValleys*(coupling_const**2.)*e_c*sqrt(e_c)/(sqrt(2.)*pi*Al_rho_GaAs(nAlBin)*ivIntact_eV)*(eff_mass(nAlBin,valf)/hbar)*sqrt(eff_mass(nAlBin,valf))/hbar    
        do i = 1, num_EnergyLevel
          initEe_eV = StepE_eV*real(i,kind = 8) ! initial energy 
          initNP_Ee_eV = initEe_eV*(1.+nonParabolPara(nAlBin,val)*initEe_eV) ! for nonparabolic 
          ivAb_Ee_eV = initEe_eV + ivIntact_eV - delta_Efi_eV
          ivEm_Ee_eV = initEe_eV - ivIntact_eV - delta_Efi_eV
          if(ivAb_Ee_eV.le.0)then  ! tempEef_eV (final energy) < 0 absorption 0
            ivAb_constTemp(nAlBin,val,valf,i) = 0.
          else
            ivAbNP_Ee_eV = ivAb_Ee_eV*(1.+ nonParabolPara(nAlBin,valf)*ivAb_Ee_eV) 
            ivAb_factor = sqrt(ivAbNP_Ee_eV)*(1.+nonParabolPara2(nAlBin,valf)*ivAb_Ee_eV)
            ivAb_constTemp(nAlBin,val,valf,i) = ivAb_factor*ivconst 
          endif
          if(ivEm_Ee_eV.le.0)then  ! tempEef_eV (final energy) < 0 Emsorption 0
            ivEm_constTemp(nAlBin,val,valf,i) = 0.
          else
            ivEmNP_Ee_eV = ivEm_Ee_eV*(1.+ nonParabolPara(nAlBin,valf)*ivEm_Ee_eV) 
            ivEm_factor = sqrt(ivEmNP_Ee_eV)*(1.+nonParabolPara2(nAlBin,valf)*ivEm_Ee_eV)
            ivEm_constTemp(nAlBin,val,valf,i) = ivEm_factor*ivconst 
          endif
        enddo
      endif ! if (nFinalValleys.eq.0) then
    enddo
  enddo 
  return
end

subroutine CombBinAssign()
  use particles
  use scattering
  use simulation
  integer :: TB,Temp_TB,CB,Temp_CB,AB,Temp_AB,Temp_AB_low,Temp_AB_up,nAlBin_LowT,nAlBin_UpT
  real ( kind = 8 ) :: LowAlBound_LowT,UpAlBound_LowT,LowAlBound_UpT,UpAlBound_UpT,TempBin_LowBound,TempBin_UpBound
  Temp_CB = 0
  CombBin(:,:) = 0
  do TB = 1, num_Bins   
    TempBin_LowBound = loc_bin(TB)-0.5*BinSize
    TempBin_UpBound = loc_bin(TB)+0.5*BinSize
    call DecideAlBin(TempBin_LowBound,nAlBin_LowT,LowAlBound_LowT,UpAlBound_LowT)
    call DecideAlBin(TempBin_UpBound,nAlBin_UpT,LowAlBound_UpT,UpAlBound_UpT)
    if (nAlBin_LowT.eq.nAlBin_UpT) then 
      Temp_CB = Temp_CB + 1
      CombBin(nAlBin_LowT,TB) = Temp_CB
    else
      Temp_CB = Temp_CB + 1
      CombBin(nAlBin_LowT,TB) = Temp_CB
      if (nAlBin_LowT.eq.(num_BinsForAl+1)) then
        Temp_AB_low = 1
      else
        Temp_AB_low = 1+nAlBin_LowT
      endif
      if (LowAlBound_UpT.eq.TempBin_UpBound) then 
        Temp_AB_up =  nAlBin_UpT - 1
      else
        Temp_AB_up =  nAlBin_UpT
      endif
      if (Temp_AB_low.le.Temp_AB_up) then
        do AB = Temp_AB_low, Temp_AB_up
          Temp_CB = Temp_CB + 1
          CombBin(AB,TB) = Temp_CB 
        enddo
      endif
    endif
  enddo
  num_BinsComb = Temp_CB    
    
  allocate ( AlBinFromCB(num_BinsComb) ) 
  allocate ( TBinFromCB(num_BinsComb) ) 
    
  do CB = 1, num_BinsComb
    do AB = 1, num_BinsForAl+1
      do TB = 1, num_Bins  
        if (CombBin(AB,TB).eq.CB) then
          Temp_AB = AB 
          Temp_TB = TB 
        endif
      enddo 
    enddo
    AlBinFromCB(CB) = Temp_AB
    TBinFromCB(CB) = Temp_TB      
  enddo
  return
end