高分子科学専攻　助教　川口辰也

はじめに

日本学術振興会の小・中・高校生のためのプログラム「ひらめき★ときめきサイエンス」に参加して

阪大理学部の化学科は毎年このプログラムに参加して，8月の頭に高校生対象の一日体験入学を開催している。
筆者は2017年度と2018年度にX線回折の話で参加した。「原子で散乱された波の足し算を観測しているんだよ！」と高校生に実感して貰いたく… 波の足し算をするプログラムをC言語で作成した。

まあ，数十万円するソフトウェア（つまりMathematicaとか）ならすぐ式を書いて計算できるだろうし，無料のソフトでも探したら適当なものがあるだろう，たぶん。でも，痒いところに手が届かないだろうし，動作がそんなに速くないだろうし，適当なソフトを探して使い方を把握して段取りを決めて…と掛かる時間を考えると（この程度なら）作る方が早い。

結果できたプログラムソースを公開する。不具合が有っても責任は取れませんが好きに使ってください。拙いのはご容赦ください。

wave.c

どんなことをするプログラムか

XY平面をX方向に進行する平面波が，原点付近を中心とする（それぞれ一個は等方的に散乱する）散乱体の格子に入射したときに発生する散乱波の振幅を計算する。入射波は波長1.0の余弦波とする。入射波の位相，散乱体の格子の情報（X方向の個数と間隔，Y方向の個数と間隔，格子の回転角）と散乱振幅を計算するグリッドの情報（Xの下限・上限・刻み個数，Yの下限・上限・刻み個数）を与えることにより，グリッドの座標（X，Y）とそこでの散乱振幅（Z）のデータファイルとグリッドの座標（X，Y）とそこでの散乱体の形状（Z:半径は0.3とする）のデータファイルの2つを出力する。
データファイルをグラフにすると上の図のようなものを描ける。条件を少しずつ変えたグラフをたくさん描けば動画（2.1MB）も出来る。

必要とする環境

Mac，あるいはWindowsまたはLinuxの動くPC（64bit）

座標の計算だけでよい場合

WindowsとMacは各環境用に用意した計算プログラムのバイナリ（gcc -Ofast -march=core2 wave.c -lm : Windows版a.out，Mac版a.out）のみでできる。特にインストールしておかなければならないものはない。
もちろん，ソースプログラムを読んで，自分の環境に悪影響が無いのを確認してからコンパイルする方がよいが，その場合は適当なCのコンパイラが必要。
Linuxは当然自分でコンパイルしよう。

3次元グラフ化は各自の好きなソフトウェアを利用のこと（開発時にはgnuplotが念頭にあるが）。

シェルスクリプトを利用して波の動画を作る場合

bash
MacとLinuxでは（ほとんどの環境で）標準のシェルなので特に何もしなくても使える。WindowsはCygwinとかMSYS2とかまずはUnixコマンドを使える環境を整えよう。Windows10なら開発者モードでBash on Ubuntu on Windowsを使う手もあり。
gnuplot
グラフの作成についてはgnuplotに全て任せるようになっているので，WindowsとMacではgnuplotのインストールと環境変数にgnuplotへのパスの登録が必要である。Macの場合はgnuplotを使用するために予めXQuartz[https://www.xquartz.org]のインストールも必要である。
Linuxではgnuplotがインストールされてなかったら入れる。普通は入れてるね？
ffmpeg
動画・音声変換用ソフトとして有名でありインストール方法はすぐに検索できる。頑張れ。

自分でコンパイルする場合　　CのコンパイラはWindowsの場合は選択肢が多いので好みのものを，ただしOpenMPを利用したい場合は対応しているものを選択する。Macの場合はXcode（とCommand Line Tools）をインストールすればよいが，OpenMPをXcodeではサポートしない（2017年7月現在）ので，OpenMPを利用したい場合はGCCを別口でインストールする。また，オプティマイズはガンガン掛けましょう。ここに書いたことの意味が分からない人はそのあたりの知識がないということなので，この件については何もしないのが正解。もちろん勉強してできるようになるのを推奨する。なぜなら，計算量にも依るが，OpenMPを使うだけでDual CoreのCore i5で計算時間が半分以下に，Quad coreのCore i7で1/4以下ぐらいになるから。

使用法

Unixでの使い方を標準として記述する。Windowsの場合はコマンドの頭の「./」は省く

プログラムの設置

作業用のディレクトリを準備し，そこに計算プログラムのバイナリ「a.out」を置く

下拵え　Cのソースファイルから自分でバイナリを作成する場合

プログラムのソースファイルは「wave.c」である。一般的なC言語で記述してあるので適当なCコンパイラ（GCCを標準として考えてはいるが…）でコンパイルして実行ファイルを作成する。 GCCでコンパイルする場合は「gcc -O3 -fopenmp wave.c -lm」と入力すれば実行ファイル「a.out」ができる。MacでOpenMPを利用するためにGCCのバージョン6をインストールした場合は「gcc」の代わりに「gcc-6」とタイプすることになる。

命令

その前に記述の仕方の確認

端末を起動してbashが入力待ちの状態では，設定によって変わるが，「$」がプロンプトの前に表示されていることが多いようだ。そこで

で入力待ちの状態を，そして

$ ls
Concat_FFmpeg.sh  PCについて.txt   a.out            wave.c          使い方.txt
MakeGraph.sh      PCについて.txt~  template.gplt    wave.c~         端が入ったグラフ.png

で「ls」を入力し，ファイルのリスト(Concat_FFmpeg.sh以下)が表示されたところを示すことにする。ほぼ筆者の環境での端末での表示そのままである。gnuplotを起動したときは

$ gnuplot

    G N U P L O T
    Version 5.0 patchlevel 3    last modified 2016-02-21 

    Copyright (C) 1986-1993, 1998, 2004, 2007-2016
    Thomas Williams, Colin Kelley and many others

    gnuplot home:     http://www.gnuplot.info
    faq, bugs, etc:   type "help FAQ"
    immediate help:   type "help"  (plot window: hit 'h')

Terminal type set to 'wxt'
gnuplot>

となるので，「>」をgnuplotの入力待ちの状態を示すこととする。記述の仕方の確認は以上。

命令

$ ./a.out （ここに引数を並べる）   
    →  ファイル「out_ScrtWave.dat」（散乱波のデータ）と「out_Atom.dat」（散乱体のデータ）を書き出す

「a.out」に引数をつけて起動すると散乱波と散乱体の2つの座標データが出力される。

散乱波を計算するのに必要なパラメータは以下の項目を順番に並べる。ここにある値は例である。

位相情報（νtの値）nu_t	実数	0.0
散乱体の個数（X）n_x	整数	10
散乱体の間隔（X）s_x	実数	1.2
散乱体の個数（Y）n_y	整数	7
散乱体の間隔（Y）s_y	実数	1.0
散乱体格子の回転角（°）phi	実数	30.0
データ範囲（X下限）x_l	実数	-10.0
データ範囲（X上限）x_u	実数	30.0
データ範囲（Y下限）y_l	実数	-15.0
データ範囲（Y上限）y_u	実数	15.0
データ分割数（X）t_x	整数	800
データ分割数（Y）t_y	整数	600

上記のパラメータで計算をする場合は

$ ./a.out 0.0 10 1.2 7 1.0 30.0 -10.0 30.0 -15.0 15.0 800 600

これで「out_SctrWave.dat」と「out_Atom.dat」の2つのファイルが出力される。オプションとして文字列を追加すると出力ファイル名の「out」の部分が置き換わるので，

$ ./a.out 0.0 10 1.2 7 1.0 30.0 -10.0 30.0 -15.0 15.0 800 600 test

とすると「test_SctrWave.dat」と「test_Atom.dat」の2つのファイルが出力される。これらの「*_SctrWave.dat」および「*_Atom.dat」にある座標データをgnuplotなど，グラフ作成ソフトで表示させれば散乱体と散乱波の関係が分かりやすい。

データファイル

「*_Atom.dat」および「*_SctrWave.dat」は

X0 Y0 Z0
X1 Y1 Z1
X2 Y2 Z2
X3 Y3 Z3
X4 Y4 Z4
...

と座標の値がフォーマット%9.5f %9.5f %9.5fで書き込まれたテキストファイルである。

gnuplotでの表示

散乱体データの床レベルを-1.1にしているため，gnuplotではzrangeの下限を-1.0以上にしておくと良い。

> set zrange [-2.0:3.0]
> splot "out_Atom.dat" w l

とした場合と

> set zrange [-1.0:3.0]
> splot "out_Atom.dat" w l

とした場合の図を並べる。

> set autoscale
> set zrange [-1.0:3.0]
> set view 15.0, 20.0
> splot "out_Atom.dat" w l , "out_SctrWave.dat" w l lw 0.1

「set view」はグラフの方位・スケールの指定。何度でも設定できる。「set view 「rot_x」, 「rot_z」, 「scale」, 「scale_z」」と入力する。rot_z以降は省略可。{{{縦軸方向にscale_z倍してから}全体でscale倍し，}縦軸方向にrot_z回転してから}横軸方向にrot_x回転した状態で表示される。
「splot」は3次元グラフ表示の命令。引数としてファイル名あるいは関数，オプションを続ける。「w l lw 0.1」は「w[ith] l[ines] l[ine]w[idth] 0.1」の省略形。線でデータ点を繋ぎ，線の太さを0.1 pointとする。「l」以外に「l[ines]p[oints]」や「p[oints]」，「lw 数値」以外に「l[ine]c[olor] 整数値」，「p[oint]t[ype] 整数値」，「p[oint]s[ize] 数値」等が指定できる。

綺麗な図にするために一手間掛ける

上のデータを用いて実際にgnuplotで線を使って表示すると端のデータ点を繋ぐ線が入ってしまう。その不具合を解消するために，xの範囲が-10.0～30.0，分割数が800とすると，データの間隔が0.05になるので，それを元に範囲の両端に二つずつデータ点を増やして分割数に4 (2+2)を加える。つまり範囲を-10.1～30.1にして分割数を804とする。yも同様にして

$ ./a.out 0.0 10 1.2 7 1.0 30.0 -10.1 30.1 -15.1 15.1 804 604

としてデータファイルを作成する。そしてgnuplotではxrangeとyrangeをデータの両端を切るように指定する。

> set xrange [-10.0:30.0]
> set yrange [-15.0:15.0]
> splot "out_Atom.dat" w l , "out_SctrWave.dat" w l lw 0.1

上の図と比較すると分かるように端のデータ点を繋げる線が入ることを防ぐことができる。データの計算の入力の時点で素直でない数値を入れることになるが，この方法でgnuplot上で綺麗に表示が出来る。

バッチによるデータ処理

例えばnu_tを少しずつ変えた1周期分のデータからgnuplotで画像ファイルを作成し，パラパラ漫画の要領でffmpegを使って動画を作成することを考える。一つ一つコマンドを入力するという手間が掛かることをする気は起こらない。このようなものはバッチ処理にするのは当然である。以下のようなスクリプトを作成して1°ずつ散乱体の格子の向きを変えた計算とグラフの作成，動画の作成を一気にさせることができる。WindowsではBATファイルで同様のことができるはずであるが…筆者はDOSに疎いのでご勘弁ください（筆者であればCygwinかMsys2を入れてこのスクリプトを走らせる）。

シェルスクリプト MakeGraph.sh パラメータは適当に編集して使うこと

#!/bin/bash

NU_T=0.0
N_X=10
S_X=2.4
N_Y=10
S_Y=1.5
PHI=0.0
X_L=-20.2
X_U=40.2
Y_L=-30.2
Y_U=30.2
T_X=604
T_Y=604

j=0
endj=90

while [ $j -le $endj ]; do

    PHI=`echo $j \+ 0.0 |bc`
    GRAPH=$(printf "%dx%d-%04.1fdegree" $N_X $N_Y $PHI)

    i=0
    endi=19

    while [ $i -le $endi ]; do

        NU_T=`echo $i \* 0.05 |bc`
        COMMAND=$(printf "./a.out %f %d %f %d %f %f %f %f %f %f %d %d out" $NU_T $N_X $S_X $N_Y $S_Y $PHI $X_L $X_U $Y_L $Y_U $T_X $T_Y)
        $COMMAND
        COMMAND2=$(printf "sed s/OUTFILE/out%02d/" $i)
        COMMAND3=$(printf "sed s/INFILE/out/g")
        COMMAND4=$(printf "sed s/GRAPHTITLE/%s/" "$GRAPH")
        cat template.gplt | $COMMAND2 > TEMP1.gplt
        cat TEMP1.gplt | $COMMAND3 > TEMP2.gplt
        cat TEMP2.gplt | $COMMAND4 > COM.gplt
        gnuplot <COM.gplt
        rm TEMP?.gplt COM.gplt
        i=$(($i+1))
    done

    IPHI=$(echo $PHI |sed 's/.0$//')
    OUTFILENAME=$(printf "Wave%dx%d-%02d.mp4" $N_X $N_Y $IPHI)
    COMMAND0='ffmpeg -r 24 -i out%02d.png -an -c:v libx264 -vf format=yuv420p -tune animation '
    COMMANDV=$(echo "$COMMAND0" $OUTFILENAME)
    $COMMANDV
    rm out??.png

    j=$(($j+1))

done

gnuplotスクリプト template.gplt これもパラメータは適当に編集すること
MakeGraph.shにより出力するファイル名が置換されたものがgnuplotに読み込まれる。

set view 15.0, 20.0
set zrange [-1:2.7]
set xrange [-20.0:40.0]
set yrange [-30:30]
set title "GRAPHTITLE"
set output "OUTFILE.png"
set terminal pngcairo enhanced size 1200, 900
set key font 'Helvetica, 18' spacing 4
set title font 'Helvetica, 28'
splot "INFILE_Atom.dat" title "Atom" w l, "INFILE_SctrWave.dat" title "Scattering Wave" w l lw 0.1

上の例のスクリプトであれば，ターミナルで

$ bash MakeGraph.sh

とすれば，散乱体の格子（10x10，間隔2.4，1.5）が0°から90°まで1°ずつ回った散乱の様子を示すビデオ（Wave10x10-00.mp4からWave10x10-90.mp4まで）が作成される。それぞれの角度で波の1周期分になっている。

作成例

1x2，s_y=1.2の散乱体格子からの散乱波: 41.9MB
1x7, s_y=1.0の散乱体格子からの散乱波　近景:76.9MB，遠景:32.5MB
5x3，s_x=1.2，s_y=1.0の散乱体格子からの散乱波: 32.3MB
10x10, s_x=1.2, s_y=1.0の散乱体格子からの散乱波　近景:66.5MB，遠景:32.3MB
10x10, s_x=2.4, s_y=1.5の散乱体格子からの散乱波　近景:78.3MB，遠景:36.6MB

詳細

コンパイル

浮動小数点の計算でループが一杯で…最適化が露骨に効いてくるプログラムである。数値の厳密さをそれほど求めるものでもないので最適化は掛けられるだけ掛けるのが吉。

Ubuntu 16.04LTS，gcc 5.4.0，Core i7-6820HQの場合，最適化を掛けないと

$ gcc -O0 wave.c -lm
$ time ./a.out 0.0 15 1.2 15 1.0 30.0 -30. 40. -30. 30. 700 600

real    0m27.264s
user    0m27.236s
sys 0m0.020s

よく使う-O3だと

$ gcc -O3 wave.c -lm
$ time ./a.out 0.0 15 1.2 15 1.0 30.0 -30. 40. -30. 30. 700 600

real    0m10.417s
user    0m10.384s
sys 0m0.032s

最適化無しでOpenMPを有功にして

$ gcc -O0 -fopenmp wave.c -lm
$ time ./a.out 0.0 15 1.2 15 1.0 30.0 -30. 40. -30. 30. 700 600

real    0m4.991s
user    0m34.372s
sys 0m0.220s

-O3でOpenMPを有功にして

$ gcc -O3 -fopenmp wave.c -lm
$ time ./a.out 0.0 15 1.2 15 1.0 30.0 -30. 40. -30. 30. 700 600

real    0m2.316s
user    0m13.412s
sys 0m0.020s

OpenMPを無効として通常の最適化を掛けられるだけ掛けると

$ gcc -Ofast -march=native wave.c -lm
$ time ./a.out 0.0 15 1.2 15 1.0 30.0 -30. 40. -30. 30. 700 600

real    0m7.289s
user    0m7.252s
sys 0m0.032s

OpenMPを有効とし，最適化を掛けられるだけ掛けると

$ gcc -Ofast -march=native -fopenmp wave.c -lm
$ time ./a.out 0.0 15 1.2 15 1.0 30.0 -30. 40. -30. 30. 700 600

real    0m1.762s
user    0m9.092s
sys 0m0.012s

「real」のところの値が端末で入力後にコマンドプロンプトが戻ってくるまでの時間になる。

表にまとめると

OpenMP	-O0	-O3	-Ofast
Off	27.3秒	10.4秒	7.3秒
On	5.0秒	2.3秒	1.8秒

1.762秒の「-Ofast -march=native -fopenmp」が最速である。最適化無しに比べるとだいたい15倍速い。

波の計算

位置P（p,q）にある散乱体から発生する散乱波を考える。入射波はX方向に進む平面波

cos[2π（νt - x/λ）]

Pにおける入射波は

cos[2π（νt - p/λ）]

散乱波はπだけ位相が変わるので，Pで発生した散乱波はPにおいて

cos[2π（νt - p/λ - 1/2）]

位置（x, y）における散乱波はPとの距離をr = sqrt（（x-p）*（x-p）+（y-q）*（y-q））として

cos[2π（νt - p/λ - 1/2 - r/λ）]

ただし，球面波なので散乱体から離れるにつれて振幅が（1/距離:2次元のため）で減衰するため

1/r * cos[2π（νt - p/λ - 1/2 - r/λ）]

となる。各散乱体に対してこの式を作成し，全てを足し合わせると散乱体の格子からの散乱波となる。

散乱体の行列

このルーチンが一番面倒だった。原点対称にするために，原点から（散乱体の個数／2）個並べていく形にした。ソースプログラムを見てもらえば分かるだろう。
もっとエレガントなやり方を誰か教えてください。

散乱波の足し算: Total Sum of Scattering Waves.