Instructions
Requirements and Specifications
Source Code
import numpy as np
import matplotlib.pyplot as plt
if __name__ == '__main__':
electron = np.loadtxt('electron_quantites.txt')
ion = np.loadtxt('ion_quantites.txt')
x = electron[:,0] # positions
rho = ion[:,1]
ki = electron[:,2]
Te = electron[:,1]
Ti = ion[:,2]
ke = ion[:,3]
n = ion[:,4]
u = ion[:,5]
p = ion[:,6]
E = ion[:,7]
# First, calculate du/dx
dudx = np.divide(u[1:]-u[0:-1], x[1:]-x[0:-1])
# Calculate dTi/dx
dTidx = np.divide(Ti[1:]-Ti[0:-1], x[1:]-x[0:-1])
# Calculate dTe/dx
dTedx = np.divide(Te[1:] - Te[0:-1], x[1:] - x[0:-1])
# From equation (1b), calculate p0u0^2 + p0
eq1b = np.multiply(p[:-1],np.power(u[:-1],2)) + p[:-1] - (4/3)*np.multiply(n[:-1], dudx)
# Equation (1c)
eq1c = np.multiply(E[:-1], u[:-1]) - np.multiply(ki[:-1], dTidx) - np.multiply(ke[:-1],dTedx) + np.multiply(p[:-1],u[:-1]) - (4/3)*np.multiply(n[:-1], np.multiply(u[:-1], dudx))
plt.plot(eq1b)
plt.plot(eq1c)
plt.show()