NASA CEAをMac OSXで実行

NASAのCEAをコマンドで実行するために必要な操作。
Windowsは比較的簡単だけど、Macはちょっと面倒。
試したのはMac book air,OSX yosemite
CEAって何かわかってる人向け。

英語読める人は下のPDFやyoutubeを参考にしてください。全部書いてあります。
http://spase.stanford.edu/CEA_files/cea_how_to_gfortran_snowleo.pdf
https://www.youtube.com/watch?v=CGpxfxmqA-c

gfortran

XCodeとCommand Line ToolsをインストールしてGCCを使えるようにする。
これでgfortranというコマンドでFortranのコンパイルが出来るようになる。
http://qiita.com/ShinichiOkayama/items/962e55ac7f8920a5f007

CEAのダウンロードとパス

NASAグレン研究センターのこのサイトに行く。
http://www.grc.nasa.gov/WWW/CEAWeb/ceaguiDownload-unix.htm
ここから下記ファイルをダウンロードして、インストールしたいフォルダに移動させる。

• CEAgui JAR (CEAgui-jar.tar.Z)
• CEA+Fortran Package (CEA+Fortran.tar.Z)
• CEAexec Package (CEAexec-mac.tar.Z)

私は~/Applications/NASAcea というフォルダを作りました。

ここで、ターミナルを起動。ターミナル上で~/Applications/NASAcea に移動して、

zcat CEA+Fortran.tar.Z | tar xvf -
zcat CEAexec-mac.tar.Z | tar xvf -
zcat CEAgui-jar.tar.Z | tar xvf -

と入れて解凍。さらにファイルのパーミッション設定として下記を入れる。

chmod a+x b1b2b3
chmod a+x syntax
chmod a+x FCEA2
chmod a+x runCEA.sh

ここで、フォルダにパスを通して、実行出来るようにする。

echo 'export PATH=$HOME/Applications/NASAcea:$PATH' >> ~/.bash_profile
bash -login

他のユーザーも使う

他のユーザーも使えるようにするためには下記をやる。個人で使う場合は必要なし。

sudo vi /etc/paths

それでこのファイルの最後の行に下記を追加。

~/Applications/NASAcea

CEAファイルのアップデート

実行ファイルをアップデートする。ターミナルで下記を入れる。

gfortran cea2.f
mv a.out FCEA2
gfortran b1b2b3.f
mv a.out b1b2b3
gfortran syntax.f
mv a.out syntax

このあと、ライブラリの更新のために順番にコマンドを入れていく。（連続でコピペするのではなく、1行づつ入れていく）。
【追記】FCEA2にパスを通してあり、別のディレクトリで実行する際にはthrmo.inp,trans.inp,cea2.inpを実行するディレクトリに事前に入れておく必要がある。以下のコマンドを実行後に、trans.lib,thermo.lib,(cea2.lib)ファイルがバイナリファイルで読めない文字が書き込まれていれば正常。からファイルの場合は上手く言っていない。

FCEA2
trans
FCEA2
thermo
FCEA2
cea2

これで完了。

サンプル

試しに、ケロシン/LOXのロケットエンジンの性能を出してみる。
sample.inpというテキストファイルを作る。拡張子.inpが大事。
中身は下記。

sample.inp

problem    o/f=2.5,
    rocket  equilibrium  frozen  nfz=2 
  p,bar=100,
  pi/p=100,
react  
  fuel=RP-1 wt=100  t,k=298.15  
  oxid=O2(L) wt=100  t,k=90.17  
output  
    plot o/f isp aeat cf t 
end

これで、ターミナル上でsample.inpのあるフォルダに移動して、

FCEA2
sample

この2行目はinpファイルの名前の部分を入れる。同じフォルダにsample.outというファイルが出てくるのでこれが結果。

CEAの.inpファイルの作り方がわからない場合はCEAのオンライン版でやりたいこと書いてポチポチして最後に結果ファイルと一緒に.inpファイルも出てくるのでそれを参考にする。
https://cearun.grc.nasa.gov/

参考

NASA-RP1311という論文中の最後の方にexampleがある。
http://www.frad.t.u-tokyo.ac.jp/public/cea/doc/xRP-1311P2.pdf
例えばロケットのことに関してはEXAMPLE8~14にある。
コメントである#のあとのproblemからendまでの部分が.inpファイルにするべき部分。

スクリプト化

CEAで出てくる結果は古臭くまとめにくい構造しているので、まとめる必要がある。試しにpythonスクリプトでやってみる。
MatlabでのCEAインターフェイスを作っている人もいたりするが
https://github.com/PurdueH2Lab/MatlabCEA
まともにインターフェイスを作るより、スクリプトで直接叩いてあげたほうが正攻法だろう。
例えば、液体酸素・液体水素のロケットエンジンの性能を計算してみる。
適宜ファイル名とか.inpファイルの部分を変えると使える。
FCEA2と同じフォルダに以下の.pyファイルを入れて、ターミナルから.pyを実行。私はipythonでFCEA2のフォルダ（~/Applications/NASAcea）に移動して

run LH2LO2.py

と実行している。

LH2LO2.py

# -*- coding: utf-8 -*-
import sys
reload(sys)
# デフォルトの文字コードを変更する．
sys.setdefaultencoding('utf-8')

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from matplotlib.font_manager import FontProperties
from scipy import optimize
import csv
import os
import subprocess

# 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()

def cea (set_of, set_chamber):
	file_name = '00temp' # 一時的に作られるファイル
	
	# ---- .inpファイル作り ----
	input_file = file_name + '.inp'
	str = """problem    o/f=%.1f,
	    rocket  equilibrium  frozen  nfz=2 
	  p,bar = %d
	  pi/p= %d
	react  
	  fuel=H2(L) wt=100  t,k=20.27  
	  oxid=O2(L) wt=100  t,k=90.17  
	output  
	    plot p pi/p o/f isp ivac aeat cf t 
	end
	""" % (set_of, set_chamber, set_chamber)
	 
	f = open(input_file, 'w') # 書き込みモードで開く
	f.write(str) # 引数の文字列をファイルに書き込む
	f.close() # ファイルを閉じる

	# ---- FCEA2コマンド ----
	# FCEA2の場所は ~/Applications/NASAcea/FCEA2 にしているが適宜書き換える
	cmd = '~/Applications/NASAcea/FCEA2'
	p = subprocess.Popen(cmd, shell=True, stdin=subprocess.PIPE, stdout=subprocess.PIPE)
	p.communicate(file_name)

	# ---- .pltファイル読み取り ----
	output_file = file_name + '.plt'
	with open(output_file, 'r') as f:
	    reader = csv.reader(f, delimiter=' ')
	    header = next(reader)  # ヘッダーの読み飛ばし

	    index = 0
	    for row in reader:
	    	while row.count("") > 0:
	    		row.remove("")
	        # print row          # 1行づつlist表示
	        if index == 0:
	        	chamber = float(row[0])
	        elif index == 2:
	        	of = float(row[2])
	        	isp = float(row[3])
	        	ivac = float(row[4])
	        	aeat = float(row[5])
	        	cf = float(row[6])
	        index += 1
	    data = [chamber, of, isp, ivac, aeat, cf]

	# ファイル削除
	os.remove(input_file)
	os.remove(output_file)

	return data

if __name__ == '__main__': # main（ここからはじまる）
	data_file = '00output.csv'

	plt.ion()
	plt.figure()

	# ---- パラメータ指定 ----
	of_array = np.linspace(3.0, 6.0, 7)
	# chamber_array = np.linspace(5, 40, 30)
	chamber_array = np.linspace(10, 200, 30)
	data_array = np.zeros((len(chamber_array), 6))
	save_array = np.zeros([1,6]) # CSVに書き込むデータのストック場所
	# ---- コマンドをループさせる ----
	# この例では燃焼室圧力を横軸、O/Fを縦軸にとってPLOTしたいので、O/F→圧力の順でループ
	for j in of_array:
		k = 0
		for i in chamber_array:
			data = cea(j, i)
			print data
			data_array[k] = np.array(data)
			k+=1
		g = 9.80665
		press = data_array[:,0] / g # 燃焼室圧力 MPa
		of = data_array[:,1] # O/F
		Isp = data_array[:,2] / g # 海上比推力 sec
		Ivac = data_array[:,3] / g # 真空中比推力 sec
		Cf = data_array[:,4] # 推力係数 Cf
		# ---- PLOT ----　例えば燃焼室圧力と海上比推力をPLOT
		plt.plot(press, Isp, label='O/F=%.1f' % (j))
		# ----      ----
		save_array = np.vstack([save_array, data_array])

	# ---- CSVで結果を保存 ----
	header = u'chamber[MPa]\tof[-]\tisp[sec]\tivac[sec]\taeat[-]\tcf[-]'
	np.savetxt(data_file, save_array, fmt='%.3f', delimiter='\t', header = header)

	# ---- PLOTの設定 ----
	plt.xlim(0,20)
	plt.xlabel('Pressure (MPa)')
	plt.ylabel('Isp (sec)')
	plt.grid()
	plt.title('LOX/LH2 100%  equilibrium  frozen  nfz=2')
	plt.legend(loc='best', fontsize=12)

こんな結果が出る。