DO i = 1, nx
          fi(i,1:nv) = SQRT(mass) * exp(-mass * vi(1:nv)**2/2.0_num/T_i(i)) &
               / SQRT(2.0_num*pi*T_i(i))
          ve_drift = 0.0_num
          f0(1:nv) = exp(-(ve(1:nv)-ve_drift)**2 / 2.0_num/T_e(i)) &
               / SQRT(pi*2.0_num*T_e(i))
          f1(1:nv) = (vp_EAW / (SQRT(2.0_num*pi)*T_e(i)**1.5_num))* &
               exp(-vp_EAW**2/2.0_num/T_e(i)) * &
               (ve(1:nv)-vp_EAW)*exp(-((ve(1:nv)-vp_EAW)**2.0)/du_EAW**2.0)
          !perterbation
          fp(1:nv) = -(a)*(sin(k_EAW*x(i))/((omega_EAW-k_EAW*ve(1:nv)))) *&
                                !due to f0
               (-ve(1:nv)*exp(-ve(1:nv)**2 / 2.0_num / T_e(i))/SQRT(2.0_num*pi)/T_e(i)**1.5_num +&
                                !due to f1
               vp_EAW*exp(-vp_EAW**2 / 2.0_num / T_e(i))/SQRT(2.0_num*pi)/T_e(i)**1.5_num * &
               exp(-(ve(1:nv)-vp_EAW)**2 / du_EAW**2) * &
               (1.0_num-2.0_num*(ve(1:nv)-vp_EAW)**2 / du_EAW**2))
          fe(i,1:nv) = f0(1:nv)+f1(1:nv)+fp(1:nv)
       END DO       DO i = 1, nx
          fi(i,1:nv) = SQRT(mass) * exp(-mass * vi(1:nv)**2/2.0_num/T_i(i)) &
               / SQRT(2.0_num*pi*T_i(i))
          ve_drift = 0.0_num
          f0(1:nv) = exp(-(ve(1:nv)-ve_drift)**2 / 2.0_num/T_e(i)) &
               / SQRT(pi*2.0_num*T_e(i))
          f1(1:nv) = (vp_EAW / (SQRT(2.0_num*pi)*T_e(i)**1.5_num))* &
               exp(-vp_EAW**2/2.0_num/T_e(i)) * &
               (ve(1:nv)-vp_EAW)*exp(-((ve(1:nv)-vp_EAW)**2.0)/du_EAW**2.0)
          !perterbation
          fp(1:nv) = -(a)*(sin(k_EAW*x(i))/((omega_EAW-k_EAW*ve(1:nv)))) *&
                                !due to f0
               (-ve(1:nv)*exp(-ve(1:nv)**2 / 2.0_num / T_e(i))/SQRT(2.0_num*pi)/T_e(i)**1.5_num +&
                                !due to f1
               vp_EAW*exp(-vp_EAW**2 / 2.0_num / T_e(i))/SQRT(2.0_num*pi)/T_e(i)**1.5_num * &
               exp(-(ve(1:nv)-vp_EAW)**2 / du_EAW**2) * &
               (1.0_num-2.0_num*(ve(1:nv)-vp_EAW)**2 / du_EAW**2))
          fe(i,1:nv) = f0(1:nv)+f1(1:nv)+fp(1:nv)
       END DO