rendered paste body 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