素人が独学でやって書いた記事です。適宜参考文献リスト等も参照下さい。

対象の読者

対象は個人ゲーム製作やCG関係で、流体の精度よりまずは動くプログラムが欲しい、学びたいという方。なので今回は一番単純で簡単な動く流体力学プログラムをつくろうと言うのが目標だ。
最終的に

こんな感じのプログラムが動くところまでやる。

導入

まず2Dでも3Dでも流体力学で避けて通れないのはナビエストークスの式（NS式）([1]wikipedia より)

(画像はwikipediaより)
左辺 - 第1項 : 時間微分項、第2項 : 移流項（対流項）
右辺 - 第1項 : 圧力項、第2項 : 粘性項（拡散項）、第3項 : 外力項


ちなみにこれでも簡単なほうらしい。これだけだと未知数が1つおおいままなので、∇・v = 0 という連続の式を連立して解けということだ（ｖは速度ベクトル場）



(画像はwikipediaより)
三角形の逆みたいな記号があるけどなんだこれは？内積の掛け算で使う「・」があるけど、ベクトルの計算式とかあるのか？と、私は最初見た時いろいろわからないとこだらけで、結局理解するのに２ヶ月くらいかかった。

実際のとこ∇をよく理解せずプログラムを組んで動いてしまったので、この記事を書いていた当時は「∇とかわからなくても四則演算でプログラムは作れる」「今はとりあえず動くプログラムを作ることが最優先」とか思っていた。(やっぱそれじゃダメだろと思ったので補足記事書きました)

粒子法か格子法か

今回の記事では格子法についてのみ解説を行う。
流体計算は大きく粒子法と格子法の2つがあり、前者は流体を粒子として計算、後者は流体の入った空間を何百という格子で分割して各々の格子の圧力や速度を計算するというもの。

どうも格子法の方が簡単だということがわかった。(粒子法専門の方には大変失礼ですが・・・)
ちなみに粒子法もたまにゲームやCGに組み込まれているのを見る。Position Based Fluidsというものらしい。

格子法で流体力学の完成形
例
[13]http://www.youtube.com/watch?v=4Ov6tguRr8E

格子法では、2Dの場合、各々の格子の中の圧力、x速度、y速度を時間ステップごとに更新していくものだ。

格子の種類

さて格子法でやると決めた後はどの種類の格子にするかを決める必要がある。大きく分けて通常格子、コロケート格子とスタガード格子がある。

通常格子

コロケート格子

スタガード格子

スタガード格子の方が少し複雑で難しい。が、プログラムを作っている途中でわかったが最終的にこっちの方が楽になる。

というわけでスタガード格子でやることに決定！

これは東工大のサイトが参考になった。
[2.1]http://www.eto.titech.ac.jp/contents/sub03/chapter08-1.html

で、次にスタガード格子の速度と圧力をNS式でどのように離散化して時間ステップを進めるかということを考える。

全体像と変数

流体格子の分割は10*10格子とする。
壁を含めると全計算領域は12*12格子である。
必要な変数はまずは圧力変数 「p」 、これは12*12格子

次にｘ速度変数 「vx」 、これは13*12
逆にｙ速度変数 「vy」 は12*13
また後々、ダイバージェンス(∇・v=0)の計算で使う各格子での「発散」値を格納するデータ 「s」 も必要。これは圧力と同じく12*12
「p」「vx」「vy」の一時的な記憶変数として「p_after」「vx_after」「vy_after」もそれぞれ同じ配列要素で用意。

そして、以降話を分り易くするためにΔxやΔyは1.0として考える　（少し衝撃的な決断だが）なので以降の式にΔxΔyが省略されているがそこはご了承を。
まともにプログラムが動いてからΔxやΔyをいじってみればいい。またΔtは1.0にすると計算が発散するので、今は0.1～0.2くらいだと思ってもらえればいい。

フローチャート

フローチャートはこうなる

上から順にみていこう

①移流粘性項

最初に移流と粘性を分けて考える。

移流項

移流とは、流れる前の速度変数 vx vy をつかって一定時間流れた後の速度を求める、ということだ。(と思ってましたがこれはあくまで陽解法の説明です。陰解法では...略)
流れ後の速度変数をそれぞれvx_after   vy_afterとすると

（i,jは配列の要素、整数）

u=vx[i][j]
v=( vy[i-1][j] + vy[i][j] + vy[i-1][j+1] + vy[i][j+1] )/4

if (u>=0.0)&(v>=0.0){
	vx_after[i][j] = vx[i][j] - u*(vx[i][j] - vx[i-1][j])*Δt - v*(vx[i][j]-vx[i][j-1])*Δt
}

if (u<0 .0)&(v>=0.0){
	vx_after[i][j] = vx[i][j] - u*(vx[i+1][j] - vx[i][j])*Δt - v*(vx[i][j]-vx[i][j-1])*Δt
}

if (u>=0.0)&(v<0 .0){
	vx_after[i][j] = vx[i][j] - u*(vx[i][j] - vx[i-1][j])*Δt - v*(vx[i][j+1]-vx[i][j])*Δt
}

if (u<0.0)&(v<0.0){
	vx_after[i][j] = vx[i][j] - u*(vx[i+1][j] - vx[i][j])*Δt - v*(vx[i][j+1]-vx[i][j])*Δt
}

ここでuとvは、代入が行なわれる変数vx_after[i][j]の位置にあるｘ速度、ｙ速度を意味する。よってｕ=vx[i][j]でいいがｖは、上のスタガード格子の性質上、ｘ速度の定義点とｙ速度の定義点がずれているため
v=( vy[i-1][j] + vy[i][j] + vy[i-1][j+1] + vy[i][j+1] )/4

と、４つの近くの格子点の情報を足して４で割る必要がある。


ここまででvxの移流の話が終わり、次はvyの移流である。
上の計算でvxとvyを置き換えるだけでいいがｕとｖは

u=( vx[i][j-1] + vx[i+1][j-1] + vx[i+1][j] + vx[i][j] )/4
v=vy[i][j]
とする。
したがって

u=( vx[i][j-1] + vx[i+1][j-1] + vx[i+1][j] + vx[i][j] )/4
v=vy[i][j]

if (u>=0.0)&(v>=0.0){
	vy_after[i][j] = vy[i][j] - u*(vy[i][j] - vy[i-1][j])*Δt - v*(vy[i][j]-vy[i][j-1])*Δt
}

if (u<0 .0)&(v>=0.0){
	vy_after[i][j] = vy[i][j] - u*(vy[i+1][j] - vy[i][j])*Δt - v*(vy[i][j]-vy[i][j-1])*Δt
}

if (u>=0.0)&(v<0 .0){
	vy_after[i][j] = vy[i][j] - u*(vy[i][j] - vy[i-1][j])*Δt - v*(vy[i][j+1]-vy[i][j])*Δt
}

if (u<0.0)&(v<0.0){
	vy_after[i][j] = vy[i][j] - u*(vy[i+1][j] - vy[i][j])*Δt - v*(vy[i][j+1]-vy[i][j])*Δt
}

これを計算領域内すべてのi jで計算する。
全計算が終わったならvx vy にvx_after vy_afterを代入する。計算領域はイメージ図１で言う緑の部分と接するか重なるところだけ。
つまり最外層の１層だけは参照されるのみで代入は行なわれない。

またまた東工大のサイトが参考になった。
[2.2]http://www.eto.titech.ac.jp/contents/sub03/chapter03.html

移流のイメージ

移流のイメージとしては例えば以下の画像の黄緑の部分の速度が（1.0,1.0）で白の部分が（0.0,0.0）、話を分り易くするためにΔt = 1.0だとしたらこうなる。

しかし実際はΔtや速度の値が 0.148 とか0.009とか不規則な数であり、これにより移流計算するごとに速度の情報が崩れまくる。
特に一次風上差分法はとてもその形が崩れやすい。黄緑色の領域がじわじわ薄くなって広がってしまうイメージだ。



数値流体力学ではこのプロファイルの崩壊をいかに防いで移流計算するかということに、そのシミュレーション精度がかかっている。
精度を上げるためには、Cubic セミ・ラグランジアン法やCIP法を用いると良いようだ。
今回は一番簡単な一次風上差分を採用した。

このサイトの座標のところがすごく理解の助けになった。(リンク切れ)
[3]http://www.eco.zaq.jp/env_univ/euler.html

粘性項

次に粘性であるがこれはただ単に、１つの点の速度が周囲４つ点に拡散するという意味である。
どのくらいの割合で拡散するかはRe（レイノルズ数）で決まる。

vx_after[i][j]=vx[i][j]-1/Re*(vx[i+1][j]+vx[i][j+1]+vx[i-1][j]+vx[i][j-1])*Δt
vy_after[i][j]=vy[i][j]-1/Re*(vy[i+1][j]+vy[i][j+1]+vy[i-1][j]+vy[i][j-1])*Δt

全計算が終わったならvx vy にvx_after vy_afterを代入する。

②外力項

マウス操作等で流体に速度をつけたいのなら、ここで速度に加算する。

また、壁と接する速度定義点は本来 速度=0.0でないといけないはずだが、前行程の移流粘性項で近隣の速度定義点から影響を受け0でなくなってしまった定義点があった場合に、ここで0に戻す。

そして流出境界などもここで設定する。

またまた東工大のサイトが参考になった！
[2.3]http://www.eto.titech.ac.jp/contents/sub03/chapter08-4.html

③ダイバージェンス計算

次に速度からダイバージェンス＝発散量を求める。

「s」はダイバージェンスを格納する配列変数。

s[i][j]=( -vx[i][j] -vy[i][j] +vx[i+1][j] +vy[i][j+1] )/Δt

ここでは1格子の空間から出ていく(or入ってくる)水の量をすべてのi,jで計算している。



スタガード格子を採用したので、この矢印を足し引きすれば発散は簡単に求まる。

この発散量を求めていく意味なのだが、記事の上の方で∇・v = 0 を連立するという話につながる。
∇・v = 0 つまり、どのi,jでも s[i][j] が0.0になる必要があると言っている。

移流や粘性、外力項などで速度情報が何度も書き換えられたため、今のままでは∇・v = 0 を満たしていないと考えられる。これでは非圧縮の条件が満たされない。

これをいっきに解決するのが次の④である。

④圧力のポアソン方程式

ここが山場だった。
上で求めた各格子のダイバージェンスから圧力を算出したい。ダイバージェンスから圧力を求めるには、以下の反復法のSOR法を使うといい。

SOR法

p[i][j]=(1.0-ω)*p[i][j] + ω/4*(p[i-1][j] + p[i+1][j] + p[i][j-1] + p[i][j+1] - s[i][j])

ωは加速係数といって自分の経験則では1.7～1.9で一番収束が早くなる。格子数が少なめなら1.7、多めなら1.9が良い(気がする)。
ΔxΔyを変えたい場合は
[2.5]http://www.eto.titech.ac.jp/contents/sub03/chapter07.html
の(7-4)式を参照のこと

これを計算領域内のi,jで、何回も反復することで収束していく。だいたい100回とか1000回とか？
当然初期p値が解から遠ければ、収束までの反復回数は増える。ここが一番計算コストのかかる部分なので、ここをいかに効率良く解くかで、全体の処理速度が決まってくるといってもいい。
反復回数は多分計算領域の大きさや許容できる誤差によるところもある。

さて、ここで２つ問題がある。

1.障害物の扱いについて

最外層以外にも任意の場所に障害物を置きたいことがある。ある流体格子の圧力を求めるとき隣の格子が壁の場合、SORの計算の前にその壁に自分の圧力をコピーしておく

薄茶色が壁、白が流体の存在するところである(さっきと色が逆ですみません)。

このとき圧力定義点の白丸はSOR法で代入も参照も行なわれる。
黄緑色の点は参照のみされる。
黒丸は参照すら行なわれない。

Bの圧力をSOR法で計算しようとする時、Bの上下左右の圧力が参照されるため、Cの圧力が参照されることになるが、Cは壁内にあるためここにはBの圧力を使う。よってBの圧力計算をする前に、CにBの圧力を代入しておく。次にAの圧力をSOR法で計算しようとする時、同じようにCにAの圧力を代入してから計算する。

結果的に黄緑色の部分の本来の圧力というものは無視される。まぁ壁の中に圧力なんてないわけなんだが。

以上の作業はSORの「ループ内」で行う必要がある。「最初に１回だけ壁に圧力をコピーしておいてあとはループ計算を流す」ではうまくいかないはず。

ここでも東工大のサイトが一番参考になった。
[2.6]http://www.eto.titech.ac.jp/contents/sub03/chapter08-2.html

この「壁」を考慮したやり方は他のサイトや書籍など調べてもあまり書いてないものの、発散から圧力を求めるときの連立一次方程式を考えると、このやり方が一番合理的だと思うので自分はこの方式を支持する。

さらにこの東工大の記事によると、格子の壁情報を格納する変数(trueなら壁 falseなら流体とか)を作ると便利だということだ。確かにそうすることで、今自分が壁なら圧力計算をskipという処理が書きやすい。

また余談になるが、SOR法は参照する４つの圧力情報のうち、計算後の情報と計算前の情報の２種類が混じり合っているように見えるがそれで良い。(ガウス・ザイデル法のアレと同じ話)

2.SOR法を何回反復するか

「p_after」に現在計算中の圧力情報を、「p」に1ループ前の圧力情報を格納しているとしてp[i][j] - p_after[i][j] の絶対値が、すべてのi,jで一定値以下なら「収束した」と判断する。

ただこんな作業を毎ループやっていると、計算量2倍になり非効率なので、10ループに1回にしたり、i,jを1つ飛ばしでやるなど工夫する。

また残差計算が面倒ならループ数を固定してしまっても良いと思われる。計算の最初のフレームだけ誤差が大きそうだから最初だけ1000回ループとして、それ以降はレンダリングやゲーム等のフレームレート維持のため10ループ固定、などという案も現実的だ。

⑤修正項

圧力により速度が修正される。ここで、修正されるというのは、∇・v = 0 になるように速度が調節されるということである。


スタガード格子の利点はここにも。圧力定義点がx速度点、y速度点を挟み込むように配置されているので修正はとても簡単だ。

vx[i][j]-=( p[i][j]-p[i-1][j] )*Δt
vy[i][j]-=( p[i][j]-p[i][j-1] )*Δt

⑤はこれだけ

これでやっと全ての処理が終わり。
vx vy p の値全てが次のタイムステップに進んだことになる。修正項のあとは、可視化ルーチンに回すなり、ランダムに配置した粒子の移動ルーチンに回すなり好きにすればいい。
そして処理は①へ

まとめ

こうして①～⑤を更新してSTEPを進めていくのだが、一発でプログラムが予想した挙動を示してくれるとは限らない。最後に注意すべき点がいくつかある。

まずコーディングミス。実装量的にそこそこな量になるし、一つでもミスると結果がアベコベな値になる。かなり慎重に書いていく必要がある。特にxとyの取違い

次に移流項。この記事で紹介した方法ではΔtを0.1～0.2と適度に細かくしておけば発散しない。（詳しくはCFL条件で検索）

次に圧力のポアソン方程式の項。圧力pを毎STEPごとに初期化せず、前STEPのを使いまわすほうが結果的に早く収束することが分かっている(warm startingとかいう手法だった気がする)。ただしその場合、何STEPか計算していくと数値誤差でpの値が正か負に発散していくことがある。pの値は前STEPのpにのみ依存するため誤差の蓄積が原因となる。pは結局その「xy方向差分」だけを計算に使うのでpが発散するまで気づかない。なのでpを使いまわすのであれば「ガス抜き」のような処理がどこかに必要。平均値が0.0になるような処理を。

もう少し突っ込んで話すと、x,yで周期境界条件を作ったとき、圧力の逃げ場がないので数値誤差が積み重なりpが正か負に発散する。境界条件にノイマン条件がどこか1か所でもあればそこが圧力の逃げ場になり発散はしない(と自分では理解している)。


ちなみに、この流体力学のサンプルソースを配布している。
コッチ(Python版)からどうぞ。

おしまい

とりあえずNS式の意味を深く考えず、ここをこうすればいいと計算の手順の部分だけ書いてきたが、一応これで動くはずである！正直、最初に流体力学をやろうと思ったときは「∇」とか「ベクトル場」とか「ラプラシアン」とか全く分からなかったが、見よう見まねでプログラムを動かそうとしているうちにだんだんと式や単語の意味が分かってきた。

私は、まずは動くプログラムを作ってから式や単語の意味を考えるのが一番効率がいいと思う。なのでとにかく動かすことを最優先で書いてきた。
2012年当時、ネットで参考になる資料はいくつかあったものの自分の理解できるレベルのものは少なく、東工大のサイトに出会わなかったらずっとCFDプログラムが作れずにいただろう。本当に感謝

参考文献

参考にしたサイトやらなにらや
[1]http://ja.wikipedia.org/wiki/%E3%83%8A%E3%83%93%E3%82%A8-%E3%82%B9%E3%83%88%E3%83%BC%E3%82%AF%E3%82%B9%E6%96%B9%E7%A8%8B%E5%BC%8F
[2]http://www.eto.titech.ac.jp/contents/
[3]http://www.eco.zaq.jp/env_univ/euler.html
[4]http://d.hatena.ne.jp/etopirika5/20110425/1303740299
[5]http://www.bee-www.com/log.htm
[6]http://journal.mycom.co.jp/articles/2007/10/09/cedec01/001.html
[7]http://www.4gamer.net/games/000/G000000/20070926041/screenshot.html?num=012
[8]http://fluid.me.seikei.ac.jp/lecture/cfd/cfd-08-NS-algorithm.pdf
[9]http://staff.aist.go.jp/kogaki.t/Dissertation/PDF_files/chap3.pdf
[10]http://homepage.mac.com/catincat/java/index.html
[11]http://oshiro.bpe.es.osaka-u.ac.jp/people/staff/imura/lecture/2007cg/fire071220.pdf
[12]http://www.youtube.com/watch?v=6rPL-QkUFf8&feature=related
[13]http://www.youtube.com/watch?v=4Ov6tguRr8E

他、個人ブログの方が書いたゲームetcに格子法流体を実装するのに参考になりそうなサイト(2020/6更新)
http://el-ement.com/blog/2015/03/16/gpu-fluid/
https://edom18.hateblo.jp/entry/2020/02/13/104744

2018/10追記
東工大の流体のページがいつからか閉鎖させておりアクセスできない状態になっていた。ただWARPという国内の教育機関などに特化したwebアーカイブサイトに保存されていた！！
見たい方はコチラで「コンピュータシミュレーションによる流体解析（発展編）」と検索
アーカイブ自体は私も保存しているのでもしほしい方はコメントを(教育目的のサイトだし著作権的に個人所有は問題ないはず・・)


次回：流体力学のプログラムを作りたい！その２

久しぶりに更新

夏休みに新しいノートパソコン買って、今まさにcore i5 450Mのスピードに感動しているところだ

前のノーパの８倍は早い！

最近暇ができてきたのでまたfortranでレイトレしてみた

パソコンも新しくしてスペックも上がり計算に２時間しかかからなかった

これもピクセルごとのRGBをfortranで出力したのをHSPで読み込んでスクリーン出力したものだ

相変わらず進歩しない風景だが一応ソースを載っけてみる

まずはfortranのソース↓（ちなみにコンパイルはg95っての使ってソースのファイル名はh.f90でした）

program IO
integer yoko,tate,kaiso,mate,jimenmate,i,k,m,mcnt1,htmate,mhtmate,hhpi,ii,kk,ll,h,jkl,jjk
real cb0,cb1,cb2,ca0,ca1,ca2,hsa0,hsa1,hsa2,hsb0,hsb1,hsb2,q0,q1,q2,qt2
real t,w,b,c,d4,vmr,vmg,vmb,vmi,kai,x,z,drtga,naiseki,light0,light1,light2
integer mtset(3,640,480)
real g(3,3)
real,allocatable,dimension(:,:) :: kyu
real,allocatable,dimension(:,:) :: ja
real,allocatable,dimension(:,:) :: jb1
real,allocatable,dimension(:,:) :: jb2
real,allocatable,dimension(:,:) :: j12
integer,allocatable,dimension(:,:,:) :: pset

print *,"grnd"
read *, jkl
allocate(ja(3,jkl))
allocate(jb1(3,jkl))
allocate(jb2(3,jkl))
allocate(j12(2,jkl))
jimenmate=jkl

print *,"circle?"
read *, jkl
allocate(kyu(5,jkl))
mate=jkl

print *,"アンチエイリアスx?"
read *, jkl
allocate(pset(3,jkl*640,480*jkl))

print *, "kaiso"
read *, kaiso

yoko=640*jkl
tate=480*jkl
t=-3.0
k=1
w=0.0
jjk=0
do i=1,mate
kyu(4,i)=0.5
kyu(1,i)=w
kyu(2,i)=4.9*mod(i,2)+0.27
kyu(3,i)=t
kyu(5,i)=kyu(4,i)**2
if (mod(i,2).eq.1) then
w=w+0.6
jjk=jjk+1
if (jjk.eq.k) then
jjk=0
w=-0.30*k
k=k+1
t=t+0.99
endif
endif

enddo
w=0.0
t=0.0
jjk=0
k=0

do i=1,jimenmate
ja(1,i)=0.01*rand()*2001-10.0
ja(2,i)=0.01*rand()*601+3.5
ja(3,i)=0.01*rand()*1500+9.25
jb1(1,i)=0.01*rand()*1000-4.95
jb1(2,i)=0.01*rand()*1000-4.95
jb1(3,i)=0.01*rand()*1000-4.95
jb2(1,i)=0.01*rand()*1000-4.95
jb2(2,i)=0.01*rand()*1000-4.95
jb2(3,i)=0.01*rand()*1000-4.95
j12(1,i)=0.7+0.01*rand()*130
j12(2,i)=0.7+0.01*rand()*130
w=sqrt(jb1(1,i)**2.0+jb1(2,i)**2.0+jb1(3,i)**2.0)
jb1(1,i)=jb1(1,i)/w
jb1(2,i)=jb1(2,i)/w
jb1(3,i)=jb1(3,i)/w
w=sqrt(jb2(1,i)**2.0+jb2(2,i)**2.0+jb2(3,i)**2.0)
jb2(1,i)=jb2(1,i)/w
jb2(2,i)=jb2(2,i)/w
jb2(3,i)=jb2(3,i)/w
enddo

light0=1.0*rand()*100-30.0
light1=1.0+rand()*50+rand()*(3.0+rand()*180)
light2=1.0*rand()*100-50.0
w=sqrt(light0**2+light1**2+light2**2)
light0=light0/w
light1=light1/w
light2=light2/w

do k=1,yoko
mcnt1=k
do ii=1,tate
ca0=0.0
ca1=3.0
ca2=-7.0
cb0=-0.004*yoko/jkl+0.008*mcnt1/jkl
cb1=0.004*tate/jkl-0.008*ii/jkl
cb2=9.3
w=sqrt(cb0*cb0+cb1*cb1+cb2*cb2)
cb0=cb0/w
cb1=cb1/w
cb2=cb2/w

hsa0=0.0
hsa1=0.0
hsa2=0.0
hsb0=0.0
hsb1=0.0
hsb2=0.0
mhtmate=-1
vmi=1.0
m=0
do hhpi=1,kaiso
if (m.eq.0) then
htmate=-1
t=99999990.0
do i=1,mate
if (mhtmate.ne.i) then

q0=ca0-kyu(1,i)
q1=ca1-kyu(2,i)
q2=ca2-kyu(3,i)
qt2=q0*q0+q1*q1+q2*q2
b=q0*cb0+q1*cb1+q2*cb2
c=qt2-kyu(5,i)
d4=b*b-c
if (d4.ge.0.0) then
kai=-sqrt(d4)-b
endif
if ((d4.gt.0.0) .and. (kai.gt.0.0) .and. (t.gt.kai)) then
t=kai
htmate=i
endif

endif
enddo
do i=1,jimenmate
if (mhtmate.ne.(i+10000)) then

q0=ca0-ja(1,i)
q1=ca1-ja(2,i)
q2=ca2-ja(3,i)
g(1,1)=cb0
g(2,1)=jb1(1,i)
g(3,1)=jb2(1,i)
g(1,2)=cb1
g(2,2)=jb1(2,i)
g(3,2)=jb2(2,i)
g(1,3)=cb2
g(2,3)=jb1(3,i)
g(3,3)=jb2(3,i)
drtga=g(1,1)*g(2,2)*g(3,3)+g(1,2)*g(2,3)*g(3,1)+g(1,3)*g(2,1)*g(3,2)
drtga=drtga-g(1,1)*g(2,3)*g(3,2)-g(1,2)*g(2,1)*g(3,3)-g(1,3)*g(2,2)*g(3,1)
if (drtga.eq.0.0) then
else
kai=-q0*(g(1,2)*g(3,3)-g(1,3)*g(3,2))+q1*(g(1,1)*g(3,3)-g(1,3)*g(3,1))
kai=kai-q2*(g(1,1)*g(3,2)-g(1,2)*g(3,1))
kai=kai/drtga

if ((kai.lt.j12(1,i)).and.(kai.ge.0.0)) then
kai=q0*(g(1,2)*g(2,3)-g(1,3)*g(2,2))-q1*(g(1,1)*g(2,3)-g(1,3)*g(2,1))
kai=kai+q2*(g(1,1)*g(2,2)-g(1,2)*g(2,1))
kai=kai/drtga
if ((kai.lt.j12(2,i)).and.(kai.ge.0.0)) then
kai=q0*(g(2,2)*g(3,3)-g(2,3)*g(3,2))-q1*(g(2,1)*g(3,3)-g(2,3)*g(3,1))
kai=kai+q2*(g(2,1)*g(3,2)-g(2,2)*g(3,1))
kai=-kai/drtga
if ((kai.gt.0.0).and.(t.gt.kai)) then
t=kai
htmate=i+10000
endif
endif
endif
endif

endif
enddo
if (cb1.lt.0.0) then

x=cb0*ca1/abs(cb1)
z=cb2*ca1/abs(cb1)
kai=sqrt(x*x+ca1*ca1+z*z)
x=x+ca0
z=z+ca2

if (kai.lt.t) then
if (mod(abs(int(x*2.0)+int(z*2.0)),2) .eq. 1) then
vmr=10.0
vmg=10.0
vmb=93.0
else
vmr=155.0
vmg=123.0
vmb=190.0
endif
w=abs(x)+abs(z)
if (w.lt.400.0) then
vmr=vmr*(400.0-w)/400.0+3.0
vmb=vmb*(400.0-w)/400.0+4.0
vmg=vmg*(400.0-w)/400.0+11.0
else
vmr=3.0
vmb=4.0
vmg=11.0
endif
ca0=x
ca1=0.0
ca2=z
x=1.0
do kk=1,mate
q0=ca0-kyu(1,kk)
q1=-kyu(2,kk)
q2=ca2-kyu(3,kk)
qt2=q0*q0+q1*q1+q2*q2
b=q0*light0+q1*light1+q2*light2
c=qt2-kyu(5,kk)
if ((b*b-c).gt.(0.0)) then
x=0.0
endif
enddo

do kk=1,jimenmate
q0=ca0-ja(1,kk)
q1=-ja(2,kk)
q2=ca2-ja(3,kk)
g(1,1)=light0
g(2,1)=jb1(1,kk)
g(3,1)=jb2(1,kk)
g(1,2)=light1
g(2,2)=jb1(2,kk)
g(3,2)=jb2(2,kk)
g(1,3)=light2
g(2,3)=jb1(3,kk)
g(3,3)=jb2(3,kk)
drtga=g(1,1)*g(2,2)*g(3,3)+g(1,2)*g(2,3)*g(3,1)+g(1,3)*g(2,1)*g(3,2)
drtga=drtga-g(1,1)*g(2,3)*g(3,2)-g(1,2)*g(2,1)*g(3,3)-g(1,3)*g(2,2)*g(3,1)
if (drtga .ne. 0.0) then
kai=-q0*(g(1,2)*g(3,3)-g(1,3)*g(3,2))+q1*(g(1,1)*g(3,3)-g(1,3)*g(3,1))
kai=kai-q2*(g(1,1)*g(3,2)-g(1,2)*g(3,1))
kai=kai/drtga
if ((kai.lt.j12(1,kk)).and.(kai.ge.0.0)) then
kai=q0*(g(1,2)*g(2,3)-g(1,3)*g(2,2))-q1*(g(1,1)*g(2,3)-g(1,3)*g(2,1))
kai=kai+q2*(g(1,1)*g(2,2)-g(1,2)*g(2,1))
kai=kai/drtga
if ((kai.lt.j12(2,kk)).and.(kai.ge.0.0)) then
kai=q0*(g(2,2)*g(3,3)-g(2,3)*g(3,2))-q1*(g(2,1)*g(3,3)-g(2,3)*g(3,1))
kai=kai+q2*(g(2,1)*g(3,2)-g(2,2)*g(3,1))
if ((kai/drtga) .lt. 0.0) then
x=0.0
endif
endif
endif
endif
enddo

if (x.eq.1.0) then
naiseki=cb0*light0-cb1*light1+cb2*light2
vmr=vmr+(1.3+naiseki)*55.0
if (naiseki.gt.0.95) then
vmr=vmr+(naiseki-0.95)*1730.0
endif
vmg=vmg+(1.3+naiseki)*35.0
if (naiseki.gt.0.95) then
vmg=vmg+(naiseki-0.95)*900.0
endif
vmb=vmb+(1.3+naiseki)*17.16
if (naiseki.gt.0.95) then
vmb=vmb+(naiseki-0.95)*450.0
endif
else
vmr=vmr/1.7
vmg=vmr/1.7
vmb=vmr/1.7
endif
htmate=-1
endif

endif
mhtmate=htmate
if (htmate.eq.-1) then
m=1
else

if (htmate.lt.10000) then
hsb0=(ca0-kyu(1,htmate)+t*cb0)/kyu(4,htmate)
hsb1=(ca1-kyu(2,htmate)+t*cb1)/kyu(4,htmate)
hsb2=(ca2-kyu(3,htmate)+t*cb2)/kyu(4,htmate)
hsa0=kyu(1,htmate)+hsb0*kyu(4,htmate)
hsa1=kyu(2,htmate)+hsb1*kyu(4,htmate)
hsa2=kyu(3,htmate)+hsb2*kyu(4,htmate)
endif

if (htmate.ge.10000) then
htmate=htmate-10000
hsa0=ca0+t*cb0
hsa1=ca1+t*cb1
hsa2=ca2+t*cb2
hsb0=jb1(2,htmate)*jb2(3,htmate)-jb1(3,htmate)*jb2(2,htmate)
hsb1=jb1(3,htmate)*jb2(1,htmate)-jb1(1,htmate)*jb2(3,htmate)
hsb2=jb1(1,htmate)*jb2(2,htmate)-jb1(2,htmate)*jb2(1,htmate)
w=sqrt(hsb0*hsb0+hsb1*hsb1+hsb2*hsb2)
hsb0=hsb0/w
hsb1=hsb1/w
hsb2=hsb2/w
htmate=htmate+10000
endif

naiseki=-(hsb0*cb0+hsb1*cb1+hsb2*cb2)
ca0=hsa0
ca1=hsa1
ca2=hsa2
cb0=2.0*naiseki*hsb0+cb0
cb1=2.0*naiseki*hsb1+cb1
cb2=2.0*naiseki*hsb2+cb2
vmi=vmi+0.7

endif
endif

enddo

if (cb1.ge.0.0) then

naiseki=(cb0*light0+cb1*light1+cb2*light2)
vmr=(1.5+naiseki)*120.0
if (naiseki.gt.0.96) then
vmr=vmr+(naiseki-0.96)*8200.0
endif
vmg=(1.5+naiseki)*120.0
if (naiseki.gt.0.96) then
vmg=vmg+(naiseki-0.96)*8200.0
endif
vmb=(1.5+naiseki)*190.0
if (naiseki.gt.0.96) then
vmb=vmb+(naiseki-0.96)*8200.0
endif

endif

ll=int(vmr/vmi)
if (ll.lt.256) then
pset(3,mcnt1,ii)=ll
else
pset(3,mcnt1,ii)=255
endif

ll=int(vmg/vmi)
if (ll.lt.256) then
pset(2,mcnt1,ii)=ll
else
pset(2,mcnt1,ii)=255
endif

ll=int(vmb/vmi)
if (ll.lt.256) then
pset(1,mcnt1,ii)=ll
else
pset(1,mcnt1,ii)=255
endif

vmr=0.0
vmg=0.0
vmb=0.0
enddo
enddo

open(11,file = 'file1.txt')
open(12,file = 'file2.txt')
open(13,file = 'file3.txt')
open(14,file = 'file4.txt')
i=0
h=0
k=tate/jkl+1
m=0
ii=0
kk=1
ll=0
do i=1,yoko
do h=1,tate
do m=1,3
mtset(m,(i-1)/jkl+1,(h-1)/jkl+1)=mtset(m,(i-1)/jkl+1,(h-1)/jkl+1)+pset(m,i,h)
enddo
enddo
enddo
m=0
h=0
i=0

hhpi=yoko*tate/jkl/jkl*3/4
do i=1,yoko*tate/jkl/jkl*3
h=mod(i,3)+1
m=mod(i,4)+1
if (h.eq.1) then
ii=mod(ii,yoko/jkl)+1
if (ii.eq.1) then
k=k-1
endif
endif
ll=ll+(mtset(h,ii,k)/jkl/jkl)*kk
kk=kk*256
if (m.eq.4) then
write(11+(i-1)/hhpi,*) ll
ll=0
kk=1
endif
enddo
close(11)
close(12)
close(13)
close(14)
print *,"end111111111111111111111"
stop
end

私は実行するときコマンドプロンプトから実行しているが、多分これを実行すると最初に

grnd?

と入力を求められる

ここに「四角形の板」の枚数を入力する。鏡みたいなやつだ。当然反射する。

上の生成画像では確か１０って設定だったかな？

次に

circle?

と入力を求められる

ここに「球」の個数を入力する

上の画像では確か15000くらい。

アンチエイリアス?

にはアンチエイリアスなしで1を、2,3,4となるごとにきめ細かくなっていくが計算時間が２乗倍になる

上の画像では確か4

kaiso?

は、最高反射計算回数。まぁ4以上の数字ならたいてい変な画像はできない

上の画像は確か19とか

これでenterキーを押せば計算開始、計算が終了する頃に「end111111111111111111111」っていう計算終了のお知らせメッセージが表示される

また出力ファイルで４つのテキストファイルが出力される

そこに全ピクセルのRGB情報が詰まっている

今度はHSPでスクリーン出力

そのソースは↓

onexit*enf
repeat 4
ww=3-cnt
exist ""+ww+""
if strsize=-1:IIIIIIIIIIDDDDDDD=ww:bsave ""+ww+"",IIIIIIIIIIDDDDDDD:break
if cnt=3:end
loop

screen 0,640,480:boxf
dim i,640*480*3/4
mref i,66
sdim k,16

h=640*480*3/4/4*IIIIIIIIIIDDDDDDD
repeat 1,IIIIIIIIIIDDDDDDD
sdim a,1
notesel a
noteload "file"+(cnt+1)+".txt"
repeat 640*480*3/4/4
noteget k,cnt
i.h=int(k)
h+
if cnt\900=0:await 1:redraw
loop
loop

bmpsave ""+IIIIIIIIIIDDDDDDD+".bmp"
delete ""+IIIIIIIIIIDDDDDDD+""
if IIIIIIIIIIDDDDDDD=1{
repeat 999999
w=0
repeat 4
exist ""+cnt+""
if strsize!-1:w+
loop
if w=0:break
wait 200
loop
screen 0,640,480:boxf
repeat 4
buffer 1+cnt,640,480:picload ""+cnt+".bmp"
loop
gsel 0
repeat 4
pos 0,0
gmode 5,640,480,256
pos 0,0
gcopy 1+cnt,0,0,640,480
loop
bmpsave "5.bmp"
delete "0.bmp"
delete "1.bmp"
delete "2.bmp"
delete "3.bmp"
stop
}else{
end
}
*enf
delete ""+IIIIIIIIIIDDDDDDD+""
end

相変わらず超汚い

こんなの誰も理解できんだろうな・・そのうち自分も分からなくなるよこれ

実はこのプログラムDualコア4スレッド対応プログラム

でもHSPはシングルスレッドなのでやっぱ４回起動する必要があるｗ

１スレッドで普通にプログラム実行すると1分くらいかかるので折角のcore i5なので「並列プログラム」を組んでみたというわけ

注意点はなるべく素早く４回起動しないとバグる！(ぇ

しかし別に原理的にはただ単に４つ起動させているだけなので１コアのCPUでも全然動くｗ

で、４つ起動させて全部の計算が終わると残り３つのプログラムは自動終了して１つのスクリーンに完成画像が表示されて「5.bmp」ってファイルに出力される。

これでやっと上の画像が完成

おしまい

っと、今のところfortranとHSPの組み合わせ研究はこんな感じで進んでいる。

今度はfortranでムービーの作成でもやってみようかと考えている！

でRGBデータ→ムービーの変換のところでHSPを使うって感じで

Perlin noiseのことを最初適当に知ろうと思ってwikipediaをのぞいてみたが、なんかあんま参考にならん・・・ページの下に関連リンクがあったので覗いてみた。

生成方 

KEN 

(2019/12現在リンク切れ)

図やグラフから文章を推測するというなんとも奇妙な解読法で、なんとか理解したような気がした。

実際プログラムするときは上の生成方のＨＰが参考になった。

ではわかったことだけでもまとめていきたいと思う。

Perlin noise と普通のnoiseの違い　（導入）

以下の話は全てパーリンノイズではなくバリューノイズというやつでした。すみませんでした。パーリンノイズは近日中にコード載せます。

定義の話とか面倒なので具体的な話から・・

ノイズ関数とパーリンノイズ関数があったとしよう。

ノイズ関数・・・noise(a)

パーリンノイズ関数・・・Pnoise(a)

 

２つの関数は-100～100までのどれかの整数値を返すとして

aに０，１，２，３，４・・・と代入していき、それぞれの関数からどんな値が返ってくるか。

noise(a)はというと・・・

61，-79，24，-10，-31，95，7，56，-88，-12，8，-47

Pnoise(a)はというと・・

36，22，-19，-14，-45，-50，-38，12，52，59，59，53，69，89

ノイズ＝完全にランダム。ばらついてる。規則性がない。

パーリンノイズ＝規則性が見えないけどちょっと連続しているような・・・。ランダムに見えて不自然にばらついてるようなランダムではない。

a(0,1,2,3…)を横に、関数の返り値を縦に取ったグラフが下だ。
 

こうやって見ると一目瞭然である！
（もしや・・・Pnoise(a)は、グラフの微分係数をランダムな値にすればできるのでは？って思うかもしれないけどそうでもない！ま、具体的なPnoise(a)関数の生成方は後にして・・・）

パーリンノイズの特徴はなんといっても、山の起伏のような、自然なばらつきがある雰囲気を出せることなのだ。

このPnoise(a)関数を使えば年輪などの模様や雲のゆらぐようなテクスチャを作れるって言う、ちょっと想像もつかないような話なのだが、本当だ。

Perlin noiseの原理！

バリューノイズの原理！

まずいきなり完成系の画像をお見せしよう！

これは以下のように作る

一番左の画像=ぼかす前＆拡大前

画像×６枚は、単にランダムな濃さの点を打って作成したノイズ画像である。

ただし６枚それぞれが、ドットの数というかスケールが違うのだ。

上から256ﾄﾞｯﾄ四方、128、64、32、16とスケールが２分の１ずつ小さくなっている。

なので一番左の画像枚数のことをオクターブと呼ぶ。

今回は６枚なので６オクターブと言うことになる。普通５～６である。

あまり少なくしても、細部のスケールが表現できなかったりするし、大きすぎても処理時間がかかるだけなので・・・

ぼかす処理はガウスぼかし等で。そしてぼかし終わった後できた256ドット四方の６枚の画像を加算合成するのだが、普通加算するときにブレンド率にpersistence(振幅の変化)を掛けて加算する必要がある。

つまり左の自家製の汚いイメージ図で言うと、上のほうの画像を弱く、下のほうの画像を強くして加算合成することで完成系ができあがるのだ。

「最初一番小さくて一番拡大する画像ほどpersistenceは大きく、一番拡大しない画像ほどpersistenceを小さくして加算することが重要！」

具体的に言うと、拡大しない画像ほどブレンド率が２分の１ずつに減っていくのが理想の形！

上の画像は1.8分の１なので約２分の１と言える。

次回：パーリンノイズアルゴリズム 後編

レイトレ実験はその②ですが、レイトレの記事はもうこれで５つ目です・・なんかタイトルのつけ方ミスった？？

実は新しいプログラム言語　「fortran」　に手を出してみました！

これは計算に特化した言語でよく科学技術用の計算で用いられる、いわば四則演算がすごく速く行える言語です！

ＨＳＰでレイトレーシングをやるとすごく時間がかかるのでfortranに移植してみました。

当然「できあがる画像」はまったく一緒です！

同じアルゴリズムですから・・・

fortranで作ったＣＧ↓





 

こんなにたくさんのオブジェクト（球）があるとＨＳＰでは何時間計算かかるか分かりません

がfortranでは２０分くらいで作れました。

なんだかよくわからないグラフィックになってしまいましたが、球が規則正しく敷き詰められているだけです。

下と上に１層ずつ、って

見ててきづいたのですが、画面の球にはさまれた部分（空と地面）の形は、本来なら長方形になるはずが、中心の部分がふっくらして見えません？？

きっと目の錯覚と言うやつですね・・・・

どうでもいいです・・・・

ところでfortranに手を出したわけは・・

これからＨＳＰをメインに使わなくなると言うわけではなく、あくまでレイトレーシングの計算をＨＳＰで完成させたあと、fortranで動かすというカンジでやっていくためです。

やはりＨＳＰはデバックが簡単だしやスクリプト画面もみやすいので。

それに一番は、もうそれに慣れてしまったところが大きいですかね

で一番自分が気になる速度比較ですが、計算速度はなんとＨＳＰの３０倍くらい早かったです！！

こんど結果をグラフにしてみたいと思います。

しかしfortranで計算が早くなったのはいいが画面出力のやり方がわからないので、計算結果→画面出力をほとんどＨＳＰでやる羽目に・・

まずfortranでレイトレ演算結果の数値（ＲＧＢ）をtxtに出力させ、それをＨＳＰで読み込み、読み込んだ文字列を数値に変換しpsetで画面に打っていくという、なんともややこしく回りくどい方法を取らざるを得なくなり、結局その変換にすごい計算時間がかかるというオチであった・・orz

そんなこんなで４、５分くらいの計算ならＨＳＰで全部やったほうが結果的に早いことが判明

もしかしたらfortranだけで画面出力できる方法があるのかもしれないがあまり書籍もないし、自分の能力も限界に近いので、数値→画像変換はＨＳＰに頼ることにした。

やっぱわたしゃＨＳＰから離れられんなぁ～

ところで↑のＣＧ見て思った
ドットが荒い！

640×480の計算結果を640×480の画面に出力させてるから、１／１でドットが目立ち細かいところは何がなんだか分からなくなっている。

なので次に作る画像はアンチエイリアシングをかけてみようと思った。

擬似的なアンチエイリアシングの技法として目的の２～３倍程度の大きい画像を最初に作り、それを縮小してドットを滑らかにする、というものがある。

２倍の例だと

1280×960の計算結果を640×480の画面に出力させる、など

これは倍率を上げれば上げるほどドットが滑らかになるがその２乗に比例して演算時間がかかってしまう。

多分だが、本当のアンチエイリアシングは浮動小数点かなんか使って↑とは全く別のアルゴリズムでドットを滑らかにしているんだろうが

ま、レイトレーシングでアンチエイリアスするアルゴリズムなんて少し考えただけでも恐ろしく難しいので擬似的な方のアンチエイリアスでやることにした。

さて前回の記事で、四角板もレイトレーシングで描画できるのではと予言しましたが、やってみたら正解！

板の全反射なので「鏡」みたいな役割をするんじゃないかと思ってたら、その通りだった

以下のＣＧは３倍の倍率でやった擬似アンチエイリアシングの、四角形板＆球混合バージョン

（これは最初からＨＳＰで計算）




うぉぉ・・

これには感動した！！（泣）

四角形の板は、球より難しいアルゴリズムであったがちゃんとピカピカ光ってる！

アンチエイリアシングも思いのほか効きが良くてうれしいー！

うんしばらくこれをデスクトップの壁紙にしよう！

四角形の板のアルゴリズムは次で解説します。

あと途中から語尾が変わってたのはまぁ文才がないってことで・・m(_ _ )m