Week 6 - Data analysis

PROGRAM

#Importing libraries
import matplotlib.pyplot as plt
import numpy as np
from numpy import trapz as trp

#Compatability check
try:
	f=open('engine_data.out','r')

except FileNotFoundError:
	print('File not recognized. Please provide a valid CONVERGE output file!')
	sys.exit()

#input
file = 'engine_data.out'
count = 0
count_header = 0
column_title = []
column_units = []
for line in open(file,'r'):
	if '#' not in line:
		count = count + 1
		columns = len(str.split(line))
	else:
		count_header = count_header +1
		if count_header == 3:
			line_content = str.split(line)
			for i in range(0,len(line_content)):
				column_title.append(line_content[i])
		elif count_header == 4:
			line_content = str.split(line)
			for i in range(0,len(line_content)):
				column_units.append(line_content[i])

if (count_header < 4 or count < 1 or columns < 8):
	print('Invalid data file')
else:
	data = np.zeros((count,columns,1))
	row = 0
	for line in open(file,'r'):
		if '#' not in line:
			column = 0
			data_line = str.split(line)
			while (column < columns):
				data[row,column] = float(data_line[column])
				column = column + 1
			row = row + 1

column_number = 1
while (column_number <= columns):
	plt.plot(data[:,column_number-1])
	plt.ylabel(column_title[column_number]+column_units[column_number])
	plt.savefig(column_title[column_number]+'.png')
	plt.figure()
	column_number = column_number + 1

x_column_number = 8
y_column_number = 2
x_data = np.array(data[:,x_column_number-1])
y_data = np.array(data[:,y_column_number-1])
plt.plot(x_data,y_data)
plt.xlabel(column_title[x_column_number]+column_units[x_column_number])
plt.ylabel(column_title[y_column_number]+column_units[y_column_number])
plt.savefig(column_title[y_column_number]+' VS '+column_title[x_column_number]+'.png')

area = trp(y = y_data.transpose(), x = x_data.transpose())
print('Area under pv diagram:',area,'(MJ)')
#Power in kw = Nm * RPM / 9550
RPM = 1500
power = area * 1e+06*RPM/(9550)
print('Power:',power,'(kw)')
fuel_per_stroke = 20*1e-06
fuel_per_sec = fuel_per_stroke*RPM*60
Spec_fuel_cons = fuel_per_sec/ (power * 1e+03)
print('Specific Fuel Consumption:',Spec_fuel_cons,'(g/J)')

PLOTS