prasinos' work memo

スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

netcdf 3.6 (2GB 超ファイルに対応)

NetCDF 3.6 というのがしばらく前に出ていたらしいことにいまごろ気がつきました。
http://my.unidata.ucar.edu/content/software/netcdf/index.html
まだ何にも試してませんが、2GB 超ファイルに対応したそうです。

とりあえず速報。
スポンサーサイト

firefox の Java プラグイン

何かイチヂルシクどうでもいい記事ばかりとなりつつある昨今なのであるが。

ラップトップに Vine Linux をインストールして 3 ヶ月くらいになるらしい。2日くらい動かしたあとで firefox を動かしているとフリーズすることがあるが、その他は順調である。えーと、サウンドとモデムが全滅という、吾輩いつもの泣き所はあるのだが、まあそれは半分あきらめている。

それでもって、今日は firefox で java アプレットが見られないことに気がついた。みてみると、何故か jdk とかいう rpm パッケージが入っており、おそらく過去の俺様がダウンロードしてインストールした気配である。じゃあなぜ動かないのヨアンタというと、firefox の plugins ディレクトリに symlink を張らなきゃいけないんだという。

# ln -s /usr/java/jdk1.5.0_01/jre/plugin/i386/ns7/libjavaplugin_oji.so /opt/mozilla-firefox/plugins/

(1) コピーしちゃダメなんだと。(2) jdk は ns7/ 以外に ns7-gcc29/ とかいうのも用意しているが、こちらもダメ。いずれも firefox が黙ってクラッシュするにつき、GUI しているとかなり焦る。

こういうドキュメントは、よくがんばってくれているのはわかるのだけど、あんまり椀ストップになっていないのでまいごになってしまう。もうちょっとなんとかならないのかと思うですな。無料でダウンsておきながら勝手なこといsてごめんなさい。

stdio を Fortranで作ろう!

stdio みたいな任意フォーマットの入出力をやるライブラリが Fortran にない
のが気に入らない。やっぱり作るしかないよね。そういうわけで考察。

・ダイレクトアクセスバイナリファイルで実現するしかない。 (HP-UX 以外のほとんどの処理系で、最終レコードの中に EOF があっても 読めてしまう)

・レコード長は、それなりに長くないと性能が悪い。でも、中途半端な レコードで EOF が来る書き込みをするためには、レコード長を短くしたほう
がいい。時に応じて切替えるべき。

CF Metadata ML での GIS 対応のドタバタ

このところ CF Metadata メーリングリストで喧嘩つーか、あまり意見がまとまりにくい話をしているようです。今まで斜め読みをしていたのですが、どうやら他人事ではないようなので、読み返してみたら結構面白いのです。

そもそも CF Metadata ML とは何者かというと、NetCDF ファイル形式でデータを保存する際にメタデータを記述する方法 (CF 規約; Climate and Forecast Conventions) を標準化(=とりまとめる)ことを目的とした集まりです。メタデータとはデータ製造元だとか単位だとかその他諸々ですね。そういうものは NetCDF では名前と値の対である属性として格納してね、といっているだけで、属性の名前、型(いろいろあります)、意味などは決めていません。XML でタグが自由に作れるのと同じですね。で、NetCDF としては利用者がタグセットじゃなかった規約を作れ、あとは知らん、となっているわけですから、NetCDF の誕生以後いろいろの規約ができたわけですけれども、所詮気象業界なんて狭いものですからバラバラに規約を作る労力を分散してもはかどらないし、互換性の問題もあるというところで、現在は CF 規約に統合されつつあるという状況であるわけです。

それで、GIS とのやりとりというのはまだ発展中の分野で、なにせ GIS では cm 単位での位置同定を行うためにカリカリやってきたのに、気象界では経緯度がわかればそれ以上気にしてこなかった、というギャップがあるので、物の見かたに差があるわけですね。

どうやらログによると話の発端は UTM ゾーン番号を書き込みたいというリクエストが出たことだったようです(参照)。

UTM とはユニバーサル横メルカトル図法のことで、横メルカトルとは軸が赤道面にある円筒を使ったメルカトル図法、それで UTM は単なるTransverse Mercator (ちなみに英語ではマーケイタみたいに発音するので戸惑います) に加えて中央経線を 6 度の半奇数倍に制約し、投影面上の x 座標値 (東距 easting といいます) に 500km 下駄をはかせて負にならないようにするもので、結局のところ中央経線を指示するゾーン番号、x座標、y座標を指定すると地球上の位置が定まるという代物です。そういえば、国土地理院の地形図も UTM 図法ですね。

それで、提案のココロはおそらく、x y 座標にメートル単位を使って記述する機構はすでにあるので、あとはゾーン番号だけ格納すれば完璧じゃんと思ったわけですね。そんな安易な提案では他の投影法を使う人が困るでしょう、とおもうのですがおいておいて、


「UTM ゾーン番号を保存する方法がないから、名前を割り当てよう!」
「もともと CF 規約では FGDC (米連邦地理情報委員会) の決めたメタデータ規
格にしたがって横メルカトル図法のパラメタを列挙して保存する方法が規定され
ているので、それを使いなさい」
「UTM と横メルカトル図法はほんとに同じ図法? 」
「そうです」


という会話になったわけですね。まあ、既に用意されているもので片がつくのにマニュアルの読みかたがよくわからない、という普通の技術系MLの質疑応答です。

ところが話はそれだけに留まらなかった。「xy 座標から地球上の位置を特定するためのパラメタって何だ」と考えると、実は UTM 投影法というだけでは決まらなくて、投影面を巻き付ける楕円体と、三角点網を楕円体上の位置に結びつける測地系を指定しないといけないわけですね。

そうすると、パラメタは結局


  • 投影法パラメタ ...

    • 投影法名称
    • 投影法依存事項 ...

      • longitude_of_central_meridian 中央経線
      • latitude_of_projection_origin 原点緯度
      • scale_factor_at_central_meridian 中央経線上の縮尺
      • false_easting 東距のゲタ
      • false_northing 北距のゲタ


  • 楕円体パラメタ ...

    • 赤道半径
    • 離心率

  • 測地系名称


てな感じになってくる。いやあ、多いですねえ、と思うのは人情なわけで、

「上記全部をひっくるめて1つの番号で指定できる、EPSG (欧州石油探査グループ ... って何だ?) (参照) が割り当てている番号を使おう。MS Access 持ってるひとは誰でもオンラインで番号から諸元が検索できるよ」

なんて意見がでてくるわけですけど、


「気象屋で GIS も Access も使う人口はほとんどいないじゃん」
「外部団体に依存するのはいかがなものか」
「登録されてない組合せは使えないじゃん」


と総スカン。中でちょっと目に止まったのが

「どこか遠くの組織に表を管理させるのって GRIB の思想だよね」

という言われかたでした。ともあれ EPSG 劣勢と思われたのも束の間、 Russ Rew 御大が登場して


「FGDC Metadata Std は古い。EPSG がデファクトだって ESRI 言ってるよ」


と言ったところで世の中大逆転の予感。

昔からの思想対立として、複雑なものを整数や浮動小数点数のような簡単なパラメタに分解して、自分で持つのがよい、というのが self-descriptivity とか self-containment を重んじる一つの立場なら、逆にそういうパラメタの組合せには実際にはたいした自由度がないので、表を参照させるのが情報の正しい圧縮になる、という表参照主義があるわけだけれど、ふつう、表がでかくなったら表参照主義はもうダメ、ということになるわけです。

ところが今時はネットワークがあたりまえになったので、膨大な表をあえて管理して、オンラインサービスをするというアプローチがあり得るようになってしまった、ということを EPSG の勝利は示している、のかなあ。いまいちよくわからん。

生卵とサルモネラ菌

欧米では生卵を食べてはいけないという。サルモネラ菌があるからだ。しかし、日本より特にあぶないのだろうか。どうも必ずしもそうでもないみたいだ。

とりあえず google してみると、さすが米国政府、CDCの一般向け解説サイトがトップにヒットする。それによると、汚染は米国北東部に多く、そのへんでは大体 10,000 個に 1 個の卵が(殻ではなく)内部でサルモネラ菌を持っているという。その他の地域ではそれほど多くない、というけれど、全米平均では 20,000 個に1個であるらしい(アメリカ卵協会 Egg Handling and Care Guideによる。このガイド、北東部では気をつけろと書いていないが、まあ業界団体だから仕方がないか)ので、桁で違うわけではないだろう。CDC は調理後2時間以内に食えとか常識的なガイドラインだが、卵協会はヴァセリンを塗ったら冷蔵しなくてもいいという迷信を信用しないように、とかいった怪しげなことも混ざっているいるので面白い。なお、卵の内部よりは殻にサルモネラ菌が多いのだそうで、そこから「割れた卵は買ってはいけない」というガイドラインになるのだろう。

さて、一方で日本においてはどうかというと、厚生労働省が発表した三重県の調査によると 24,000 個中 6 個が汚染されていた... というと、この調査が卵内部の汚染についてならば、 4000 個に1個だからアメリカ北東部より多いことになる(実のところこれ以上の資料がみつからないのでよくわからない)。

Principal Component Analysis and EOF analysis

Meteorology people use EOF analysis. Remote sensing people use Principal Component Analysis. I realized that those two are the same thing.

In PCA, eigenmatrix analysis is performed for covariance matrix C. In EOF, singular vector decomposition is performed for the data matrix Z (usually not square). But according to a MIT lecture note, singular vector decomposition is just taking eigenmatrix of ZTZ, which should be called covariance matrix (although not normalized).

conversion equations for Oblique Lambert Projection

地図投影法の本 (英語)

地図投影法の本は日本語ではいいのがない。とりあえず、身の回りでよく使われているのでメモ。

John Parr Snyder, 1987: "Map Projections - A Working Manual". US Geological Survey Professional Paper 1395, Washington DC.

この本は球面じゃなくて回転楕円体上での投影変換の式をさぼらずに書いているので利用価値が高いわけだけれど、斜軸ランベルト正角図法については、USGS での利用例がないらしく、球面での式しか載っていない。

うーん困った。どこかに国土が斜めに延びている国はないだろうか。

Win32 API チュートリアル

ペゾルドの本なんかいまどき本屋に行っても売っていないので、ぐぐってみつけたもの。Goodbye, cruel world! だってさ。

[AOS650] 客観解析の手法

数値天気予報 (NWP; Numerical Weather Prediction) は初期値として現在の大気の状態を与え、時間微分を含んだ微分方程式として与えられる初期値問題を解いて将来の大気の状態を予測するものである。 電子計算機によって微分方程式を数値的に解くために規則的に並んだ格子の上での初期値が必要となる。ところが気象観測は(1)正確に格子の上で行われるわけではないし(2)観測がない格子もあるところから、初期値をどうやって作るかが問題となる。 人間がデータを眺めながら適当に初期値を推定することも可能(昔リチャードソンがやったように)であるが、実用的な速度を得るためには機械的に推定するアルゴリズムが必要である。 この作業を客観解析(Objective Analysis)という。過去の数値予報結果と観測値を併用する場合をデータ同化 (Data Assimilation)という。

クレスマン法

まず簡単に重みつき平均を行うことを考える。それぞれの格子に対して近い観測値ほど大きく影響するだろうと考えると、 解析値は解析方程式 analysis equation

xaj = xbj + {[Σi=1M Wi,j (yi - xbj)] / (Σi Wi,j)} ... (1)
のように与えることになるだろう。ここで iは観測点 (M 個) の番号、 jは格子点 (N 個) の番号、 xa は解析値、 xb は推定値 guess (数値予報のためには過去の数値予報から得るが、 目的によっては気候値でも何でもいい)、 y は観測値、 Wi,jは 観測点iと格子点jの位置関係で決まる重み係数である。 ついでだから用語を導入する。上式の { ... } を 解析インクリメント analysis increment, その分母の中にある [yi - xbj] はイノベーション innovation と呼ばれる。

ここで、

Wi,j = max[0, (R2 - di,j) / (R2 + di,j)] ... (2)
のように定義するものをクレスマン法 Cressman scheme という。 di,jは観測点iと格子点jの距離、 Rは影響半径 influence radius (つまり観測値が格子点に影響を及ぼす範囲) である。 Rはデータの配置から経験的に決めるほかなく、 クレスマン法というとは解析結果を見ながらくり返し R を小さくしていく過程を含める。

クレスマン法はそれなりにうまく働くが、問題点は (a) 推定値の品質の良否にかかわらず、観測がありさえすれば それに引き摺られてしまうこと、 (b) 観測だけで推定値がないと使えないこと、 (c) 解析値が地衡流平衡・静水圧平衡などの「望ましい関係」から離れていって しまっても止める術がないこと、である。

一応参照: G.P. Cressman, 1959: An operational objective analysis system. Mon. Wea. Rev., 87, 367-374.

バーンズ法

数値予報からみると脇道にそれるが、観測だけで何とかすることを考える (前節問題 (a))。 どうすればいいかというと、推定値も観測から作ってしまえばいいのである。

xbj = {Σi W0i,j yi} / {Σi W0i,j} ... (3)
ここで、重み関数は (2) 式と違って遠距離でもかすかには効力を持ってほしいので (本当の理由は後述) 指数関数型
W0i,j = exp[- di,j2/κ] ... (4)
とする。ここでκは長さの二乗の次元を持つ定数である。 次に、この推定値と一緒に再度観測値を使って (2) を使う、 ただし重み関数は短距離に集中するように
W1i,j = exp[- di,j2/(γκ)] ... (5)
(0 < γ < 1) というのがバーンズ法 Barnes objective analysis scheme である。

一応参照: S.L. Barnes, 1964: a technique for maximizing details in numerical weather map analysis. J. Appl. Meteor., 3, 396-409.

周波数応答の検討

さて、どうして2回も補間のような式を使うのだろうか。 重み関数を工夫したら1回で済まないのだろうか。 この問題に答えるために、世界が1次元であると考えて (ここまでのところ、空間が2次元であるか3次元であるかについては 何も語っていないことに注意されたい。何次元でも使えるのだ) 正弦波型の場

f(x) = Asin(2πx/λ) ... (6)
が等間隔に観測された場合を考えてみる。 するとまず推定値は正弦波型となり、その振幅は
D0(λ) = exp(-π/λ)2 ... (7)
となる。補間法によってデータが「ぼやける」ので、 必然的にローパスフィルタであるわけだ。 次のステップを含めると
D0+1(λ) = D0 + (1 - D0) D0γ ... (8)
のようになる。そこで、γの値を調整 (0.2 < γ < 1) すれば結構な周波数応答プロファイルが 作れるというわけである。 詳細は S.E. Koch et al., 1983: an interactive Barnes objective map analysis scheme for use of satellite and conventional data. J. Climate Appl. Meteor. 22, 1487-1503. の 1492 ページに載っている図を見ていただきたい。 どうも、バーンズは本来 (2) 式をくり返し使うというスキームを発表したのだが、 Koch たちが GEMPAK で実装するために 1回だけ使うように変えたということらしい。

最適内挿法

さて、数値予報の文脈に戻ってクレスマン法の問題 (a) を考える。 要するに最良の解析値が欲しいということなので、 何が「良」であるかを定量的に決めさえすれば最小二乗法で答えが求まる。 観測・推定値にバイアスがないとすれば、 最小化したいのは観測・推定値からの差の二乗を 観測・推定値の誤差の分散で割ったものの和であろう。

まず簡単のため、1地点で1変数の観測・推定が行われることを考える。 このとき最小化したい関数 (費用関数 cost function) は

J = (y - xa)2o2 + (xb - xa)2b2 ... (9)
である。ここで σb が推定値の、 σo が観測値の誤差の標準偏差である。 最小の必要条件は
dJ/dxa = 2(y - xa)/σo2 + 2(xb - xa)/σb2 = 0
から
xa = (σb2y + σo2xb) / (σb2 + σo2)
となり、 クレスマン解析方程式 (2) 風に書き直すことができて
xa = xb + [σb2 / (σb2 + σo2)] (y - xb) ... (10)
である。

実際的な多変数、多格子、格子点と一致しない観測点の場合も 変数が行列になるだけで本質的にはこれと同様である。

まず xaxbN 要素のベクトル、 yM 要素のベクトルとする。 今まで y - x と書いて来たところは 観測点が格子点と一致しないので補間が必要となり y - Hx となる。 ここで H は NM 列のベクトルで、 i 番目の列が 観測点 i の値を得るための補間係数のベクトルで ΣjN hij = 1 となるようなものである (なお、たとえば衛星からみた輝度といった予報変数でない観測を使う場合、 補間行列 H は診断的モデルになるであろう)。 すると費用関数は

J = (y - Hxa)T R-1 (y - Hxa) + (xa - Hxb)T B-1 (xa - Hxb) ... (11)
のように書ける。 ここで R は観測の誤差の共分散行列 (rij が観測点 i, j の誤差の 共分散)、 B は推定の誤差の共分散行列である。 すると
dJ/dxa = HTR-1Hxa + [(Hxa)T R-1H]T - HTR-1y - (yT R-1H)T + B-1xa + (xaTB-1)T - B-1xb - (xbTB-1)T = 0
から
(HTR-1H + B)xa + [(Hxa)T R-1H + xaTB-1]T = HTR-1y + B-1xb + (yT R-1H + xbTB-1)T
となる。ここで HTH, R, B が対称行列であることを利用して
(HTR-1H + B)xa = HTR-1y + B-1xb
であることが示せて、
xa = xb + (HTR-1H + B)-1 HTR-1 (y - Hxb) ... (12)
または
xa = xb + BHT (HTBH + R)-1 (y - Hxb) ... (13)
を得る。 これを最適内挿法 optimal interpolation という。

変分法

力学的バランスをなるべく保った解析を行いたい場合、 費用関数にバランスからのずれに応じて大きくなる項を追加する。 すると解析的に解くのが難しくなるので反復法で いろいろ xa を変えてみて、ボチボチ J が 小さくなったらよしとする。 これを三次元変分法 three-dimensional variational method という。 さらにいろいろの時間の観測値を用いたい場合、 費用関数の観測値の項を時間積分の形にする。 これを四次元変分法 four-dimensional variational method という。

参照: T.W.Schlatter, 2000: Variational assimilation of meteorological observations in the lower atmosphere: a tutorial on how it works. J. Atmos. Solar-Terrestrial Phys., 62, 1057-1070.

g77 プロジェクト終焉

だいぶん前からフリーの Fortran 95 のコンパイラ g95 が使えるようになってきたのだけれど、これとは別物の gfortran というのがあって、長年親しまれた FORTRAN 77 コンパイラ g77 の代わりに今度の gcc 3.5 から GCC (Gnu Compiler Collection) の Fortran 95 コンパイラとして用いられるようになるということです。



ついに g77 頼みの時代が終わったというのは時代の進歩ですが、このニュースを聞いてちょっと不安になりました。というのは、g77 は意地悪で数値型の種別型パラメタがバイト数ではなくてワード数になっていたのです。
たとえば世間のコンパイラで

INTEGER(2):: I
INTEGER(4):: J
INTEGER(8):: K

と書くところを

INTEGER(5):: I ! 0.5 ワードなんだって ! ちなみに i386 の g77 では使えない
INTEGER(1):: J
INTEGER(2):: K

と書かせるようになっていたのです。種別型パラメタが規格上バイト数とは別物であることを判ってもらうためとはいえ、単に意地悪としか言いようがない仕様ですが、ひょっとして、gfortran がこの意地悪を継承しているのだとすると、バイト数を決め打ちしなくてはならないアプリケーションは悲惨、と思ったのですが、まだほんとうに作りかけのマニュアルを斜め読みする限りさすがにそんなことはないようです。

つまり、今回の g77 滅亡で、「種別型パラメタはバイト数」というデファクトスタンダードが確立したことになります。

ところで、g95 プロジェクトの人は見放されて落胆しているのじゃないか、とおもって見ていたら、g95 やってる Andy ちゃっかり gfortran の創始者しているのではないですか。なんだこれは。

[Fortran] Winsock worked!



module winsock_types
implicit none

type wsadata
integer(2):: version
integer(2):: highversion
character(257):: szDescription
character(129):: szSystemStatus
integer(2):: maxSockets
integer(2):: maxUdpdg
integer(4):: lpVenderinfo
end type

type sockaddr_in
sequence
integer(2):: sin_family
integer(2):: sin_port
integer(4):: sin_addr
integer(4):: sin_zero(2)
end type

type socket_t
integer(4):: ival
end type

type hostent_t
integer(4):: h_name
integer(4):: h_aliases
integer(2):: h_addrtype
integer(2):: h_length
integer(4):: h_addr_list
end type

end module

module winsock
use winsock_types
implicit none
interface

integer(4) function wsastartup(reqver, wsainfo)
!dec$ attributes dllimport, alias: '_WSAStartup@8' :: WSAStartup
use winsock_types
integer(2), intent(in):: reqver
!dec$ attributes value:: reqver
type(wsadata), intent(out):: wsainfo
end function

integer(4) function wsacleanup()
!dec$ attributes stdcall, dllimport, alias: '_WSACleanup@0' :: WSACleanup
end function

integer(4) function wsagetlasterror()
!dec$ attributes stdcall, dllimport, alias: '_WSAGetLastError@0' :: WSAGetLastError
end function

function socket(af, type, protocol) result(result)
!dec$ attributes stdcall, dllimport, alias: '_socket@12' :: socket
use winsock_types
integer(4), intent(in):: af
integer(4), intent(in):: type
integer(4), intent(in):: protocol
type(socket_t):: result
end function

integer(4) function closesocket(s)
!dec$ attributes stdcall, dllimport, alias: '_closesocket@4' :: closesocket
use winsock_types
type(socket_t), intent(in):: s
end function

end interface

integer(4), parameter:: AF_UNIX = 1
integer(4), parameter:: AF_INET = 2
integer(4), parameter:: AF_INET6 = 23
integer(4), parameter:: SOCK_STREAM = 1
integer(4), parameter:: SOCK_DGRAM = 2
type(socket_t), parameter:: INVALID_SOCKET = socket_t(not(0))
integer(4), parameter:: SOCKET_ERROR = -1
type(sockaddr_in), parameter:: SOCKADDR_IN_ANY = sockaddr_in(AF_INET, 0, &
& 0, (/2 * 0/))

interface operator(==)
module procedure socket_eq
end interface

interface

integer(4) function bind(s, name, namelen)
!dec$ attributes stdcall, dllimport, alias: '_bind@12' :: bind
use winsock_types
type(socket_t), intent(in):: s
type(sockaddr_in), intent(in):: name
!dec$ attributes reference:: name
integer(4), intent(in):: namelen
end function

function accept(s, addr, addrlen) result(result)
!dec$ attributes stdcall, dllimport, alias: '_accept@12' :: accept
use winsock_types
type(socket_t), intent(in):: s
type(sockaddr_in), intent(out):: addr
!dec$ attributes reference:: addr
integer(4), intent(out):: addrlen
type(socket_t):: result
end function

! GUISE: return value is pointer to HOSTENT structure
integer(4) function gethostbyname(szname)
!dec$ attributes stdcall, dllimport, alias: '_gethostbyname@4' :: gethostbyname
character(*), intent(in):: szname
!dec$ attributes reference:: szname
end function

integer(4) function inet_addr(szname)
!dec$ attributes stdcall, dllimport, alias: '_inet_addr@4' :: inet_addr
character(*), intent(in):: szname
!dec$ attributes reference:: szname
end function

integer(4) function connect(s, name, namelen)
!dec$ attributes stdcall, dllimport, alias: '_connect@12' :: connect
use winsock_types
type(socket_t), intent(in):: s
type(sockaddr_in), intent(in):: name
!dec$ attributes reference:: name
integer(4), intent(in):: namelen
end function

end interface

interface send

integer(4) function sendc(s, buf, len, flags)
!dec$ attributes stdcall, dllimport, alias: '_send@16' :: sendc
use winsock_types
type(socket_t), intent(in):: s
character(*), intent(in):: buf
!dec$ attributes reference:: buf
integer(4), intent(in):: len, flags
end function

integer(4) function send4(s, buf, len, flags)
!dec$ attributes stdcall, dllimport, alias: '_send@16' :: send4
use winsock_types
type(socket_t), intent(in):: s
integer, intent(in):: buf(*)
!dec$ attributes reference:: buf
integer(4), intent(in):: len, flags
end function

end interface

interface recv

integer(4) function recvc(s, buf, len, flags)
!dec$ attributes stdcall, dllimport, alias: '_recv@16' :: recvc
use winsock_types
type(socket_t), intent(in):: s
character(*), intent(out):: buf
!dec$ attributes reference:: buf
integer(4), intent(in):: len, flags
end function

integer(4) function recv4(s, buf, len, flags)
!dec$ attributes stdcall, dllimport, alias: '_recv@16' :: recv4
use winsock_types
type(socket_t), intent(in):: s
integer, intent(out):: buf(*)
!dec$ attributes reference:: buf
integer(4), intent(in):: len, flags
end function

end interface

contains

logical function socket_eq(s1, s2) result(result)
type(socket_t), intent(in):: s1, s2
result = s1%ival == s2%ival
end function

type(sockaddr_in) function make_sockaddr_in(name, port) result(result)
character(*), intent(in):: name
integer, intent(in):: port
character(2048):: namebuf
integer:: namelen
integer(4):: addr, ival
type(HOSTENT_T):: hostent
integer(1):: bval(4)
integer:: i
pointer(addr, hostent)
pointer(addr, ival)
pointer(addr, bval)
namebuf = name
namelen = min(len(namebuf), len_trim(namebuf) + 1)
namebuf(namelen:namelen) = char(0)
addr = inet_addr(namebuf)
if (addr /= not(0_4)) then
result%sin_addr = addr
else
addr = gethostbyname(namebuf)
if (addr == 0) then
result%sin_addr = 255
result%sin_port = 0
return
endif
addr = hostent%h_addr_list
addr = ival
result%sin_addr = ival
endif
result%sin_family = AF_INET
! result%sin_port = port
result%sin_port = ior(ishft(iand(255, port), 8), iand(ishft(port, -8), 255))
result%sin_zero(:) = 0
end function

character(60) function wsastrerror(errno) result(result)
integer, intent(in), optional:: errno
integer:: error_code
if (present(errno)) then
error_code = errno
else
error_code = wsagetlasterror()
endif
select case(error_code)
case(10004); result = "interrupted function call"
case(10009); result = "EBADF"
case(10013); result = "permission denied"
case(10014); result = "bad address"
case(10022); result = "invalid function argument"
case(10024); result = "too many open files"
case(10035); result = "resource temporarily unavailable"
case(10036); result = "operation now in progress"
case(10037); result = "operation already in progress"
case(10038); result = "socket operation on non-socket"
case(10039); result = "destination address required"
case(10040); result = "message too long"
case(10041); result = "protocl wrong type for socket"
case(10042); result = "bad protocol option"
case(10043); result = "protocol not supported"
case(10044); result = "socket type not supported"
case(10045); result = "operation not supported"
case(10046); result = "protocol family not supported"
case(10047); result = "address family not supported by protocol family"
case(10048); result = "address already in use"
case(10049); result = "cannot assign requested address"
case(10050); result = "network is down"
case(10051); result = "network is unreachable"
case(10052); result = "network dropped connection on reset"
case(10053); result = "software caused connection abort"
case(10054); result = "connection reset by peer"
case(10055); result = "no buffer space is available"
case(10056); result = "socket is already connected"
case(10057); result = "socket is not connected"
case(10058); result = "cannot send after socket shutdown"
case(10059); result = "ETOOMANYREFS"
case(10060); result = "connection timed out"
case(10061); result = "connection refused"
case(10062); result = "ELOOP"
case(10063); result = "ENAMETOOLONG"
case(10064); result = "host is down"
case(10065); result = "no route to host"
case(10066); result = "ENOTEMPTY"
case(10067); result = "too many processes"
case(10068); result = "EUSERS"
case(10069); result = "EDQUOT"
case(10070); result = "ESTALE"
case(10071); result = "EREMOTE"
case(10091); result = "network subsystem is unavailable"
case(10092); result = "WINSOCK.DLL version out of range"
case(10093); result = "successful WSAStartup not yet performed"
case(10094); result = "graceful shutdown in progress"
case(10101); result = "DISCON"
case(10102); result = "NOMORE"
case(10103); result = "CANCELLED"
case(10104); result = "invalid procedure table from service provider"
case(10105); result = "invalid service provider version number"
case(10106); result = "unable to initialize a service provider"
case(10107); result = "system call failure"
case(10108); result = "SERVICE_NOT_FOUND"
case(10109); result = "class type not found"
case(10110); result = "NO_MORE"
case(10111); result = "CANCELLED"
case(10112); result = "REFUSED"
case(11001); result = "authoritative DNS answer: host not found"
case(11002); result = "non-authoritative: host not found; or server fail"
case(11003); result = "non-recoverable DNS error FORMERR, REFUSED, or NOTIMP"
case(11004); result = "no DNS data record of requested type"
case default
result = ""
write(unit=result, fmt="('unknown Winsock error ',i12.1)") error_code
end select
end function

character(1024) function sz2char(sz) result(result)
character(*), intent(in):: sz
character:: c
integer:: i, j
j = 1
do, i = 1, len(sz)
if ((j + 4) > len(result)) exit
c = sz(i:i)
select case(ichar(c))
case(0)
exit
case(92)
result(j:j+1) = c // c
j = j + 2
case(10)
result(j:j+1) = char(92) // 'n'
j = j + 2
case(13)
result(j:j+1) = char(92) // 'r'
j = j + 2
case(1:9, 11, 12, 14:31, 127:255)
write(result(j:j+3), "(A2, Z2.2)") char(92) // 'x', ichar(c)
j = j + 4
case default
result(j:j) = c
j = j + 1
end select
enddo
if (j < len(result)) then
result(j: ) = ''
endif
end function

end module

subroutine main
use winsock
implicit none

type(wsadata):: wsainfo
integer(2):: reqver
integer(4):: i
type(socket_t):: s
type(sockaddr_in):: sa
character(1024):: buf

reqver = 514
i = wsastartup(reqver, wsainfo)
print *, "WSAStartup =", i, wsainfo%version, wsainfo%highversion, &
& "'"//trim(sz2char(wsainfo%szDescription))//"'("//trim(sz2char(wsainfo%szSystemStatus))//")"
if (i /= 0) then
stop 16
endif

s = socket(AF_INET, SOCK_STREAM, 0)
print *, 'socket', s
if (s == INVALID_SOCKET) then
print *, 'socket:', wsastrerror()
goto 900
endif

sa = make_sockaddr_in("www.asahi.com", 80)
print *, ibits(sa%sin_addr, 0, 8), ibits(sa%sin_addr, 8, 8), &
& ibits(sa%sin_addr, 16, 8), ibits(sa%sin_addr, 24, 8), &
& ibits(sa%sin_port, 0, 8), ibits(sa%sin_port, 8, 8)
if (sa%sin_port == 0) goto 900

i = connect(s, sa, 16)
print *, 'connect', i
if (i == SOCKET_ERROR) then
print *, 'connect:', wsastrerror()
goto 910
endif

i = send(s, "GET /" // char(13) // char(10), 7, 0)
if (i == SOCKET_ERROR) then
print *, 'send:', wsastrerror()
goto 910
endif

do
i = recv(s, buf, len(buf), 0)
if (i == SOCKET_ERROR) then
print *, 'recv:', wsastrerror()
goto 910
endif
if (i == 0) exit
print *, buf(1:i)
enddo

910 continue
print *, 'closesocket', closesocket(s)
print *, "WSACleanup =", wsacleanup()
print *, "OK"
return

900 continue
print *, "WSACleanup =", wsacleanup()
print *, "OK"

end subroutine

call main
end

整数除算による単位の換算

範囲や精度があらかじめ判っている実数を保存するのに、整数をつかうことがあります。たとえば気温が 1℃の位でしか測れない状況で、まあ常識的にいって -100 .. 100 くらいの範囲をサポートしておけばいいだろうというと、1バイト(1オクテット)の整数で保存できます。これが float だったら4倍のファイル容量を食ってしまうところだから、大だすかり。

さて、そんなデータが幾種類かあって、単位が違っている(温度の摂氏華氏とか、風速のメートル毎秒と海里毎時とか、気圧のヘクトパスカル対水銀インチとか....決して特定の国を名指しして批難したいわけではないのですが、でも迷惑なものは迷惑だ)という悲惨なときに、単位の換算をどうするか。どうするかって、普通は何も考えないで定義式そのまんま

integer(1):: c, f

c = (dble(f) - 32.0) / 1.8

とかするわけです。しかし、整数データから整数データを作るのに、浮動小数点演算はないだろう、と思うわけです。PowerPC だったら知りませんが Intel だったら時間がかかって仕方がない。それで、素朴には

c = (f - 32) * 10 / 18

でよさそうな気がします。が、ダメです。整数除算の切捨てがあるので、多くの場合1度低い温度を出してしまいます。華氏温度が1度単位99度以下であるとして、最適なのは、

c = (f - 32) * 51 / 92

です。華氏温度が 0.1 度単位であるとすると、

c = (f - 320) * 551 / 992

です。

そういう係数を決めるためのプログラムです。(FORTRAN 77)
FACTOR に実数だったら使いたかった係数、MAXVAL に変換対象の最大の数値、をいれます。
負値には対応していません。


program convprx
implicit none
double precision factor
integer maxval, try1st
data factor / 3.1415927 /
data maxval / 999 /
data try1st / 0 /
namelist /param/ factor, maxval, try1st
integer iostat
write(*, nml=param)
read(*, nml=param, iostat=iostat)
write(*, nml=param)
if (factor <= 0.0) stop 'FACTOR MUST BE POSITIVE'
if (maxval <= 0.0) stop 'MAXVAL MUST BE POSITIVE'
call optima(factor, maxval, try1st)
end program

subroutine optima(factor, maxval, try1st)
implicit none
double precision factor
integer maxval, try1st
integer div, kd, ko, i
dimension kd(4), ko(4)
if (try1st .ne. 0) then
div = try1st
call trydiv(factor, maxval, div, ko)
print '(4(I7.1))', ko(1), ko(2), ko(3), ko(4)
if (ko(4) .eq. 0) goto 4001
else
ko(4) = maxval + 1
endif
div = 0
4000 continue
div = div + 1
if (div .eq. 1 000 000) goto 4001
call trydiv(factor, maxval, div, kd)
if (kd(4) .le. min(ko(4), maxval - 1)) then
print '(4(I7.1))', kd(1), kd(2), kd(3), kd(4)
do, i = 1, 4
ko(i) = kd(i)
enddo
endif
if (kd(4) .eq. 0) goto 4001
goto 4000
4001 continue
if (ko(4) .ne. 0) then
print *, '--- give up ---'
endif
print *, 'RESULT: Y = (', ko(1), ' * X +', ko(3), ' )/ ', ko(2)
end subroutine

subroutine trydiv(factor, maxval, div, kd)
implicit none
double precision factor
integer maxval, div, kd
integer i, mul, ofs, err
dimension kd(4)
mul = int(factor * div + 0.5)
ofs = 0
do, i = 1, 4
kd(i) = maxval + 1
enddo
4100 continue
call geterr(factor, maxval, mul, div, ofs, err)
if (err .lt. kd(4)) then
kd(1) = mul; kd(2) = div; kd(3) = ofs; kd(4) = err
if (err .eq. 0) return
else if (kd(4) .le. maxval .and. err .gt. kd(4)) then
return
endif
if (err .gt. maxval) return
if (ofs .ge. (div / 2)) return
ofs = ofs + 1
goto 4100
end subroutine

subroutine geterr(factor, maxval, mul, div, ofs, err)
implicit none
double precision factor
integer maxval, mul, div, ofs, err
integer i, appro, true, dif
err = 0
do, i = 0, maxval
appro = (mul * i + ofs) / div
true = int(dble(i) * factor + 0.5d0) ! 訂正したよ !
dif = appro - true
if (abs(dif) .ge. 2) then
err = maxval + 1
return
endif
if (dif .ne. 0) then
err = err + 1
endif
enddo
end subroutine


(追記 2/9 03:21) プログラムにバグがあったので訂正した。

DF から Win32 API を呼ぶ

わかった。引数が値渡しになっていなかったのが原因。

こうすれば動作する。


module winsock_types
implicit none

type wsadata
integer(2):: version
integer(2):: highversion
character(257):: szDescription
character(129):: szSystemStatus
integer(2):: maxSockets
integer(2):: maxUdpdg
integer(4):: lpVenderinfo
end type

type sockaddr_in
integer(2):: sin_family
integer(2):: sin_port
integer(1):: sin_addr(4)
integer(1):: sin_zero(8)
end type

end module

module winsock
use winsock_types
implicit none
interface

integer(4) function wsastartup(reqver, wsainfo)
!dec$ attributes dllimport, alias: '_WSAStartup@8' :: WSAStartup
use winsock_types
integer(2), intent(in):: reqver
!dec$ attributes value:: reqver
type(wsadata), intent(out):: wsainfo
end function

integer(4) function wsacleanup()
!dec$ attributes stdcall, dllimport, alias: '_WSACleanup@0' :: WSACleanup
end function

integer(4) function wsagetlasterror()
!dec$ attributes stdcall, dllimport, alias: '_WSAGetLastError@0' :: WSAGetLastError
end function

integer(4) function socket(af, type, protocol)
!dec$ attributes stdcall, dllimport, alias: '_socket@12' :: socket
use winsock_types
integer(4), intent(in):: af
integer(4), intent(in):: type
integer(4), intent(in):: protocol
end function

integer(4) function closesocket(s)
!dec$ attributes stdcall, dllimport, alias: '_closesocket@4' :: closesocket
integer(4), intent(in):: s
end function

end interface

integer(4), parameter:: AF_UNIX = 1
integer(4), parameter:: AF_INET = 2
integer(4), parameter:: AF_INET6 = 23
integer(4), parameter:: SOCK_STREAM = 1
integer(4), parameter:: SOCK_DGRAM = 2
integer(4), parameter:: INVALID_SOCKET = not(0)
integer(4), parameter:: SOCKET_ERROR = -1

interface

integer(4) function bind(s, name, namelen)
!dec$ attributes dllimport, alias: '_bind@12' :: bind
use winsock_types
integer(4), intent(in):: s
type(sockaddr_in), intent(in):: name
integer(4), intent(in):: namelen
end function

integer(4) function accept(s, addr, addrlen)
!dec$ attributes dllimport, alias: '_accept@12' :: accept
use winsock_types
integer(4), intent(in):: s
type(sockaddr_in), intent(out):: addr
integer(4), intent(out):: addrlen
end function

end interface

contains

character(1024) function sz2char(sz) result(result)
character(*), intent(in):: sz
character:: c
integer:: i, j
j = 1
do, i = 1, len(sz)
if ((j + 4) > len(result)) exit
c = sz(i:i)
select case(ichar(c))
case(0)
exit
case(92)
result(j:j+1) = c // c
j = j + 2
case(10)
result(j:j+1) = char(92) // 'n'
j = j + 2
case(13)
result(j:j+1) = char(92) // 'r'
j = j + 2
case(1:9, 11, 12, 14:31, 127:255)
write(result(j:j+3), "(A2, Z2.2)") char(92) // 'x', ichar(c)
j = j + 4
case default
result(j:j) = c
j = j + 1
end select
enddo
if (j < len(result)) then
result(j: ) = ''
endif
end function

end module

subroutine main
use winsock
implicit none

type(wsadata):: wsainfo
integer(2):: reqver
integer(4):: s

reqver = 514
s = wsastartup(reqver, wsainfo)
print *, "WSAStartup =", s, wsainfo%version, wsainfo%highversion, &
& "'"//trim(sz2char(wsainfo%szDescription))//"'("//trim(sz2char(wsainfo%szSystemStatus))//")"
print *, ''

s = socket(AF_INET, SOCK_STREAM, 0)
print *, 'socket', s
if (s == INVALID_SOCKET) then
print *, 'error', wsagetlasterror()
goto 900
endif

print *, 'closesocket', closesocket(s)
print *, "WSACleanup =", wsacleanup()
print *, "OK"
return

900 continue
print *, "WSACleanup =", wsacleanup()
print *, "OK"

end subroutine

call main
end

CEE556 Lec2

Ratio Transformation



  • Numeric
  • Angle
  • Difference
  • Normalized
  • Log -- loge(Sn/Sn+1)

Difference subtle - not so much.

use - atmospheric effect


Absorption (A) and Scattering (B).
The atmospheric effects actually depends on the wavelength.

S(λ) = A S0(λ) + B(λ)

Scattering can be corrected by examining the darkest pixels in schene.

color display of angle ratios etc.


2 vs 1 -> ch1, 3 vs 2 -> ch2, 1 vs 3 -> ch3.

Difference Ratio


In order to avoid to handle negative values, usually form

scale * (1 + [(IR - Red) / (IR + Red)])

is used.

FC2Ad

上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。