2018年07月16日

GMT講座12】pscoast その4

引き続いて、今回は異なる絵になります。
pscoastを使います。

円錐図法の一つのアルバート正積図法になるようです。
そのため、-Jオプションは
-JB経度/緯度1/緯度2/幅
を指定するようです。

とにかく実行してみると、
$ gmt pscoast -R-130/-70/24/52 -JB-100/35/33/45/6i -Ba -B+t"Conic Projection" -N1/thickest -N2/thinnest -A500 -Ggray -Wthinnest -P > GMT_tut_4.ps



となります。

国境線や州境も入ってますね。
オプションは次回くわしく見てみます。  


Posted by 和歌山サイエンスカフェインフィニティ at 21:16Comments(0)GMT

2018年07月12日

【GMT講座11】pscoast その3

エクササイズを見てみると、
gmt.confファイルのFORMAT_GEO_MAPを変えて
試してみましょう。
とあります。

gmt.confファイルはGMTの環境設定をする
ファイルで、これで、いろんなことを設定できるようです。

設定を変更するには gmtsetコマンドを使います。

では、
$ gmt gmtset FORMAT_GEO_MAP ddd:mm
と打ち込んでみます。

同じ範囲で pscoastを使ってみますが
違いを見やすくするために -B オプションを -Ba2.5として
2.5度づつ目盛りを打つようにします。

$ gmt pscoast -R270/290/0/20 -JM6i -P -Ba2.5 -Gchocolate > GMT_tut_3b.ps

を実行すると、

となります。

度分表示になります。
0.5度ずつなので、30’(分)刻みで表示されています。

度だけで、1度未満が小数点以下で表示されるようにするには
$ gmt gmtset FORMAT_GEO_MAP D
とします。

同じコマンド(出力ファイルだけ変えてます)
$ gmt pscoast -R270/290/0/20 -JM6i -P -Ba2.5 -Gchocolate > GMT_tut_3c.ps
を実行すると


違いが分かると思います。

gmt.confは奥が広そうなので、私も勉強しつつ、今後も取り上げるようにします。
  


Posted by 和歌山サイエンスカフェインフィニティ at 22:58Comments(0)GMT

2018年07月11日

【GMT講座10】pscoast その2

前回実行したコマンドは
$ gmt pscoast -R-90/-70/0/20 -JM6i -P -Ba -Gchocolate > GMT_tut_3.ps
です。

-Rオプションは、psbasemapと同じですね。
x軸、y軸の範囲を設定していますが、
ここでは緯度、経度になっています。
また、マイナスは西経、南緯を表します。
プラスの場合は東経、北緯です。
つまり、西経90度から西経70度と
北緯0度(赤道)から北緯20度の範囲を
指定しています。

-Jオプションは投影の仕方です。
Mでメルカトル図法を指定しています。
6iは6インチで、グラフの大きさを指定しています。

-Bオプションは境界の描き方です。
aで目盛りを描くと指定しています。

-Gは塗りつぶし方を指定していて、
chocolate色に指定しています。

-Pは用紙を縦置きにするためのオプションです。  


Posted by 和歌山サイエンスカフェインフィニティ at 21:48Comments(0)GMT

2018年06月26日

【GMT講座9】pscoast その1

今回から、地図を描いていきます。
pscoastというコマンドを使います。


では、さっそく、チュートリアルを参考に
$ gmt pscoast -R-90/-70/0/20 -JM6i -P -Ba -Gchocolate > GMT_tut_3.ps
と打ち込んでください。

もしかすると、エラーが発生する場合があります。

pscoast: GSHHG version 2.2.0 or newer is needed to use coastlines with GMT.
Get and install GSHHG from ftp://ftp.soest.hawaii.edu/gshhg/.
pscoast: Could not find file [GSHHG low resolution shorelines]
pscoast: No GSHHG databases available - must abort

海岸線のデータであるGSHHGというものがインストールされていないときに
発生します。

$ sudo apt install gmt-gshhg
と実行するとデータがインストールできます。
より高精度なものや、粗いものもあります。

その上で、再度、先程のpscoastコマンドを実行してください。


地球のあるところの海岸線が描かれます。
縦軸、横軸が緯度、経度を表していますので、
調べてみてください。

オプションは次回に説明していきます。


  


Posted by 和歌山サイエンスカフェインフィニティ at 22:16Comments(0)GMT

2018年06月24日

【GMT講座8】psbasemap その6

psbasemapも今回で終わりです。
チュートリアルの
http://gmt.soest.hawaii.edu/doc/5.4.3/pdf/GMT_Tutorial.pdf
エクササイズを試してみましょう。

今回はg3オプションを-Bオプションにつけてみます。

$ gmt psbasemap -R1/10000/1e20/1e25 -JX9il/6il -Bxa2f3g3+l"wevelength(m)" -Bya1pf3g3+3+l"Power(w)" -BWS > GMT_tut_2.ps

すると、グリッドが引かれます。

psbasemap8
  


Posted by 和歌山サイエンスカフェインフィニティ at 16:34Comments(0)GMT

2018年06月22日

【GMT講座7】psbasemap その5

gmt psbasemap -R1/10000/1e20/1e25 -JX9il/6il -Bxa2+l"Wavelength (m)"
-Bya1pf3+l"Power(W)" -BWS > GMT_tut_2.ps

-Rオプションでx軸とy軸のそれぞれの範囲を指定します。x軸は1から10,000まで、
y軸は10の20乗から10の25乗までと指定しています。1e20は指数表記でeの前後で
仮数部、指数部を表しています。

-Jオプションで、投影方法と縦横の大きさを指定します。前回と同様、大文字X
で直角座標を指定しています。横幅が9インチ、縦幅が6インチで、対数グラフの
場合は小文字のlを付けます。

-Bオプションで、軸の目盛等を指定します。
今回は小文字のxとyをオプションに付けることで
x軸とy軸を分けて設定しています。
また、小文字lをオプションに追加して軸毎のラベルを設定しています。

最後の -BWS も -Bオプションの一つですが、
WはWestで左側の縦軸、SはSouthで下側の横軸を表しています。
大文字にすると、軸と目盛などが表示されます。
小文字にすると、軸だけが表示されます。
指定しないと表示されません。
ここでは右側の縦軸E(e)と上側の横軸N(n)が指定されていないので
表示されない設定となっています。


  


Posted by 和歌山サイエンスカフェインフィニティ at 21:54Comments(0)GMT

2018年06月21日

【GMT講座6】psbasemap その4

チュートリアルを参考に、次にいきましょう。

今回も、psbasemapですね。
両対数のグラフになります。


$ gmt psbasemap -R1/10000/1e20/1e25 -JX9il/6il -Bxa2+l"Wavelength (m)" -Bya1pf3+l"Power(W)" -BWS > GMT_tut_2.ps

を打ち込み、実行すると、



横置きの紙に
下と左のみ軸が描かれて、それぞれが対数表示になっていますね。
またそれぞれの軸にラベルがつきました。

次回で、オプションの解説をします。  


Posted by 和歌山サイエンスカフェインフィニティ at 23:11Comments(0)GMT

2018年06月20日

【GMT講座5】psbasemap その3

前回のpsbasemapのオプションを少し変えてみて
どうなるか見てみましょう。

グラフの範囲をx軸について、−100から100まで
y軸を-60から60までにします。
-Rオプションは -R-100/100/-60/60にします。

目盛りを20毎に数値を入れて、小さい目盛りを10毎に入れ、
20毎にグリッドの線も入れてみます。
-Bオプションに f と gを加えます。また数値を入れてそれぞれ指定します。
-Ba20f10g20とします。

グラフの大きさは横に5インチ、縦に3インチにします。
-Jオプションは-JX5i/3iとします。

描画域の色を青色に
タイトルも「My Second Plot」に変更してみましょう。
-B+gblue+t"My Second Plot"とします。

紙は横置きのため -P オプションは抜きました。

コマンドは
$ gmt psbasemap -R-100/100/-60/60 -JX5i/3i -Ba20f10g20 -B+gblue+t"My Second Plot" > GMT_tut_1b.ps
になります。
出力ファイル名も変更しています。

出来上がる図は次のようになります。


  


Posted by 和歌山サイエンスカフェインフィニティ at 22:49Comments(0)GMT

2018年06月19日

【GMT講座4】psbasemap その2

前回では、
$ gmt psbasemap -R10/70/-3/8 -JX4i/3i -Ba -B+glightred+t"My first plot" -P
> GMT_tut_1.ps
を実行してみました。

それぞれのコマンドを詳しく見ていきましょう。

GMTにおいて、コマンドのことを調べるには、
man pagesを参照するのが一番良いです。
http://gmt.soest.hawaii.edu/doc/5.3.2/man_pages.html
ただ、英語なのでいろいろ困難を伴うところもあります。

ちなみに、psbasemapコマンドは
http://gmt.soest.hawaii.edu/doc/5.3.2/psbasemap.html
を見てください。


まず、冒頭に
psbasemapは
Postscript(ポストスクリプト)
の白地図を描くものです
と書かれています。

グラフの大きさの設定や、
横軸(x軸)・縦軸(y軸)について範囲や目盛の付け方、
グリッドを描くか、タイトルやラベル等の
設定をするために、まず呼び出すコマンドと
考えるといいと思います。

グラフの中身については、psbasemapで描いた白地図(軸など)
に上書きしながら付け加えていくイメージで大丈夫です。


記述の仕方は、
psbasemap -Jparameters -Rwest/east/south/north[/zmin/zmax][r] [ -B[p|s]
parameters ] [ -A[file] ] [ -Dinsert box ] [ -Fbox ] [ -K ] [ -Jz|
Zparameters ] [ -Lscalebar ] [ -O ] [ -P ] [ -U[stamp] ] [ -Trose ] [ -
Tmag_rose ] [ -V[level] ] [ -Xx_offset ] [ -Yy_offset ] [ -ccopies ] [ -
fflags ] [ -pflags ] [ -ttransp ]
となっています。

たくさんオプションがありますが、
-J -Rオプションは必須のようです。
[ ]で閉じられたオプションは必要に応じて設定します。

前回の例では、
-R10/70/-3/8
-JX4i/3i
-Ba
-B+glightred+t"My first plot"
-P
の5つオプションを使用しています。

-Rオプションは
http://gmt.soest.hawaii.edu/doc/5.3.2/gmt.html#r-full
に詳しい説明があります。
このオプションでx軸、y軸の範囲を設定します。
ここでは、x軸を10から70までの範囲に、
y軸を-3から8までの範囲に設定しています。


-Jオプションは
http://gmt.soest.hawaii.edu/doc/5.3.2/gmt.html#j-full
に詳しい説明があります。
-Jオプションは、地図やグラフの投影の仕方を設定します。
Jに続けて文字を加えることで設定します。
たとえば、-Jmとするとメルカトル図法で描く等の設定があります。
ここでは -JX4i/3i となっていますので、Xで直角座標系の投影法を
指定し、4i/3iでx軸、y軸のそれぞれの長さを指定しています。
iはインチを表すので、x軸が4インチ、y軸が3インチとなります。
注意すべきは小文字を用いた-Jxでは、違う設定になる点です。
小文字xを用いた場合、1単位当たりの長さを指定することになります。



-Bオプションは
軸等の設定をするオプションです。
まず-Baのaは「anotation」の頭文字から取られています。
目盛の数値を描くかを設定します。aに続けて数値を入れると
目盛の数値の間隔を設定できます。
また、-B+glightred+t"My first plot" は
+gに続いて描画領域の色を設定し、
+tに続く文字列でタイトルを設定します。

最後の-Pオプションは
紙を縦置きに指定するものです。
「Portrait」のPです。
指定しないと横置き(Landscape)になります。

では、もう一度出力を確認してみましょう
psbasemap2

  


Posted by 和歌山サイエンスカフェインフィニティ at 20:31Comments(0)GMT

2018年06月18日

【GMT講座3】psbasemap その1

では、さっそくGMTを始めましょう。

本家GMT wiki
http://gmt.soest.hawaii.edu/

チュートリアルを参考にすすめるのがいいと思います。
http://gmt.soest.hawaii.edu/doc/5.4.3/pdf/GMT_Tutorial.pdf

pdf版がダウンロードできます。

まずは、基本コマンドである
psbasemap
を試してみましょう。

GMTのバージョン5からは、
$ gmt psbasemap
のように
gmtというコマンドに続けて、入力します。

これだけ入力してリターンすると、

psbasemap1

となってエラーがたくさん出てきちゃいます。

GMTでは、コマンドに引き続いて、いくつかのオプションをつけなければいけません。
-J や -R と言ったものです。
他にも、出力される画像(ポストスクリプト形式のもの)も指定しないといけません。

tutorialでは、
$ gmt psbasemap -R10/70/-3/8 -JX4i/3i -Ba -B+glightred+t"My first plot" -P > GMT_tut_1.ps
となっているので、
入力すると GMT_tut_1.ps
というファイルができあがります。

psbasemap2

こんな感じ。
次回はこの説明で。

  


Posted by 和歌山サイエンスカフェインフィニティ at 21:20Comments(0)GMT

2018年06月17日

【GMT講座2】GMTをインストールする

第2回目は、GMTをインストールします。

ubuntu18.04では、ターミナル(端末)を立ち上げ
$ sudo apt install gmt
と打ち込み、続いてパスワードを入力すれば、
インストールされる。

海岸線のデータもインストールしたい場合は、
$ sudo apt install gmt-gshhg gmt-gshhg-high gmt-gshhg-low
を実行すればよい。

随分、簡単になりました。

$ gmt
と打ち込んで、エラーがでなければOKです。

$ gmt --version
と打ちこむと、バージョン情報がでます。
  


Posted by 和歌山サイエンスカフェインフィニティ at 12:07Comments(0)GMT

2018年06月15日

【GMT講座1】GMTとは

GMTとはGeneral Mapping Toolsの略で、
ハワイ大学で開発されている。
オープンソースの約80のコマンドラインで
扱う命令群である。
高度な地図やグラフが作れ、無料で利用が可能である。
GSHHGの海岸線のデータやDCWの国境線のデータなども
簡便に用いることができる。

GMTは現在、世界各地のボランティアやNationalScienceFoundationの
サポートのもと、Paul Wessel氏、Walter H. F. Smith氏、
Remko Scharroo氏、Joaquim Luis氏とFlorian Wobbe氏により、
開発・メンテナンスが行われている。

GMTに関する情報は、ハワイ大学のサイト
http://gmt.soest.hawaii.edu/
で得ることができる。

WindowsやMacでも利用可能であるが、Linuxでの利用が一番利用しやすく
情報も多い。
当サイトではUbuntu(18.04LTS)上でGMTのバージョン5.4.2を
用いる形で当面は進めていく。  


Posted by 和歌山サイエンスカフェインフィニティ at 19:12Comments(0)GMT

2017年06月30日

第1回 arduino勉強会

6月25日に、和歌山市内の北ぶらくり丁会館にて
arduino勉強会を開催しました。

arduino(アルドゥイーノ)とは、オープンソースハードウエアの
代表的なデバイスで、
コンピュータや電子機器の理解、プログラミング言語の
習得などに最適な教材です。

arduinoはいくつかの種類がありますが、
もっとも一般的なarduino unoを使って、
しばらく基礎的な勉強を重ねていく予定です。




こんな感じの、手のひらサイズのかわいいやつです。
日本なら、3千円くらいで手に入ります。



これをパソコンに繋いで、C言語風のプログラムを
arduinoに流し込むと、いろんな動作が可能になります。
arduinoには、センサーやLED,モーターなどの電子部品を
接続すれば、センサーからの情報を取得したり、
LEDに出力したり、モーターを動かすことができます。

これらを組み合せば、いろんなことが可能になるのです!

今回は初めての勉強会だったので、
arduinoの基本事項を勉強して、
・arduino上のLEDを光らせる(Lチカ)
・ブレッドボードを接続して、ブレッドボードに
 付けたLEDを光らせる
・タクトスイッチを使って、パソコン画面に
 タクトスイッチのオン・オフ状態を表示する
・タクトスイッチを使って、LEDをオン・オフさせる
作業をしてみました。






今後も、定期的に開催予定です。
(次回は7月16日(日曜日)を予定)

電子工作に興味がある方や、
お子さんの夏休みの自由研究に使えそうって思われる方、
プログラミングの練習をしてみたい方
など、
幅広く楽しめるツールです。


ご興味ある方は、ぜひご連絡ください。
wsc.infinity@gmail.com
  


Posted by 和歌山サイエンスカフェインフィニティ at 21:07Comments(0)arduino勉強会

2017年06月25日

わかやま路上観察学 無用扉

超芸術トマソンというものがありますが、
これらは和歌山でも散見されます。

都市の間で生まれる超芸術ですが、
見ようと思って探さないと、なかなか
見つけることができない面白いオブジェクトです。

トマソンには、いくつかのタイプに分類されますが、
今回は「無用扉・無用窓」と「原爆タイプ」を取り上げます。

1.無用扉と無用窓
建築物は基本的に壁と柱から構成されますが、
人間やモノの出入りのため、また、室内環境を
快適にするため、扉や窓といった開口部が設けられます。

その開口部を、なにかの理由で閉じたもの。
無用となった開口部がトマソンとして鑑賞できます。



2.原爆タイプ
都市には密接した建物が多数あります。
その隣の建物が、何かの理由で忽然となくなってしまうと、
隣の建物に現れるトマソンの一種。

在りし日の姿がしのばれれます。


  


Posted by 和歌山サイエンスカフェインフィニティ at 11:52Comments(0)わかやま路上観察学

2017年06月20日

第1弾 とある駅のベンチ

新シリーズの
『わかやま路上観察学』

和歌山の街で見つけた、美しいもの
面白いものをアップしていくシリーズ。

超芸術トマソンや、なにげないものを
見つけていきます。

評価は読者におまかせ。

第1弾
「とある駅で見つけたベンチ」

電線を架線するコンクリート柱の
ための気遣いと手仕事の丁寧さを感じる






  


Posted by 和歌山サイエンスカフェインフィニティ at 21:44Comments(0)わかやま路上観察学

2017年02月26日

2017初頭のウルトラブック比較

2017年初頭のモバイルノート(ウルトラブック)な
パソコンを比較検討中。
それなりのCPUパワーと、16GB、512GB(できれば1TB)のSSD
で、1.3キロくらいより軽くて、コンパクトな電源アダプターなのが
必要条件だとする。
お金を出せれば候補も増えるので、20万円以下くらいが候補。

1.MacBookPro 13インチ(タッチバーなし)
  メモリを16GB、SSDを512GBにカスタマイズして、
  179,800円(税別、学割適用)
  http://www.apple.com/jp_edu_1460/shop/buy-mac/macbook-pro?product=MLL42J/A&step=config#
  いまなら、8,500円分のApple Storeギフトカードのプレゼント付き。

  良い点
  ・英語キーボードが選べる。
  ・Retinaディスプレイになる。
  ・タッチバーなしモデルは、コスパは良好
  悪い点
  ・CPUがskylake。秋くらいに Kabylakeになるとも言われている。)
・type-C に対応した、SDカードリーダとHDMIソケットのアダプターを買う必要が
   自分の場合ある。
  

2.東芝 dynabook V82/B (東芝ダイレクトだと VZ72/B)
Core i7 7500U で8GBの512GBのSSD(Serial ATA)
214、000円(税別) (東芝ダイレクトの会員になると割引あり)
https://dynabook.com/2in1-mobile-notebook-tablet/v-series/v82b-2016-winter-model-onyx-metallic-12-5-inch-high-end-convertible-pv82bmp-nja.html

 良い点
 ・国内メーカーという点
 ・いま所有してるR632は、非常に良いPCだった
 ・ノングレアな画面で、タッチ対応
 ・CPUはKabylake
 ・HDMIもVGA接続も可能なアダプター付き(電源供給も可能)

 悪い点
 ・英語キーボードなし
 ・メモリ8GB設定なし
 ・SSDに1TB設定なし。
 ・ちょっとベゼル幅が大きめ

3.HP spectre x360(パフォーマンスモデル)
  CPUはCore i7 7500U、メモリは16GBで、1TBのSSD(NVMe接続)
  161,800円(税別、学割適用、電話申し込みの場合)
HPの学割サイト http://jp.ext.hp.com/campaign/personal/notebooks/2017student_pc_collection/
http://jp.ext.hp.com/directplus/personal/?jumpid=ps_74m86jxecn

 良い点
 ・CPUはKabylake
 ・16GBも1TBも満たす。しかも、NVMe接続
 ・コストパフォーマンつ高め
 ・2月28日までに申し込むと、Office Businessが無料

 悪い点
 ・やや重め(1.3kg)
 ・英語キーボードなし
 ・電源アダプターのコードが太め
 ・HPにやや悪い印象


他にもDELL(XPS13)とかLenovo(Carbon)なども候補にあがるかも。
隠し玉はXiaomiですが、一応、法的な問題もあるので。

spectrex360とXPS13との比較動画


spectre360とMacBookPro13との比較動画

  


Posted by 和歌山サイエンスカフェインフィニティ at 19:08Comments(0)スタッフ日記

2017年02月21日

Javaを使った有限体積法による数値計算(2)

第2回目では、さっそくソースをアップする。

今回は、有限体積法の教科書としては、詳しい解説と、例題の丁寧な
回答が掲載されている
数値流体力学 第2版
の例題6.2をコード化した。

定常の移流方程式を解いている。
この例題を解いてみたのは、これから解きたい浅水流方程式に
近い形であること。つまり、速度場と密度場が未知であり、
これらを求めるためにSIMPLEアルゴリズムを適用するものであるためである。
また、断面形状も一定でなく、変化しているところも
今後の応用性が高く、参考になると思う。

Javaについては、1か月程度の知識しかなく、
静的メソッドのみの使用と、mainメソッドにほとんどの必要な処理を書いた。
初めて作ったJavaコードなので、至らないところが多い。

勉強しながら、現在はクラス等を用いたコードに修正を加えている。
また、多次元配列の初期化については、
Arrays#fillメソッドの多次元版の一例
を利用させていただいた。

連立方程式の解法には、トーマス法を用いている。
格子には等間隔のスタッガード格子を採用している。
また、例題の解答を直接コーディング化ではなく、最後の方に触れられている
不足緩和も係数計算の際に、組み込んでいる。

格子の数や不足緩和の係数を、いろいろ試してみたり、
質量フラックスを計算したりしてみて、有限体積法の面白さを
感じることができると思う。

(注)コードの利用による損害等は一切責任を負いません。

/*
* 有限体積法の計算
* 数値流体力学(第2版)
*
* 2017/2/11
* KZ
*
*/

package FVM;

import java.util.Arrays;

public class FvmTest007 {

public static void main(String args[]){

int nn = 5; // 圧力ポイントの数
int it_max = 200; // 最大回数
int count = 0; // 計算回数の記憶用
double r_limit = 10e-5; // 残差の制限


double m_0 = 1.0; // 初期仮定の質量フラックス
double p_0 = 10.0; // 上流端の圧力
double rho = 1.0; // 密度
double alpha = 0.7; // 不足緩和の係数 (nn,alpha)=(10,0.7),(20,0.5)
double residual = 0.0; // 残差

//double[] A_p = {0.5 ,0.4 ,0.3 ,0.2 ,0.1}; //圧力ポイントの面積
double[] A_p = new double[nn];
for(int i=0; i < nn; i++){
A_p[i] = 0.5 - (0.5-0.1)/(nn-1) * i;
}
//double[] A_u = {0.45,0.35,0.25,0.15};
double[] A_u = new double[nn-1];
for(int i=0; i < nn-1; i++){
A_u[i] = ( A_p[i] + A_p[i+1] ) / 2.0;
}

double[] x_u = new double[nn-1]; // 速度場
double[] u_star = new double[nn-1]; // 補正後の速度
double[] u_old = new double[nn-1]; // 速度の古いもの保存(残差計算用)
//double[] x_p = new double[nn]; // 圧力場
//double[] x_p = {10.0, 7.5, 5.0, 2.5, 0.0}; // 圧力場と仮定初期値
double[] x_p = new double[nn];
for(int i = 0; i < nn; i++){
x_p[i] = 10.0 - (10.0 - 0.0)/(nn-1) * i;
}
double[] p_prime = new double[nn]; // 圧力の補正値
double[] p_star = new double[nn]; // 補正後の圧力
double[] d = new double[nn-1]; // 圧力補正のパラメータ

double u_in = 2.0; // 上流端の流速(境界条件)

//方程式の係数
double F_w = 0.0;
double F_e = 0.0;
double a_W = 0.0;
double a_E = 0.0;
double a_P = 0.0;
double S_u = 0.0;

double b_prime = 0.0;

//連立方程式の係数行列
double[][] matrixA = new double[nn-1][nn-1];
double[] matrixB = new double[nn-1];


double[][] matrixC = new double[nn][nn];
double[] matrixD = new double[nn];


//配列の初期値
Arrays.fill(x_u,0.0);
Arrays.fill(u_star, 0.0);
Arrays.fill(p_prime, 0.0);
Arrays.fill(p_star, 0.0);
Arrays.fill(d, 0.0);
ArraysfillEx.fill(matrixA, 0.0); //多次元配列の初期値
Arrays.fill(matrixB, 0.0);
ArraysfillEx.fill(matrixC, 0.0); //多次元配列の初期値
Arrays.fill(matrixD, 0.0);



//速度場の初期仮定の設定
for(int i = 0; i < nn-1 ; i++){
x_u[i] = m_0 / ( rho * A_u[i] );
}

System.out.println("x_u : "+Arrays.toString(x_u)); // 行列の標準出力
System.out.println("x_p : "+Arrays.toString(x_p)); // 行列の標準出力
System.out.println("A_u : "+Arrays.toString(A_u));
System.out.println("A_p : "+Arrays.toString(A_p));



//計算の開始
for(int it = 0; it < it_max ; it++){

System.out.println("Count : " + it);

// 速度の更新(残差計算用)
for(int i = 0; i < nn-1 ; i++){
u_old[i] = x_u[i];
}
residual = 0.0; // 残差の初期化


//運動量保存則 始まり
//速度の格子点//
//格子点 1(0)
u_in = x_u[0] * (A_u[0] / A_p[0]); System.out.println("u_in : " + u_in);
F_w = rho * u_in * A_p[0]; //System.out.println(F_w);
F_e = rho * (x_u[0] + x_u[1])/2 * A_p[1]; //System.out.println(F_e);

a_W = 0.0;
a_E = 0.0;
a_P = F_e + F_w * 0.5 * Math.pow((A_u[0]/A_p[0]),2.);


S_u = (x_p[0]-x_p[1])*A_u[0] + F_w*(A_u[0]/A_p[0])*x_u[0] + (1-alpha)*a_P/alpha * x_u[0];
a_P = a_P/alpha; //不足緩和

matrixA[0][0] = a_P;
matrixA[0][1] = -1.0 * a_E;

matrixB[0] = S_u;

d[0] = A_u[0]/ a_P * alpha;//不足緩和


//格子点 2,3(1 〜 nn-3)
for(int i=1 ; i < nn-2;i++){
F_w = rho * (x_u[i-1] + x_u[i])/2 * A_p[i];
F_e = rho * (x_u[i] + x_u[i+1])/2 * A_p[i+1];
a_W = Max(F_w,0.0);
a_E = Max(0.0, -1.0*F_e);
a_P = a_W + a_E + (F_e - F_w);

S_u = (x_p[i]-x_p[i+1])*A_u[i] + (1-alpha)*a_P/alpha * x_u[i];//不足緩和
a_P = a_P/alpha; //不足緩和

matrixA[i][i-1] = -1.0 * a_W;
matrixA[i][i] = a_P;
matrixA[i][i+1] = -1.0 * a_E;

matrixB[i] = S_u;

d[i] = A_u[i]/ a_P * alpha;
}


//格子点 4(nn-2)
F_w = rho * (x_u[nn-3] + x_u[nn-2])/2 * A_p[nn-2]; //System.out.println(F_w);
//F_e = rho * (x_u[nn-2] + x_u[1])/2 * A_p[1]; //System.out.println(F_e);
if(it == 0) {
F_e = m_0; // 最初の計算だけ 要検討!!!
}else{
F_e = rho * x_u[nn-2] * A_u[nn-2];
}

a_W = F_w;
a_E = 0.0;
a_P = a_W + a_E + (F_e - F_w);

S_u = (x_p[nn-2]- x_p[nn-1])*A_u[nn-2] + (1-alpha)*a_P/alpha * x_u[nn-2];
a_P = a_P/alpha; //不足緩和

matrixA[nn-2][nn-3] = -1.0 * a_W;
matrixA[nn-2][nn-2] = a_P;

matrixB[nn-2] = S_u;

d[nn-2] = A_u[nn-2]/ a_P * alpha;



//System.out.println("matA: "+Arrays.deepToString(matrixA));
//System.out.println("matB: "+Arrays.toString(matrixB));
//System.out.println(" d : "+Arrays.toString(d));

//連立方程式を解く

u_star = thomas(nn-1,matrixA,matrixB);

System.out.println("u_star: " + Arrays.toString(u_star));

//運動量保存則 終わり

//圧力の補正式 始まり
//圧力の格子点 1,A (0)
matrixC[0][0] = -1.0;

//圧力の格子点 2から4,A (1 - nn-2)
for(int i = 1; i < nn-1 ; i++){
a_W = rho * d[i-1] * A_u[i-1];
a_E = rho * d[i] * A_u[i];
F_w = rho * u_star[i-1] * A_u[i-1];
F_e = rho * u_star[i] * A_u[i];
a_P = a_W + a_E;
b_prime = F_w - F_e;

matrixC[i][i-1] = -1.0 * a_W;
matrixC[i][i] = a_P;
matrixC[i][i+1] = -1.0 * a_E;

matrixD[i] = b_prime;
}

//圧力の格子点 5,E (nn-1)
matrixC[nn-1][nn-1] = -1.0;


//System.out.println("matC: "+Arrays.deepToString(matrixC));
//System.out.println("matD: "+Arrays.toString(matrixD));

//連立方程式を解く
p_prime = thomas(nn,matrixC,matrixD);

System.out.println("p_prime: " + Arrays.toString(p_prime));

//圧力の補正式 終わり

//圧力の補正 始り
for(int i=0;i < nn;i++){
x_p[i] = x_p[i] + alpha * p_prime[i];
}
System.out.println("calculated x_p : " + Arrays.toString(x_p));
//圧力の補正 終わり

//速度の補正 始まり
for(int i=0;i < nn-1;i++){
x_u[i] = u_star[i] + d[i] * (p_prime[i] - p_prime[i+1]);
}
System.out.println("calculated x_u : " + Arrays.toString(x_u));
//速度の補正 終わり
//格子点Aの圧力の補正
x_p[0] = p_0 - 0.5 * rho * Math.pow(x_u[0], 2.0) * Math.pow(A_u[0]/A_p[0],2.0);
System.out.println("calculated2 x_p : " + Arrays.toString(x_p));
// 終わり

// 残差の計算
for(int i = 0; i < nn-1 ; i++){
residual += Math.pow(u_old[i]-x_u[i],2.0);
}
System.out.println("residual : " + residual);

if(residual < r_limit) {
count = it;
break;
}

}

System.out.println("¥n Count: " + count);
System.out.println("x_u : " + Arrays.toString(x_u));
System.out.println("x_p : " + Arrays.toString(x_p));
System.out.println("A_u : "+Arrays.toString(A_u));
System.out.println("A_p : "+Arrays.toString(A_p));



}



//以下、メソッド

//大きい数字を返す
public static double Max(double a,double b){
if(a>b){
return a;
} else {
return b;
}
}


//トーマス法
public static double[] thomas(int n,double[][] A,double[] B){

double[] x = new double[n];

double[] g = new double[n];
double[] s = new double[n];

//配列の初期値
Arrays.fill(x,0.0);
Arrays.fill(g, 0.0);
Arrays.fill(s, 0.0);

g[0] = A[0][0];
s[0] = B[0];

for(int i = 1 ; i < n ; i++){
s[i] = B[i] - A[i][i-1]*s[i-1]/g[i-1];
g[i] = A[i][i] - A[i][i-1]*A[i-1][i]/g[i-1];
}

x[n-1] = s[n-1]/g[n-1];

for(int i = n-2 ; i >= 0; i--){
x[i] = (s[i]-A[i][i+1]*x[i+1])/g[i];
}

return x;
}
}

  


Posted by 和歌山サイエンスカフェインフィニティ at 23:53Comments(0)気になる水理学

2017年02月18日

Latexをwindowsにインストール

windows7ノートに
Latexをインストールした。

Texインストーラ3を利用
http://www.math.sci.hokudai.ac.jp/~abenori/soft/abtexinst.html

参考にしたサイト
http://did2memo.net/2016/04/24/easy-latex-install-windows-10-2016-04/  

Posted by 和歌山サイエンスカフェインフィニティ at 21:28Comments(0)スタッフ日記

2017年02月16日

Javaを使った有限体積法による数値計算(1)

ここ数ヶ月、流体の数値計算の勉強を続けている。
2次元の河床変動プログラムを作るのが、当面の最終目標。

河床変動の計算の基礎式は、
1:流体の質量保存則(連続の式)
2;流体の運動量保存則(ナビエストークス方程式)
3;流砂の連続式
4;流砂量式
である。これら4つの方程式を解くことで
流れの場を一般的には求めることができる。

これらの式は複雑に絡み合い、特に2の運動量保存則は、
非線形の方程式(移流方程式)であり、一般に解析的に解くことが
できない。
そこで、数値計算が多用される。

数値計算においては、上記の方程式を離散化して、
解く方法がとられている。

その離散化の方法の違いにより、大きく
1: 格子法
2: 粒子法
に区別される
(他に個別要素法があるが、どう分類されるのかはまだ、勉強不足)

最近は2の粒子法も注目を集めているが、
計算範囲が大きい場合や、時間スケールを長くするには、
計算パワーが必要なことから、まだまだ一般的ではない。
粒子法については
粒子法入門 流体シミュレーションの基礎から並列計算と可視化まで C/C++ソースコード付
が詳しくて分かりやすい。

このため、一般には1の格子法が採用される。
格子法は、さらに、
A:差分法
B:有限体積法
C;有限要素法
などに分類される。

差分法は、微分方程式をテイラー展開を利用して離散化する。
また、有限要素法は、重み付きの関数を用いて離散化する。
(有限要素法も勉強不足)
一方で、有限体積法は、コントロールボリュームと呼ばれる
有限な体積の要素に積分することで、その収支を用いて離散化する。
有限体積法では、コントロールボリュームの質量収支等を考慮するので、
質量が保存されやすいという特があり、この点が他の計算方法との違いと
言われている。

連続体としての質量の保存を重視する流体力学の数値計算では、
有限体積法が人気を集めているようである。
他の工学分野(構造力学等)では、有限要素法が一般的なようである。

以上のような理由から、有限体積法を用いて、
数値計算(数値シミュレーション)を進めていくこととした。


なお、使用言語は、
Javaを使っていくことにする(途中で変更するかもしれない)
もともと数値計算でよく用いられているFortranも使っていたが、
Androidアプリ開発のため、Javaを勉強し始めたことと、
オブジェクト指向を理解するためにも、Javaをプログラミング言語として
学んでいくことは、今後に役立つと考えたためである。


有限体積法の勉強に用いた参考文献を以下に示す。








  


Posted by 和歌山サイエンスカフェインフィニティ at 21:07Comments(0)スタッフ日記気になる水理学

2016年12月12日

GMT5.3.1で陰影図や鳥瞰図を作る

GMT5.3.1を使い、国土地理院の基盤地図情報サービスの
10mメッシュのDEMから、地形の陰影図や鳥瞰図を作る。

まずはデータのダウンロード
http://www.gsi.go.jp/kiban/
(今回は紀伊半島の一部のデータを利用)

登録してダウンロードし解凍してできる
FG-GML-5035-37-dem10b-20161001.xml
を利用します。

GMTで扱えるgrd形式に変換する必要がある。
grd形式に変換するスクリプトを
産業技術研究所の野田さんが公開されているので
利用させていただく。ありがとうございます。
https://staff.aist.go.jp/a.noda/programs/jpgis2grd/jpgis2grd.html

このスクリプトを動かすとgrd形式に変換できるが
一部環境に合わせて変更する。
(commandを「gmt xyz2grd」とか出来上がるファイル名などの変数を)

$ perl jpgis2grd.pl
を実行
しかし、エラーで止まる。
途中のどこかに3列目の数値がないところがあるようだ。
(今回は 843751行目)

このスクリプトのglobal optionに
my $outfile2 = "out2.xyz"
を追記。

grdに変換するあたりを

print STDOUT "Making a grd file.\n";
# my $command="gmt xyz2grd $outfile -R$xmin/$xmax/$ymin/$ymax -I$xint/$yint -G$grdfile -N-9999 -Dlongitude/latitude/height/1/0/$opts{title}/$opts{remark} ";
my $command_sed="sed -e '843751d' $outfile > $outfile2";
my $command="gmt xyz2grd $outfile2 -R$xmin/$xmax/$ymin/$ymax -I$xint/$yint -G$grdfile -di-9999 -Dlongitude/latitude/height/1/0/$opts{title}/$opts{remark} ";
if ($opts{verbose}) { print STDERR "Running \"$command\"\n"; };
system("$command_sed");
system("$command");
system("gmt grdinfo $grdfile");

と一部修正。
(perlに詳しくないので、あんまりいいやり方ではなさそうですが)

再度、
$ perl jpgis2grd.pl -v -n
を実行
(実行状況や中間ファイルを残すのに -v -nオプションをつけた)

すると、out.grdファイルができるので
とりあえず、これを
FG-GML-5035-37-dem10b-20161001.grd
と名前を変更。

次はgmtで陰影図を描く
以下のスクリプトで
#!/bin/bash
#
# DEM10mメッシュのgrdファイルから陰影図を作る。
# 2016/12/12
#
# perl jpgis2grd.pl で出力される範囲をコピー
# xmin=135.87500000, xmax=136.00000000, ymin=33.58333333, ymax=33.66666667

grd_file=FG-GML-5035-37-dem10b-20161001.grd
shade_file=5035-37.shade
eps_file=5035-37_inei.eps
cpt_file=color.cpt
ctb=-1000/1000/50 #color table min/max/interval
region=135.875/136.00/33.5833333/33.6666667

gmt makecpt -Cglobe -T$ctb -Z > $cpt_file

gmt gmtset FORMAT_GEO_MAP D

gmt grdgradient $grd_file -A90 -G$shade_file -Ne0.6 -V
gmt grdimage $grd_file -R$region -JM15 -C$cpt_file -I$shade_file -X4 -Y4 -K > $eps_file
gmt psscale -D7.5/-1/12/0.3h -C$cpt_file -Ba200f100:"altitude(m)": -I -V -K -O >> $eps_file
gmt psbasemap -R$region -JM15 -Ba0.05f0.025g0.01 -O -V >> $eps_file

こんな図ができる



次は鳥瞰図
#!/bin/bash
#
# DEM10mメッシュのgrdファイルから陰影図を作る。
# 2016/12/12
#
# perl jpgis2grd.pl で出力される範囲をコピー
# xmin=135.87500000, xmax=136.00000000, ymin=33.58333333, ymax=33.66666667

grd_file=FG-GML-5035-37-dem10b-20161001.grd
shade_file=5035-37.shade
d3_file=5035-37_3d.eps
cpt_file=color.cpt
ctb=-1000/1000/50 #color table min/max/interval
region=135.875/136.00/33.5833333/33.6666667
boundary_cap=wSnEZ
gmt makecpt -Cglobe -T$ctb -Z > $cpt_file
gmt gmtset FORMAT_GEO_MAP D
gmt grdgradient $grd_file -A90 -G$shade_file -Ne0.6 -V
gmt grdview $grd_file -Jm130 -Jz0.002 -C$cpt_file -I$shade_file -p150/25 -Ba0.05f0.025g0.01 -Bza400 -B$boundary_cap -Qi300 -N0+ggray -V > $d3_file

を実行すると、鳥瞰図ができる。



角度とかは適宜設定する。

  
タグ :gmt5.3.1grdmac


Posted by 和歌山サイエンスカフェインフィニティ at 15:55Comments(0)スタッフ日記