RASPフォーマット（.engファイル）の解析

showEng.py

#!/usr/bin/env python
import pandas as pd
import matplotlib.pyplot as plt
from scipy import integrate
import sys
import warnings
warnings.filterwarnings("ignore")

filename = sys.argv[1]
metadata = pd.read_csv(filename, delimiter='\s').columns.to_list()
df = pd.read_csv(filename, delimiter='\s', header=None)[2:].iloc[:,0:2].dropna()
df = df.rename(columns={0:'time', 1:'thrust'}).astype('float')
motor_name = metadata[0]
diameter = metadata[1]
length = metadata[2]
delays = metadata[3]
prop_weight = metadata[4]
total_weight = metadata[5]
manufacturer = metadata[6]

max_thrust = df['thrust'].max()
burn_time = df['time'].iloc[-1]
thrusts = df['thrust'].to_list()
times = df['time'].to_list()
total_impulse = integrate.cumtrapz(thrusts, times, initial=0)[-1]
ave_thrust = total_impulse / burn_time

print("(Metadata)")
print(f"motor name: {motor_name}")
print(f"diameter: {diameter} [mm]")
print(f"length: {length} [mm]")
print(f"delays: {delays} [s]")
print(f"propellant weight: {prop_weight} [kg]")
print(f"total weight: {total_weight} [kg]")
print(f"manufacturer: {manufacturer}")
print("(Calculated result)")
print(f"max thrust: {max_thrust} [N]")
print(f"burn time: {burn_time} [s]")
print(f"total impulse: {total_impulse} [N.s]")
print(f"average thrust: {ave_thrust} [N]")

plt.plot(times,thrusts)
plt.title('Thrust curve')
plt.xlabel('Time [sec]')
plt.ylabel('Thrust [N]')
plt.show()

追記：コメントを含む.engファイルに対応し、また複数プロットを重ねられるように修正を行った。(2022/6/23)

#!/usr/bin/env python
import pandas as pd
import matplotlib.pyplot as plt
from scipy import integrate
import sys
import warnings
warnings.filterwarnings("ignore")

class bcolors:
    HEADER = '\033[95m'
    OKBLUE = '\033[94m'
    OKCYAN = '\033[96m'
    OKGREEN = '\033[92m'
    WARNING = '\033[93m'
    FAIL = '\033[91m'
    ENDC = '\033[0m'
    BOLD = '\033[1m'
    UNDERLINE = '\033[4m'

def analyze(filename):
    lines = open(filename).readlines()
    comment_lines = 0
    metadata = ""
    for line in lines:
        comment_lines += 1
        if line[0] != ';':
            metadata = line
            break
    metadata = metadata.split() 
    df = pd.read_csv(filename, delimiter='\s', header=None, index_col=False)[comment_lines:].iloc[:,0:2].dropna()
    other_comment_lines = list(df[df[0]==';'].index)
    if len(other_comment_lines) >0:
         df = df[:other_comment_lines[0]-comment_lines]
    df = df.rename(columns={0:'time', 1:'thrust'}).astype('float')

    max_thrust = df['thrust'].max()
    burn_time = df['time'].iloc[-1]
    thrusts = df['thrust'].to_list()
    times = df['time'].to_list()

    motor_name = metadata[0]
    diameter = metadata[1]
    length = metadata[2]
    delays = metadata[3]
    prop_weight = metadata[4] 
    total_weight = metadata[5] 
    manufacturer = metadata[6] if len(metadata)<6 else ""
    
    total_impulse = integrate.cumtrapz(thrusts, times, initial=0)[-1]
    ave_thrust = total_impulse / burn_time

    print(f"{bcolors.WARNING}Motor name: {motor_name}{bcolors.ENDC}")
    print(f"    diameter: {diameter} [mm]")
    print(f"    length: {length} [mm]")
    print(f"    delays: {delays} [s]")
    print(f"    propellant weight: {prop_weight} [kg]")
    print(f"    total weight: {total_weight} [kg]")
    print(f"    manufacturer: {manufacturer}")
    print(f"    {bcolors.OKGREEN}max thrust: {max_thrust} [N]{bcolors.ENDC}")
    print(f"    {bcolors.OKGREEN}burn time: {burn_time} [s]{bcolors.ENDC}")
    print(f"    {bcolors.OKGREEN}total impulse: {total_impulse} [N.s]{bcolors.ENDC}")
    print(f"    {bcolors.OKGREEN}average thrust: {ave_thrust} [N]{bcolors.ENDC}")

    return (times, thrusts, motor_name)


def plot_thrust(times, thrusts, name):
    plt.plot(times,thrusts)
    plt.legend([name])
    plt.title('Thrust curve')
    plt.xlabel('Time [sec]')
    plt.ylabel('Thrust [N]')
    plt.show()

def plot_multiple_thrust(times1, thrusts1, name1, times2, thrusts2, name2):
    plt.plot(times1,thrusts1)
    plt.plot(times2,thrusts2)
    plt.legend([name1,name2])
    plt.title('Thrust curve')
    plt.xlabel('Time [sec]')
    plt.ylabel('Thrust [N]')
    plt.show()

if __name__=="__main__":
    num_files = len(sys.argv) - 1
    if num_files == 1:
        # only one file
        file1 = sys.argv[1]
        times, thrusts, name = analyze(file1)
        plot_thrust(times, thrusts, name)
    elif num_files == 2:
        # two files
        file1 = sys.argv[1]
        file2 = sys.argv[2]
        times1, thrusts1, name1 = analyze(file1)
        times2, thrusts2, name2  = analyze(file2)
        plot_multiple_thrust(times1, thrusts1, name1, times2, thrusts2, name2)