Missile DATCOMのプロット

Missile Datcomの基礎の続き。

出力ファイルを可視化したくなったときに出力ファイルであるfor0042.csvファイルから可視化を行なうと効率的になる。
出力結果の例がネット上に落ちていたので、それを可視化してみる。

# -*- coding: utf-8 -*-
# Plotting script for Missile DATCOM.
# It is used output file "for042.csv"
# 2016/06/15　ina111
import sys
reload(sys)
import platform
# default encoding
sys.setdefaultencoding('utf-8')

import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.font_manager
from matplotlib.font_manager import FontProperties
from matplotlib.backends.backend_pdf import PdfPages

if 'Windows' == platform.system():
    fp = FontProperties(fname=r'C:\WINDOWS\Fonts\MSGothic.ttf')

if 'Darwin' == platform.system(): # for Mac
    font_path = '/Library/Fonts/Osaka.ttf'
    font_prop = matplotlib.font_manager.FontProperties(fname=font_path)
    matplotlib.rcParams['font.family'] = font_prop.get_name()
    matplotlib.rcParams['pdf.fonttype'] = 42

plt.ion()
plt.close()
plt.style.use('ggplot')

# output PDF file
pdf = PdfPages('plot.pdf')

input_file = u"for005.dat"
output_file = u"for042.csv"

# Reading CSV file for dateframe
df = pd.read_csv(output_file, skiprows=1)

# Clean column name
temp = list(df.columns)
for (i,x) in enumerate(temp):
    temp[i] = temp[i].lstrip().replace("'","")
df.columns = temp

# Grouping
grouped = df.groupby("ALPHA", sort=False)

# making X-CP_NOSE that is position of C.P. from nose top.(aft-)
df["X-CP_NOSE"] = -(df["XMRP"] - df["X-C.P."]) * 100

def plot_MD(key,label):
    plt.figure()
    keylist = sorted(grouped.groups.keys())
    for selected in keylist:
        df_temp = grouped.get_group(selected)
        # df_temp.plot(x="MACH", y=key, label="alpha=%d deg" % (selected))
        plt.plot(df_temp["MACH"], df_temp[key], label="alpha=%d deg" % (selected))
    plt.xlabel(u"Mach number")
    plt.ylabel(label)
    plt.legend(loc="best")
    # plt.savefig("plot_" + key.replace(".","") + ".png")
    pdf.savefig()

plot_MD("X-C.P.", u"重心からのX-C.P. (-)")
plot_MD("X-CP_NOSE", u"ノーズからのX-C.P. (%)");plt.ylim(ymax=0)
plot_MD("CN", u"法戦力係数 C_N")
plot_MD("CL", u"揚力係数 C_L")
plot_MD("CA", u"軸力係数 C_A")
plot_MD("CD", u"抗力係数 C_D")
pdf.close()

# Dateframe INDEX
# ['CASE', 'RAD?', 'TRIM?', 'MACH', 'RE', 'ALT', 'Q',
# 'BETA', 'PHI', 'SREF', 'XCG', 'XMRP', 'LREF', 'LATREF',
# 'ALPHA', 'CN', 'CM', 'CA', 'CA_0B', 'CA_FB', 'CY', 'CLN',
# 'CLL', 'CL', 'CD', 'CL/CD', 'X-C.P.', 'CNA', 'CMA',
# 'CYB', 'CLNB', 'CLLB', 'Unnamed: 32']