高分子科学専攻　　講師　川口辰也　(revised 2026.03.25)

Cプログラムの改定とページの微修正

ChatGPTとGeminiにプログラムの改良を相談してみると，まあ大体同じことを指摘されたので修正してみると...
散乱体の行列というか，並びを作る部分が簡潔に。
あれ，そうすると格子ベクトルのなす角を90度に固定にせずにパラメータとして入れるのも簡単やね。じゃ，入力部分に角度を入れようか。
他にもループ内での計算を減らすとか構造体の配列ではなく配列の構造体にしたほうが速く動くなどなど色々なアドバイスが。 手を加えだしたら面白くなってきてしまった。

はじめに

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

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

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

結果できたプログラムソース（改定前）...を生成AIのアドバイスで少しブラッシュアップしたものを公開する。 不具合が有っても責任は取れませんが好きに使ってください。拙いのはご容赦ください。

• wave_2.5.c

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

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

必要とする環境

Mac，あるいはWindowsまたはLinuxの動くPC（64bit）とそれぞれにおける開発環境（C言語）
古いバージョンではコンパイルしたバイナリを用意してありましたが，今回はソースのみの提供ですので...

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

• Windows，Mac，Linuxすべて共通してCのコンパイラがあれば取り敢えずOK。ただし，OpenMPを利用しないと計算にまあそこそこ時間が掛かる。
  開発環境はありがたくも全て無料で揃うが...ストレージの容量は喰う。
  Windowsはwsl2でUbuntuでも入れてGCCとOpenMPのライブラリをインストールすればOKでしょう。
  LinuxはGCCとOpenMP（以下略）。
  MacはXcode（Command Line Tools）を入れて，Homebrewを入れて，HomebrewのOpenMPのライブラリ"libomp"を入れればOK。
  XcodeのコンパイラはOpenMPに対応しているのだそうな（ただ，ライブラリは入っていないのでHomebrewから引用。HomebrewでGCCをがっつり入れる必要はない）。 まあOpenMP無しでも計算規模によって極端に遅くなるわけでもないのでXcodeだけでもOK。

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

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

• bash
  Windowsはwsl2でUbuntuを入れたならそのまま使える。
  Macの標準シェルはzshになったけど，ここのシェルスクリプトはzshでもそのまま使えるからbashに入れ替える必要はない。
  Linuxでは（ほとんどの環境で）標準のシェルなので特に何もしなくても使える。

• gnuplot
  グラフの作成についてはgnuplotに全て任せるようになっているので，WindowsとMacではgnuplotのインストールが必要である。
  Linuxではgnuplotがインストールされてなかったら入れる。普通は入れてるね？

• ffmpeg
  動画・音声変換用ソフトとして有名でありインストール方法は検索ですぐ見つかる。頑張れ。

使用法

Unixでの使い方を標準として記述する。

プログラムの設置

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

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

プログラムのソースファイルは「wave_2.5.c」である。一般的なC言語で記述してあるので適当なCコンパイラ（GCCを標準として考えてはいるが…）でコンパイルして実行ファイルを作成する。 GCCでコンパイルする場合は「gcc -O3 -ffast-math -march=native -fopenmp wave_2.5.c -lm」と入力すれば実行ファイル「a.out」ができる。
MacでXcode + libompだと「cc -Xpreprocessor -fopenmp -lomp -I"$(brew --prefix libomp)/include" -L"$(brew --prefix libomp)/lib" wave_2.5.c」となるらしい。

命令

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

端末を起動して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つの座標データが出力される。

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

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

$ ./a.out 0.0 10 7 1.2 1.0 90.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 7 1.2 1.0 90.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 7 1.2 1.0 90.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°ずつ散乱体の格子の向きを変えた計算とグラフの作成，動画の作成を一気にさせることができる。bashやzshで使える。

シェルスクリプト MakeGraph2.sh パラメータは適当に編集して使うこと
書き換えたパラメータでいきなりフルに走らせず，まずは「endj=90」を「endj=0」とし，「rm out???.png」を「# rm out???.png」とコメントアウトする， 等のループを減らす・記録を残す手はずをして動作を確認するのが良い

#!/bin/bash

NU_T=0.0
N_a=10
N_b=10
L_a=2.4
L_b=1.5
GAMMA=90.0
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 %d %f %f %f %f %f %f %f %f %d %d out%02d" $NU_T $N_a $N_b $S_a $S_b $GAMMA $PHI $X_L $X_U $Y_L $Y_U $T_X $T_Y $i)
        $COMMAND
        COMMAND2=$(printf "sed s/OUTFILE/out%02d/" $i)
        COMMAND3=$(printf "sed s/INFILE/out%02d/g" $i)
        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-%03d.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 MakeGraph2.sh ⏎

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

作成例

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

詳細

コンパイル

浮動小数点の計算でループが一杯で…最適化が露骨に効いてくるプログラムである。
コンパイラが最適化しやすい書き方をするのも重要とのことで，今回はそこにも気を使っている。
数値の厳密さをそれほど求めるものでもないので最適化は掛けられるだけ掛けるのが吉。

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

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

real	0m4.929s
user	0m4.845s
sys	0m0.032s

よく使う-O3だと

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

real	0m3.658s
user	0m3.553s
sys	0m0.040s

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

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

real	0m1.483s
user	0m5.836s
sys	0m0.068s

-O3でOpenMPを有功にして

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

real	0m1.316s
user	0m4.544s
sys	0m0.052s

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

$ gcc -O3 -ffast-math -march=native wave_2.5.c -lm ⏎
$ time ./a.out 0.0 15 15 1.2 1.0 90.0 30.0 -30. 40. -30. 30. 700 600 ⏎

real	0m1.221s
user	0m1.125s
sys	0m0.036s

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

$ gcc -O3 -ffast-math -march=native -fopenmp wave_2.5.c -lm ⏎
$ time ./a.out 0.0 15 15 1.2 1.0 90.0 30.0 -30. 40. -30. 30. 700 600 ⏎

real	0m0.926s
user	0m1.496s
sys	0m0.049s

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

表にまとめると

0.9秒の「-O3 -ffast-math -march=native -fopenmp」が最速である。最適化無しに比べるとだいたい5倍速い。

ちなみに改定前のwave.c（バージョン1）のときは以下のとおり

なんとまあ，プログラムのブラッシュアップ(とGCCの差？)で最適化無しでの時間が27.3秒から4.9秒に短縮されてますよ！
Ubuntu 16.04LTS，gcc 5.4.0，Core i7-6820HQ

波の計算について

位置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/sqrt(距離):エネルギーは1/(距離)）するため

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

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

散乱体の行列

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

何のことはない，並べてから散乱体の行列の中心を原点に平行移動すれば良いだけだった。
しかもこの方法なら散乱体の並びが直交してなくても簡単に作れるではないか！

シェルスクリプトが並列化されてないのでは？

気付いてしまったね。こちらへ