!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!   INITIALIZE CARRIER ENERGY AND WAVEVECTOR ACCORDING TO THE MAXWELL-BOLTZMANN STATISTICS
!   Assumption:  All carriers are initially in valley 1 for GaAs
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine initial_distribution()
    use particles
    use scattering
    use simulation    
    real ( kind = 8 ):: k,kx,ky,kz,e, sumEeV, aveEeV, rr, tc, binFreq, fai, ct, st, err1 !! local 
    
!    real ( kind = 8 ), dimension(DIM) :: Temp_wvk, sq_wvk ! carrier position and wave vector k 
!    integer :: Temp_valley_i !, Temp_valley_f
!    real ( kind = 8 ) :: Temp_ene !, Temp_tau
!    real ( kind = 8 )::  sq_wvk_sum, abs_wvk, gk
    
    integer valley, i ,j, ebin
    sumEeV = 0.
    call initial_energy() 
    do i = 1, num_Carriers
        call random_number(harvest = RandomNum)
    !    e = -(1.5*Tsys_eV)*log(RandomNum) ! carrier energy  
       !      Initial valley index
        do while(RandomNum.eq.0)
            call random_number(harvest = RandomNum)
        enddo
        ebin = 1
        do j = 2, num_EnergyLevel
            if ((RandomNum.gt.RelCumFreqDist(j-1)).and.(RandomNum.le.RelCumFreqDist(j))) then
                ebin = j                
            endif
        enddo
        if (ebin.eq.1) then 
            binFreq =  RelCumFreqDist(ebin)
            e = StepE_eV*(RandomNum/RelCumFreqDist(ebin)) !(log(RelCumFreqDist(ebin))/log(RandomNum))
        else
            binFreq =  RelCumFreqDist(ebin) - RelCumFreqDist(ebin-1)
            e = StepE_eV*real((ebin-1),kind = 8) + StepE_eV *(1.-(log(RelCumFreqDist(ebin))-log(RandomNum))/(log(RelCumFreqDist(ebin))-log(RelCumFreqDist(ebin-1))))
        endif                 
                
       ! call random_number(harvest = RandomNum)
        call random_number(harvest = RandomNum)
        rr = 3.*RandomNum	! random number
        if(rr.le.1.)then
          valley=1
        elseif(rr.le.2.)then  
          valley=2
        elseif(rr.le.3.)then
          valley=3
        endif
        valley = 1 		! this is for GaAs materials
        !      Initial wavevector
        k=semh(valley)*sqrt(e*(1.+ nonParabolPara(valley)*e))
        call random_number(harvest = RandomNum)
        fai=2.*pi*RandomNum
        call random_number(harvest = RandomNum)
        ct=1.-2.*RandomNum
        st=sqrt(1.-ct*ct)
        kx=k*st*cos(fai)
        ky=k*st*sin(fai)
        kz=k*ct
        !      Initial free-flight
103     call random_number(harvest = RandomNum)
        if(RandomNum.le.1.e-5)go to 103
        tc = -(log(RandomNum))*tau_max(valley) ! random free flight time.. 
        !      Map particle atributes
        Par_wvk(i,1) = kx             
        Par_wvk(i,2) = ky
        Par_wvk(i,3) = kz
        Par_tau(i) = tc
        Par_valley(i) = valley
        Par_ene(i) = e
        
        

! Test
!        sq_wvk_sum = 0.
!        do j = 1,DIM !col_len = RCutoff_SiIn        
!            Temp_wvk(j) = Par_wvk(i,j)
!            sq_wvk(j) = Temp_wvk(j)*Temp_wvk(j)        
!            sq_wvk_sum = sq_wvk_sum + sq_wvk(j)
!            ! need to consider pos
!        enddo
!        abs_wvk = sqrt(sq_wvk_sum)
!        gk = hhem(valley)*sq_wvk_sum
!        Temp_ene = 2*gk/(1.+sqrt(1.+nonParabolPara4(valley)*gk))
!        if (Temp_ene.ne.e) then
!                err1 = 1.0
!        endif
        
        
        if (Par_ene(i).ge.1) then
            err1 = 1.0
        endif
    enddo
    return
end


subroutine initial_energy()
    use scattering  
    use simulation
    integer i,valley
    real ( kind = 8 )::totalFreq, NonparaCoeffG, D_eeV, fe, Temp_EeV  
    totalFreq = 0.
    valley = 1 ! initially located in Gamma
    CumFreqDist(:) = 0.
    do i = 1, num_EnergyLevel
        Temp_EeV = StepE_eV*real(i,kind = 8)
        NonparaCoeffG = (1+nonParabolPara2(valley)*Temp_EeV)*(sqrt(1+nonParabolPara(valley)*Temp_EeV))
        D_eeV = NonparaCoeffG*e_c*primV*sqrt(2*Temp_EeV*e_c)*eff_mass(valley)*sqrt(eff_mass(valley))/(pi*pi*hbar*hbar*hbar)
        exp_e = exp(((Temp_EeV-E_F)*e_c)/(kB*Te)); 
        fe = 1/(exp_e + 1);
        FreqDist(i) = fe*D_eeV*StepE_e;
        if (i.eq.1) then 
            CumFreqDist(i) = FreqDist(i)
        else
            CumFreqDist(i) = FreqDist(i)+CumFreqDist(i-1)
        endif
        totalFreq = totalFreq + FreqDist(i)
    enddo 
    do i = 1, num_EnergyLevel
        RelCumFreqDist(i) = CumFreqDist(i)/totalFreq
    enddo
    return
end