2008 年度 特別研究

u-w 法を用いた樹木の成長解析ソフトウェア開発と検証

2008 年１月 23 日
龍谷大学理工学部
環境ソリューション工学科
T050592
青木 将晃
宮浦 富保 教授

目次
要旨・・・・・・・・3P
1.はじめに・・・・・4P
2.u-w 法について・・4P
3.目的・・・・・・・6P
4.方法・・・・・・・9P
5.結果・・・・・・・8P
6.考察・・・・・・・12P
7.引用文献・・・・・14P
8.謝辞・・・・・・・14P
9.付録・・・・・・・15P2

u-w 法を用いた樹木の成長解析ソフトウェア開発と検証
T050592
宮浦

青木
富保

将晃
教授

１．はじめに
樹木の成長はとてもゆっくりとしたものであるが、成長の記録を年輪として幹の内側に残しながら着実に大きくなっている。樹木の年輪を調べ、その成長の仕方を曲線で近似することにより、
ある程度過去の森林の姿を再現できたり、森林の将来の状態を予測したりできると考えられる。
現在までに数多くの成長曲線が提唱されてきたが、どれを選べばよいのかの判断基準は存在せず、それらの成長曲線では実際のデータをうまく近似出来ない場合が多かった。
Hozumi(1985)が提唱した u-w 法で得られる成長曲線は従来の成長曲線の多くを包含しており、実際の成長データのあてはめが従来の成長曲線式よりも柔軟に行える。u-w 法では以下のように定義される u と w の関係をプロットして近似曲線を求める。
w m
−c
d 1
1 dw
) ) ただし u =u = Aw (１―(w:幹の材積（㎥）=− ( )
Ww 2 dtdt w
しかし、u-w 法の適用には大変な手間を必要とするため、u-w 法による計算を効率よく正確に行い、成長曲線の解析に役立てることを目的として u-w 法の計算ソフトウェアを作成し、その検証
を行った。
２．方法
Visual Basic 2008 というプログラム言語を使って、Windows アプリケーションとして成長解析
ソフトウェアを作成した、また 2007 年に「龍谷の森」のコナラ二次林で伐採して得られた実際の年輪データを使用して成長解析ソフトウェアの有効性を検証した。
３．結果と考察
図 1 は、開発した成長解析ソフトウェアの実行例である。データから u と w の値を計算し、それを元に自動的に近似曲線を引くことができる。
今回開発したソフトウェアでは A,c,W それぞれを一定量ずつ変化させて u の推定値と実測値の残差平方和が最小になる値を探すという方法で近似曲線を探した。多くの樹木で近似できたが、一部うまく近似できないものもあった。
定数を解析した結果、c は 1～1.4 の範囲をとることが多かった。A と c には強い負の相関が見られた。
実際の樹木の成長は複数の成長曲線で近似されることが多く、この成長曲線が変わる現象を成長曲線の乗り換えと呼ぶ。成長曲線の乗り換えが起こるごとに W が大き
図 1,今回開発した u-w 図を描いて計算する
くなる傾向がほとんどの樹木で見られた。また成長曲線
ソフトウェアの乗り換えの年代を調べた結果、樹種ごとに特徴がある
ことが分かった。コナラでは 1970 年代に乗り換えが集中していた。
４．結論
ソフトウェア開発により、u-w 法による樹木の成長解析が効率的に行えるようになった。
u-w 曲線式は定数の変動幅が大きく、今回開発したソフトウェアで使ったアルゴリズムでは全ての曲線式に対応できなかったため今後の課題となった。
c は 1～1.4 の範囲をとることが多かったが、これは多くの場合ロジスティック曲線、指数関数
曲線、ミッチャーリッヒ曲線などの有名な成長曲線式のいずれにも当てはまらないと言えるので、
u-w 法を使った方がより正確な成長曲線を導き出せると考えられる。
u-w 図において複数の成長曲線で分けたほうがきれいに近似できるので、多くの樹木は成長曲
線の乗り換えを何度も行いながら成長していると考えられる。
3

1.はじめに樹木の成長はとてもゆっくりとしたものであるが、成長の記
録を年輪として幹の内側に残しながら着実に大きくなっている。
樹木の年輪を調べ、その成長の仕方を曲線で近似することにより、ある程度過去の森林の姿を再現できたり、森林の将来の状態を予測したりできると考えられる（宮浦,2002）。
個体サイズの変化の時間的経過を定量的にあらわす曲線を、一般的に成長曲線とよぶ。樹木の成長のしかたを記述する方法

図１.成長曲線の例

として、現在までに数多くの成長曲線が提唱されてきた。例え

日本生態学会（2004）より

ば実際によく使われるものとしては図１のように指数関数的に成長する指数曲線や、途中まで指数関数的に成長し上限値に近づくと成長が止まっていくロジスティック曲線などがある。図 1 で
t は時間、N は個体数や個体サイズなど、K は N の上限値である。しかし数多くある成長曲線のうちどれを選べばよいのかという点についての合理的な判断基準は存在せず、それらの成長曲線
では実際のデータをうまく近似出来ない場合が多かった（宮浦,1990）。

2.u-w 法について
u-w 法は一本の木の幹材積の成長解析方法として Hozumi(1985)によって提案された。
ここでいくつかの変数を定義する。
t：時間[year]

w：幹材積(幹部分の体積)[㎥]

W：w の上限値[㎥]

v：成長速度、

dw
dt

−1

]

・year

−1

[㎥・year

1 dw
w dt
s
v
d 1
u：
= 2 = − ( )w
dt ww
s：相対成長速度、

[year
[m

−3

−1

]
]

u-w 法では、u をｗの関数として次式で表わす（Hozumi,1985）

u = Aw − c (１―(w m) ) ・・・(1)W

u と w の関係を表す u-w 曲線は従来の成長曲線の多くを包含するため、実際の成長データのあてはめが従来の成長曲線式よりも柔軟に行える。また生物的な根拠をもった方法であり、あてはめの結果得られる係数や定数はなんらかの生物的な意味を持っていることが期待される(宮浦,1990)。
この u-w 曲線式の A、c、m は定数である。また W は w の上限値である。(1)式の m,c,W の取り方によっては指数曲線、ロジスティック曲線、ミッチャーリッヒ曲線などの成長方程式を得る
ことができる。

1

実際の使い方の例としては以下の通りである。すなわち（１）式において例えば W=1 に固定しておき、m=1～２、c=-2～2 程度の範囲で w と c の値を何通りか変えて、両対数グラフ上に u-w
曲線を描く、この曲線をトレーシングペーパーに転写したものを定規として用いる（図 4）。次に表１のように年輪の樹齢・幹材積のデータから算出した u と w を両対数グラフ上にプロットし、作成した定規を重ねて最もよく適合する曲線を選ぶ。そして選んだ曲線から残りの定数である Wと A を算出して成長曲線式を決定する。
図 2 は u-w 図をプロットした例である、実際にはこのように一本だけでなく複数の曲線で近似することが多いが、この曲線が変わる現象は成長曲線の乗り換えと呼ばれる。図 3 は幹材積の成
長経過を片対数グラフにプロットしたものである。図 2 の方が成長曲線の乗り換えがはっきりとわかる。
u-w 法を用いて得られる係数や定数の値を、気象条件や林分構造の側面から吟味することにより、樹木の成長に関する内的な要因と外的要因の効果をある程度分離して理解できると期待されている（宮浦,1990）。

u

w

t

w

図 3.時間ごとの幹材積

図 2.u-w 図

表 1.u と w の計算例
樹齢

幹材積[㎥]

15

0.00092

20
25
30

u [m

−3

・year

−1

w [㎥]

]

139.3

0.0026121.9

0.01097

6.42

0.03450.00567
0.018850.05706

図 4.c=1.1、1.2、1.3…と変化させたグラフを描いた紙

u と w の値は年齢を t とすると t2-t1 の間曲線の形を決めるときに使う。
2

の近似値であり、u≒u、w≒w である
3.目的
u-w 法を上記のように手作業で適用する方法は、データが大量にある場合、合うグラフを探すのにとても時間がかかってしまい、また目視でグラフを選ぶため良い判断基準がないという欠点がある。そのため u-w 図の作成を効率よく正確に行い、成長曲線の解析に役立てることを目的として u-w 法の計算ソフトウェアを作成する。次に実際にそれを使って年輪データの解析を行い、成長曲線を正しく計算できるか検証し、u-w 法について考察する。

4.方法
4-1 ソフトウェア開発
本ソフトウェア開発には Microsoft Visual Basic 2008 Express Edition を使った。またこれを使用して実際に成長解析を行った。Visual Basic は Windows 用のソフトウェア（GUI アプリケーション）を比較的簡単に作れるという特徴がある。 図 5 はソフトウェアのソースコードの一部である。

図 5.Microsoft Visual Basic 2008 による開発画面
3

4

4-2.プログラムの構造
Excel(CSV)ファイルのデータ読み取り

読み取り完了uとwをグラフ上にプロット

当てはめ範囲を選択して曲線を描く

曲線の選び方曲線描画

手動・・・・定数を変えながら残差平方和を参考に曲線を探す

曲線を追加する

自動描画

手動描画

モード

モード

自動・・・・定数を自動的に決めて残差平方和が最小のときの曲線を引く

正しい曲線と判断

計算結果と定数(A,c,m,W)を出力

図 6.状態遷移図

図 6 はソフトウェアの動作を示した状態遷移図である。まず Excel ファイル（CSV ファイルでも可）から年輪データを読み取り、u と w を計算してプロットする、次に曲線を選び、最後に計算で得られた結果とグラフを保存するしくみとなっている。
Excel ファイルには 1 列目に樹齢、2 列目に樹齢に対応した w があるようにする。曲線を選ぶには、手動で選ぶ方法と自動で選ぶ方法がある。手動で選ぶ方法は定数の値をグラフの変化を見ながら入力し、よく当てはまるグラフを選ぶ。
自動的に曲線を引く方法は、A,c,w を一定量ずつ変化させ残差平方和が最小のときの A,c,w を返して曲線を引くというしくみとなっている。つまり S が最小のときが最も当てはまる曲線であるというようにした。キザミ幅は表 2 のようにした。ただし m はほとんどの場合 1 に近い値をとるので、ここでは常に 1 として計算している。S は u の実測値と推定値の残差平方和であり、以下の（2）式で表す。
nS = ∑ (un − f (w)n )

2

(2)

i =1

5

u n ：表 1 で計算した uf(w)：(1)式に u n のときの w と定数を適当にあてはめた u の推定値
表 2.曲線を自動で選ぶときの定数のキザミ幅
定数

初期値

加算値

終了値

A

0.01

0.01

3

c

-1

0.05

2

W

10^( log 10 wmin )10^( log 10 wmin + i×s)s=log10 wmax − log10 wminni = 0、１、２ … n

表 2 の n は S×10 として、つまり 10^n から 10^n+1 までを等間隔で 10 分割するようにプログラムを作成した。n の値を変えれば W についてのキザミ幅を変えることができる。また u-w
図は図 2 のように複数の曲線で近似できることが多いので、複数本引けると判断した場合、別の曲線を追加できるようにもした。開発したソフトウェアのソースコードは付録に載せてある。

5.結果
5-1 開発したソフトウェア
図 7 はソフトウェアを起動した直後の状態である。[ファイル]ボタンを押すと使うファイルを選ぶ画面が表示される、選んだ Excel ファイルから樹齢と幹材積のデータを読み取り u と w が計算される。図 8 はデータファイルの例である。読み込んだデータは図 9 のようにグラフにプロットされる。
右の[Draw]ボタンを押すと曲線が引かれる、右上のテキストボックスに定数 A,c,m,W それぞれを入力するとそれに対応する曲線が引ける。[u-w 曲線]から[自動描画]を押すと自動的に近似曲線が引かれる、[曲線追加]を押すと曲線を追加できる。あてはめるデータの範囲は右のテキストボ
ックスで選択できる、何も変えなければ全ての範囲が選択される。図 10 の右下の 0.0149 は残差
平方和 S であり、この数値が小さくなるような曲線を選べばうまく近似することができる、自動描画では S の最小値を探して自動的に曲線を引くことができる。
グラフ画面の大きさはフォームの大きさを変えると自動で変わるようになっている。右下の[ズーム]ボタンを押すとグラフの一部を拡大できる、データが狭い範囲に集まっている場合はこれを使えば作業がしやすくなる。
[計算結果]を押すと図 11 のように現在計算したデータの u,w,A,c,m,W の値が Excel に出力される。また[ファイル]から[保存]を選ぶとグラフ画面も画像として保存することができる。

6

図 8.データファイルの例

図 7.ソフトウェア起動時

図 9.ファイルからデータを読み込んでプロット

図 10.自動描画機能で曲線を引いたところ

図 11.計算結果を出力

7

5-2.成長解析
ソフトウェアの検証には、2007 年に龍谷の森で伐採して年輪解析した合計 45 本の年輪データを使用した。それぞれ u-w 曲線と定数を求めた後、w の成長曲線である w-t（樹木の幹材積と時間）のグラフを作成した。

コナラ294
1.0E-01
1.0E-02w (㎥)1.0E-03
1.0E-04w-t
近似曲線

1.0E-05
1.0E-06
1.0E-07
1.0E-08
1950

1960

1970

1980

1990

2000

2010

t (year)
タカノツメ67
1.0E-03w (㎥)1.0E-04

1.0E-05w-t
近似曲線

1.0E-06

1.0E-07

1.0E-08
1995

2000

2005

2010

t(year)

図 12.コナラとタカノツメを解析した結果
左が u-w 図で右は w の成長である、このコナラ 294 では 3 段階に分けて成長曲線をあてはめた。

図 12 のコナラの場合データは u-w 関係でも w-t 関係でもほぼ満足できる程度に曲線近似された。一方の場合,u-w 関係はおおむねよく近似されたが、w-t 関係につては、充分な近似とは言い難しい。

8

5-3.R による検証
R とは統計解析を行うためのオープンソースのソフトウェアである。R の非線形最小二乗法の機能を使って u と w から各定数を導き出した。初期値には今回開発したソフトウェアで計算した
A,c,W の値を用いた。表 3 は開発したソフトウェア（初期値）と R で求めた推定値である（m の値はすべて 1 で統一してある）。
表 3,開発したソフトウェアと R による推定値の定数の比較

樹種 No.

年齢[yr]

635
391
660
660
660
35

0～6
0～10
0～23
24～46
47～57
0～13

開発したソフトウェア（初期値）
A
c
W
0.02
1.35 0.0000700
0.01
1.35 0.0000800
0.01
1.4
0.00200
0.27
10.0900
0.56
0.35
0.300
0.01
1.35 0.0000800

A
0.03300.000700
0.0201
0.711
0.205
0.406

R
c
W
1.320 0.0000531
1.537 0.000320
1.322
0.001480.847
0.05450.732
0.427
1.0600 0.0000473

今回開発したソフトウェアにより、u-w 曲線のデータのあてはめがかなり効率よく行えるようになった。このソフトウェアによって得られる係数値 A,c,W は第一近似的なものである。データへの近似の精度をあげるために、このソフトウェアで推定された A,c,W の値を初期値として、統計解析ソフト R の非線型最小二乗法を用いて、さらにあてはめを行った。なお、この場合も m=1とした。45 本のデータのうち 41 本については、今回開発したソフトウェアで推定された係数値よりも適合度の高い推定値をえることができなかった。4 本については、R の非線形最小二乗法
により表 3 に示すような推定値を得た。すなわち 45 本中 41 本については、今回開発したソフト
ウェアにより、十分な精度の A,c,W の推定値を得ることができたと言える。

9

６.考察
今回開発したソフトウェアでは、手動に頼ることなく全てのプロットデータを一度に自動で曲線近似できることが一番望ましかったが、そのようなアルゴリズムの関数は大変困難であったため実現できなかった。また定数のキザミを細かくすると計算時間が膨大な量になってしまうため、今回は表 2 で示したようにキザミ幅をある一定量変化させて定数を探す方法をとった。一部この
キザミ幅で対応できない場合もあったが、手動で定数を選ぶ方法を使えば全ての数値に対応することは可能である。よって自動的に曲線を求め、それがうまくいかなかった場合は手動であてはめるという使い方が良いと思われる。成長曲線の乗り換えについてははっきりとした基準がなくプログラムで判断するのは困難なのでユーザーの判断に任せる形となった。
今回開発したソフトウェアを使用して 2007 年に龍谷の森で伐採して年輪解析した 42 本の年輪データを解析した結果、多くの曲線はうまく近似することができたが、一部のデータではこの方法でうまく近似できないものもあった。w-t 曲線に直した場合も同様に多くの場合はうまく近似できたが、一部プロットデータのずれたものがあった。しかし u-w 図で少ない曲線であてはめた場合と多くの曲線であてはめた場合では、後者の方が w-t に直したときにきれいな近似曲線を引けることが多かった。図 13 と図 14 は同じコナラを使ってあてはめ方を変えた時のグラフである。
コナラ11
1.0E+00
1.0E-01
1.0E-02w (㎥)1.0E-03
1.0E-04w-t
近似曲線

1.0E-05
1.0E-06
1.0E-07
1.0E-08
1.0E-09
1950

1960

1970

1980

1990

2000

2010

t(year)

図 13.u-w 図で曲線を 4 段階に分けた場合
コナラ11
1.0E+00
1.0E-01
1.0E-02w (㎥)1.0E-03
1.0E-04w-t
近似曲線

1.0E-05
1.0E-06
1.0E-07
1.0E-08
1.0E-09
1950

1960

1970

1980

1990

t(year)10
図 14.u-w 図で一本の曲線で近似した場合

2000

2010

2020

定数 c は 1～1.4 の値をとることが多かったが、その範囲を大きく超えるものもまれにあった。
特に A については変動幅が大きかったので上記の方法では課題が残る形となった。A と W については、A の値が大きくなると、u-w 曲線は全体が上方向に移動する。また W の値が大きくなると、u-w 曲線は右斜め下方向に移動する特徴がある（宮浦,1990）。これをうまく利用するアルゴ
リズムを組めば多くの数値に対応でき、より正確に成長解析ができるのではないかと考えられる。
また今回計算した定数では c は 1～1.4 の範囲をとることが多かったため、多くの場合ロジステ
ィック曲線、指数関数曲線、ミッチャーリッヒ曲線などの有名な成長曲線式のいずれにも当てはまらないと言える。そのため u-w 法を使った方がより正確な成長曲線を導き出せると考えられる。
成長曲線の乗り換え時期をまとめた結果、若い樹種は一本の曲線で近似できたものが多かったが、50 年生を超えるような大きなものはほとんどの樹木に乗り換えが発生しており、90 年生を超えるアカマツでは何度も乗り換えが発生していた。また成長曲線の乗り換えが起きるごとに W が大きくなる傾向がほとんどの樹木で見られたことと、u-w 図でも多くの曲線をあてはめたほうがうまく近似できたことから、多くの樹木は成長曲線の乗り換えを何度も行いながら、つまり成長曲
線を何度も変えながら成長していくと考えられる。
図 15 はウワミズザクラとコナラの成長曲線の乗り換えが発生した年を一本ごとに示した図である。線の始まるところが樹木の発生した年である。乗り換えの発生した年をまとめた結果、全体ではどの年代でもまんべんなく乗り換えの発生があったが、コナラとウワミズザクラでは図 15のように 1970 年代に乗り換えが多く発生しているという特徴が見られた。1978 年には松枯れ病の被害が全国ピークを示しており、滋賀県では 1970 年頃から被害が増加したという報告があった（引用）。そのため 70 年代に多く発生した乗り換えは松枯れ病によるアカマツの衰退によって、他の樹種の光条件等が改善されたためだと考えられる。今後 u-w 法を用いた今回開発したソフトウェアを使ってより多くの樹種を解析することにより、乗り換えの発生状況がさらに正確に分かってくるのではないかと期待できる。

（本数）

ウワミズザクラ

コナラ（年）
図 15.コナラとウワミズザクラの
年代ごとの乗り換え発生状況
11７.引用文献
宮浦富保（1990）：新しい成長解析法「u-w 法」の紹介と応用例.北方林業 Vol.42 No5 129-135
宮浦富保（2002）：樹木の成長について.龍谷理工ジャーナル 14 巻 2 号 4-11
Kazuo Hozumi（1985）

Phase Diagrammatic Approach to the Analysis of Growth Curve

Using the u-w Diagram. Bot Mag Tokyo 98 239-250
日本生態学会（2004）：生態学入門、東京化学同人出版
松枯れ病（2009 年１月閲覧）：http://ajisime.web.fc2.com/Matugare/Matugare.htm

８.謝辞
本論文の作成にあたり、方針と助言を与えてくださった龍谷大学理工学部環境ソリューション工学科の宮浦富保教授に心より感謝します。

12

付録１
左が今回開発したソフトウェアで描画した u-w 図、右がそれを w と t の関係に直したグラフ。

アカマツ
1.0E+00
1.0E-01
1.0E-02w (㎥)1.0E-03w-t
近似曲線

1.0E-04
1.0E-05
1.0E-06
1.0E-07
1.0E-08
1900

1920

1940

1960

1980

2000

2020

t (year)アラカシ347
1.0E-02

1.0E-03w (㎥)1.0E-04w-t
近似曲線

1.0E-05

1.0E-06

1.0E-07
1994

1996

1998

2000

2002

2004

2006

2008

t(year)アラカシ379
1.0E-02
1.0E-03w (㎥)1.0E-04
1.0E-05w-t
近似曲線

1.0E-06
1.0E-07
1.0E-08
1990

1995

2000

t(year)

13

2005

2010

アラカシ522w (㎥)1.0E+01

1.0E+00w-t
近似曲線

1.0E-01
1990

1995

2000

2005

2010

t(year)アラカシ542
1.0E-02

1.0E-03w (㎥)1.0E-04w-t
近似曲線

1.0E-05

1.0E-06

1.0E-07
1994

1996

1998

2000

2002

2004

2006

2008

t(year)アラカシ635
1.0E-02
1.0E-03w (㎥)1.0E-04
1.0E-05w-t
近似曲線

1.0E-06
1.0E-07
1.0E-08
2001

2002

2003

2004

2005

t(year)

14

2006

2007

2008

ウワミズザクラ107
1.0E-01

1.0E-02w (㎥)1.0E-03w-t
近似曲線

1.0E-04

1.0E-05

1.0E-06
1970

1980

1990

2000

2010

t (year)

ウワミズザクラ116
1.0E+00
1.0E-01w (㎥)1.0E-02
1.0E-03w-t
近似曲線

1.0E-04
1.0E-05
1.0E-06
1.0E-07
1970

1980

1990

2000

2010

t(year)

ウワミズザクラ858
1.0E+00
1.0E-01w (㎥)1.0E-02
1.0E-03w-t
近似曲線

1.0E-04
1.0E-05
1.0E-06
1.0E-07
1980

1985

1990

1995

t(year)

15

2000

2005

2010

ウワミズザクラ750
1.0E+00
1.0E-01w (㎥)1.0E-02
1.0E-03w-t
近似曲線

1.0E-04
1.0E-05
1.0E-06
1950

1960

1970

1980

1990

2000

2010

t(year)

カナメモチ391
1.0E-02
1.0E-03w (㎥)1.0E-04
1.0E-05w-t
近似曲線

1.0E-06
1.0E-07
1.0E-08
1996

1998

2000

2002

2004

2006

2008

t(year)

カナメモチ857
1.0E-03w (㎥)1.0E-04

1.0E-05w-t
近似曲線

1.0E-06

1.0E-07
1990

1995

2000

t(year)

16

2005

2010

コシアブラ90
1.0E-01
1.0E-02w (㎥)1.0E-03
1.0E-04w-t
近似曲線

1.0E-05
1.0E-06
1.0E-07
1.0E-08
1990

1995

2000

2005

2010

t(year)コシアブラ118
1.0E-02

1.0E-03w (㎥)1.0E-04w-t
近似曲線

1.0E-05

1.0E-06

1.0E-07
1996

1998

2000

2002

2004

2006

2008

t(year)コシアブラ570
1.0E-02

1.0E-03w (㎥)1.0E-04w-t
近似曲線

1.0E-05

1.0E-06

1.0E-07
1998

2000

2002

2004

t(year)

17

2006

2008

コシアブラ804
1.0E-01
1.0E-02w (㎥)1.0E-03
1.0E-04w-t
近似曲線

1.0E-05
1.0E-06
1.0E-07
1.0E-08
1985

1990

1995

2000

2005

2010

t(year)コシアブラ519
1.0E-02

1.0E-03w (㎥)1.0E-04w-t
近似曲線

1.0E-05

1.0E-06

1.0E-07
1998

2000

2002

2004

2006

2008

t(year)

コナラ11
1.0E+00
1.0E-01
1.0E-02w (㎥)1.0E-03
1.0E-04w-t
近似曲線

1.0E-05
1.0E-06
1.0E-07
1.0E-08
1.0E-09
1950

1960

1970

1980

1990

t(year)

18

2000

2010

2020

コナラ22
1.0E+00
1.0E-01w (㎥)1.0E-02
1.0E-03w-t
近似曲線

1.0E-04
1.0E-05
1.0E-06
1950

1960

1970

1980

1990

2000

2010

t(year)

コナラ294
1.0E-01
1.0E-02w (㎥)1.0E-03
1.0E-04w-t
近似曲線

1.0E-05
1.0E-06
1.0E-07
1.0E-08
1950

1960

1970

1980

1990

2000

2010

t (year)

コナラ343
1.0E+00
1.0E-01w (㎥)1.0E-02
1.0E-03w-t
近似曲線

1.0E-04
1.0E-05
1.0E-06
1.0E-07
1950

1960

1970

1980

t(year)

19

1990

2000

2010

コナラ590
1.0E+00
1.0E-01w (㎥)1.0E-02
1.0E-03w-t
近似曲線

1.0E-04
1.0E-05
1.0E-06
1960

1970

1980

1990

2000

2010

t(year)

コナラ6601.0E+00
1.0E-01w (㎥)1.0E-02
1.0E-03w-t
近似曲線

1.0E-04
1.0E-05
1.0E-06
1.0E-07
1940

1950

1960

1970

1980

1990

2000

2010

t(year)コバノミツバツツジ82
1.0E-02

1.0E-03w (㎥)1.0E-04w-t
近似曲線

1.0E-05

1.0E-06

1.0E-07
1990

1995

2000

t(year)

20

2005

2010

コバノミツバツツジ797
1.E-03w (㎥)

1.E-04

1.E-05w-t近似曲線
1.E-06

1.E-07
1980

1985

1990

1995

2000

2005

2010

t (year)

コバノミツバツツジ576_1
1.0E-02

1.0E-03w (㎥)1.0E-04w-t
近似曲線

1.0E-05

1.0E-06

1.0E-07
1985

1990

1995

2000

2005

2010

t(year)コバノミツバツツジ576_2
1.0E-02

1.0E-03w (㎥)1.0E-04w-t
近似曲線

1.0E-05

1.0E-06

1.0E-07
1985

1990

1995

2000

t(year)

21

2005

2010

ソヨゴ532
1.0E-01

1.0E-02w (㎥)1.0E-03

1.0E-04w-t
近似曲線

1.0E-05

1.0E-06

1.0E-07
1950

1960

1970

1980

1990

2000

2010

t (year)

タカノツメ35
1.0E-02

1.0E-03w (㎥)1.0E-04w-t
近似曲線

1.0E-05

1.0E-06

1.0E-07
1990

1995

2000

2005

2010

t(year)

タカノツメ67
1.0E+00
1.0E-01
1.0E-02w (㎥)1.0E-03
1.0E-04w-t
近似曲線

1.0E-05
1.0E-06
1.0E-07
1.0E-08
1994

1996

1998

2000

2002

t(year)

22

2004

2006

2008

2010

タカノツメ223
1.0E-01
1.0E-02w (㎥)1.0E-03
1.0E-04w-t
近似曲線

1.0E-05
1.0E-06
1.0E-07
1994

1996

1998

2000

2002

2004

2006

2008

2010

t(year)

タカノツメ257
1.0E-01
1.0E-02w (㎥)1.0E-03
1.0E-04w-t
近似曲線

1.0E-05
1.0E-06
1.0E-07
1.0E-08
1990

1995

2000

2005

2010

t(year)

タカノツメ205
1.0E-02

1.0E-03w (㎥)1.0E-04w-t
近似曲線

1.0E-05

1.0E-06

1.0E-07
1998

2000

2002

t(year)

23

2004

2006

ヒサカキ4
1.0E-02

1.0E-03w (㎥)1.0E-04w-t
近似曲線

1.0E-05

1.0E-06

1.0E-07
1990

1995

2000

2005

2010

t(year)

ヒサカキ757
1.0E-01
1.0E-02w (㎥)1.0E-03
1.0E-04w-t
近似曲線

1.0E-05
1.0E-06
1.0E-07
1.0E-08
1960

1970

1980

1990

2000

2010

t(year)

ヤマザクラ71
1.0E-01
1.0E-02w (㎥)1.0E-03
1.0E-04w-t
近似曲線

1.0E-05
1.0E-06
1.0E-07
1970

1980

1990

t(year)

24

2000

2010

ヤマザクラ148
1.0E+00
1.0E-01
1.0E-02w (㎥)1.0E-03
1.0E-04w-t
近似曲線

1.0E-05
1.0E-06
1.0E-07
1.0E-08
1960

1970

1980

1990

2000

2010

2020

t(year)

ヤマザクラ305
1.0E-03

1.0E-04w (㎥)1.0E-05w-t
近似曲線

1.0E-06

1.0E-07

1.0E-08
1985

1990

1995

2000

2005

2010

t (year)

ヤマザクラ523
1.0E+00
1.0E-01w (㎥)1.0E-02
1.0E-03w-t
近似曲線

1.0E-04
1.0E-05
1.0E-06
1950

1960

1970

1980

t(year)

25

1990

2000

2010

リョウブ392
1.0E-01
1.0E-02w (㎥)1.0E-03
1.0E-04w-t
近似曲線

1.0E-05
1.0E-06
1.0E-07
1.0E-08
1996

1998

2000

2002

2004

2006

2008

t(year)

リョウブ423
1.0E-02w (㎥)1.0E-03

1.0E-04w-t
近似曲線

1.0E-05

1.0E-06
1994

1996

1998

2000

2002

2004

2006

2008

t(year)

リョウブ470
1.0E-02w (㎥)1.0E-03

1.0E-04w-t
近似曲線

1.0E-05

1.0E-06

1.0E-07
1996

1998

2000

2002

t(year)

26

2004

2006

2008

ウワミズザクラ141
1.0E+00
1.0E-01w (㎥)1.0E-02
1.0E-03w-t
近似曲線

1.0E-04
1.0E-05
1.0E-06
1.0E-07
1970

1980

1990

2000

2010

t(year)

ウワミズザクラ515
1.0E+00
1.0E-01w (㎥)1.0E-02
1.0E-03w-t
近似曲線

1.0E-04
1.0E-05
1.0E-06
1.0E-07
1960

1970

1980

1990

2000

2010

t(year)

コナラ231
1.0E+00
1.0E-01w (㎥)1.0E-02
1.0E-03w-t
近似曲線

1.0E-04
1.0E-05
1.0E-06
1.0E-07
1950

1960

1970

1980

t(year)

27

1990

2000

2010

付録 2
以下は開発したソフトウェアのソースコードである。
'グラフスケールの基礎設定
Public Class Graph
Protected xOLd, yOLd As Single
Protected kx, ky As Single
Protected MyPen As New Pen(Color.Black, 1)
Protected MyBrush As New SolidBrush(Color.Black)Protected g As Graphics
Protected bmp As Bitmap'コンストラクタ
Public Sub New(ByRef pic As System.Windows.Forms.PictureBox, ByVal xmin As Single,
ByVal xmax As Single, ByVal ymin As Single, ByVal ymax As Single, ByVal dx As Single,
ByVal dy As Single)
pic.Image = New Bitmap(pic.Width, pic.Height)
g = Graphics.FromImage(pic.Image)g.FillRectangle(Brushes.White, 0, 0, pic.Width, pic.Height) '背景色を白に
kx = (pic.Width - 2 * dx) / (xmax - xmin)
ky = (pic.Height - 2 * dy) / (ymax - ymin)g.TranslateTransform(dx - xmin * kx, pic.Height - dy + ymin * ky)g.ScaleTransform(1, -1)
End Sub'グラフの主軸
Public Sub Axis(ByVal x1 As Double, ByVal x2 As Double, ByVal y1 As Double, ByVal y2
As Double, ByVal cL As Color, ByVal index As Integer)
Select Case index
Case 1
MyPen.Color = cL
g.DrawLine(MyPen, CInt(kx * x1), CInt(ky * 0), CInt(kx * x2), CInt(ky *y1))g.DrawLine(MyPen, CInt(kx * x1), CInt(ky * y1), CInt(kx * x1), CInt(ky *y2))Case 2
MyPen.Color = cL
g.DrawLine(MyPen, CInt(kx * x1), CInt(ky * y1), CInt(kx * x2), CInt(ky *
28

y1))
g.DrawLine(MyPen, CInt(kx * x1), CInt(ky * y2), CInt(kx * x1), CInt(ky *y1))Case 3
MyPen.Color = cL
g.DrawLine(MyPen, CInt(kx * x1), CInt(ky * y2), CInt(kx * x2), CInt(ky *y2))g.DrawLine(MyPen, CInt(kx * x1), CInt(ky * y2), CInt(kx * x1), CInt(ky *y1))Case 4
MyPen.Color = cL
g.DrawLine(MyPen, CInt(kx * x1), CInt(ky * 0), CInt(kx * x2), CInt(ky * 0))
g.DrawLine(MyPen, CInt(kx * 0), CInt(ky * y1), CInt(kx * 0), CInt(ky * y2))
End SelectEnd Sub'線を引く
Public Sub Line(ByVal x1 As Double, ByVal y1 As Double, ByVal x2 As Double, ByVal y2
As Double, ByVal cL As Color, ByVal index As Integer)
MyPen.Color = cL
Select Case index
Case 1
MyPen.DashStyle = Drawing2D.DashStyle.Solidg.DrawLine(MyPen, CInt(kx * x1), CInt(ky * y1), CInt(kx * x2), CInt(ky *y2))Case 2
MyPen.DashStyle = Drawing.Drawing2D.DashStyle.Dash
g.DrawLine(MyPen, CInt(kx * x1), CInt(ky * y1), CInt(kx * x2), CInt(ky *y2))
End SelectxOLd = x2 : yOLd = x2
End Sub'続けて線を引く
Public Sub LineTO(ByVal x As Double, ByRef y As Double, ByVal cL As Color)
MyPen.Color = cL
g.DrawLine(MyPen, CInt(kx * xOLd), CInt(ky * yOLd), CInt(kx * x), CInt(ky * y))xOLd = x : yOLd = y
End Sub
Public Sub PutPoint(ByVal x As Double, ByVal y As Double, ByVal cL As Color)
MyBrush.Color = cL
g.FillRectangle(MyBrush, CInt(kx * x), CInt(ky * y), 1, 1)
29

xOLd = x : yOLd = y
End Sub'テキストの初期設定
Public Sub Drawtext(ByVal s As String, ByVal x As Double, ByVal y As Double, ByVal cL
As Color, ByVal fsize As Integer)
Dim myFont As New Font("MS 明朝", fsize, FontStyle.Regular)
MyBrush.Color = cL
g.ScaleTransform(1, -1)
g.DrawString(s, myFont, MyBrush, kx * x, -ky * y)g.ScaleTransform(1, -1)
End Sub
Public Sub Drawtext_y_Axis1(ByVal s As String, ByVal x As Double, ByVal y As Double,
ByVal cL As Color, ByVal fsize As Integer)
Dim myFont As New Font("MS 明朝", fsize, FontStyle.Regular)
MyBrush.Color = cL
g.ScaleTransform(1, -1)
g.DrawString(s, myFont, MyBrush, (kx * x) - 35, -ky * y)g.ScaleTransform(1, -1)
End Sub
Public Sub Drawtext_y_Axis2(ByVal s As String, ByVal x As Double, ByVal y As Double,
ByVal cL As Color, ByVal fsize As Integer)
Dim myFont As New Font("MS 明朝", fsize, FontStyle.Regular)
MyBrush.Color = cL
g.ScaleTransform(1, -1)
g.DrawString(s, myFont, MyBrush, (kx * x) - 50, -ky * y)g.ScaleTransform(1, -1)
End Sub
Public Sub Drawtext_x_Axis1(ByVal s As String, ByVal x As Double, ByVal y As Double,
ByVal cL As Color, ByVal fsize As Integer)
Dim myFont As New Font("MS 明朝", fsize, FontStyle.Regular)
MyBrush.Color = cL
g.ScaleTransform(1, -1)
g.DrawString(s, myFont, MyBrush, kx * x, (-ky * y) + 20)
g.ScaleTransform(1, -1)
End Sub'プロットのマーカー設定のプロージャー
Public Sub Plot(ByVal x As Double, ByVal y As Double, ByVal cL As Color, ByVal size As
Integer, ByVal index As Integer)
Try
MyBrush.Color = cL
MyPen.Color = cL
30

Select Case index
Case 1 '×記号
g.DrawLine(MyPen, CInt(kx * x) - size, CInt(ky * y) - size, CInt(kx * x)
+ size, CInt(ky * y) + size)g.DrawLine(MyPen, CInt(kx * x) - size, CInt(ky * y) + size, CInt(kx * x)
+ size, CInt(ky * y) - size)
Case 2 '○記号
g.DrawEllipse(MyPen, CInt(kx * x) - size, CInt(ky * y) - size, 2 * size, 2* size)Case 3 '●記号
g.FillEllipse(MyBrush, CInt(kx * x) - size, CInt(ky * y) - size, 2 * size, 2* size)End SelectCatch
MsgBox("計算が正しくありません。")
End Try
End Sub
End Class'フォーム
Public Class MainForm
Private g1 As Graph
Const dx As Integer = 55 '横の余白Const dy As Integer = 45 '縦の余白
Public y_min As Double = 0.001 'ｘの最小値
Public y_max As Double = 1000 'ｘの最大値
Public x_min As Double = 0.001 'ｙの最小値
Public x_max As Double = 100 'ｙの最大値
Private log_y_min As Double = Math.Log10(y_min)
Private log_y_max As Double = Math.Log10(y_max)
Private log_x_min As Double = Math.Log10(x_min)
Private log_x_max As Double = Math.Log10(x_max)

31

Public W_index(200) As Double 'W の値を格納する配列
Public U_index(200) As Double 'U の値を格納する配列
Public Plot As Integer 'プロットの数
Public k1, k2 As Integer 'プロットの範囲
Public line As Integer = 0 '曲線の数をカウント
Public consts(20, 3) As Double '曲線ごとの定数を格納(曲線の数の上限、4 つの定数)
Private k(10) As Integer 'プロットの範囲を記憶
Private ret As DialogResult 'ダイアログの結果を格納する変数'初期化
PrivateSub

MainForm_Load(ByVal

senderAs

System.Object,

ByVal

e

As

System.EventArgs) Handles MyBase.Load
log_y_min = Math.Log10(y_min)log_y_max = Math.Log10(y_max)
log_x_min = Math.Log10(x_min)log_x_max = Math.Log10(x_max)
g1 = New Graph(pcbox1, log_x_min, log_x_max, log_y_min, log_y_max, dx, dy)
GraphBase(g1, x_min, x_max, y_min, y_max)
End Sub'Excel ファイルの場合<System.Runtime.InteropServices.ComVisible(True)> Private Sub Excel フ ァ イ ル
EToolStripMenuItem_Click(ByVal sender As System.Object, ByVal e As System.EventArgs)
Handles Excel ファイル EToolStripMenuItem.Click'Excel 関連の変数
Dim xlsSheet As Object
Dim xls As Object
Dim xlsBook As Object
Dim x_max2, x_min2, y_max2, y_min2 As Double'ダイアログの設定
OpenFileDialog1.FileName = ""OpenFileDialog1.Filter = "Excel ファイル(*.xls)|*.xls"
OpenFileDialog1.InitialDirectory = "C:¥temp"

32

ret = OpenFileDialog1.ShowDialog()
Try'ダイアログを開いたら以下の処理
If ret = Windows.Forms.DialogResult.OK Thenxls = CreateObject("Excel.Application")
xlsBook = xls.Workbooks.Open(OpenFileDialog1.FileName)
xlsSheet = xlsBook.sheets("Sheet1")'UW の計算
figure(xlsSheet, W_index, U_index, Plot)x_max2 = W_index(Plot - 1)
x_min2 = W_index(0)y_max2 = U_index(0)y_min2 = U_index(Plot - 1)
'グラフメモリの再計算
max_min(x_min, x_max, y_min, y_max, x_min2, x_max2, y_min2, y_max2)'グラフの書き込み()
MainForm_Load(sender, e)'U と W をプロット
PlotPoint(g1, W_index, U_index, Plot)Label12.Text = Plot
'Excel の終了
xls.Quit()
xlsSheet = Nothingxls = Nothing
xlsBook = Nothing
End IfCatch
MsgBox("Sheet1 に正しい内容を入力してください")
End Try
End Sub'CSV ファイルの場合
Private Sub CSV ファイル CToolStripMenuItem_Click(ByVal sender As System.Object,
ByVal e As System.EventArgs) Handles CSV ファイル CToolStripMenuItem.Click'CSV 読み込み33

Dim sr As System.IO.StreamReader
Dim x_max2, x_min2, y_max2, y_min2 As Double
OpenFileDialog1.FileName = ""OpenFileDialog1.Filter = "CSV ファイル(*.csv)|*.csv"
OpenFileDialog1.InitialDirectory = "C:¥temp"
ret = OpenFileDialog1.ShowDialog()
Try'ダイアログを開いたら以下の処理
If ret = Windows.Forms.DialogResult.OK Then
'CSV ファイルの読み込み
sr

=

New

System.IO.StreamReader(OpenFileDialog1.FileName,

System.Text.Encoding.Default)'計算
figure(sr, W_index, U_index, Plot)x_max2 = W_index(Plot - 1)
x_min2 = W_index(0)y_max2 = U_index(0)y_min2 = U_index(Plot - 1)
'グラフメモリの再計算
max_min(x_min, x_max, y_min, y_max, x_min2, x_max2, y_min2, y_max2)'グラフの書き込み
MainForm_Load(sender, e)'U と W をプロット
PlotPoint(g1, W_index, U_index, Plot)Label12.Text = Plot
sr.Close()
End IfCatch ex As Exception
MsgBox("このファイルの内容では計算できません。")
End Try
End Sub'フォームのサイズを変更したとき
Private Sub frmGraphsample1_SizeChanged(ByVal sender As Object, ByVal e As
34

System.EventArgs) Handles Me.SizeChanged
Dim A, c, W, m As Double
g1 = New Graph(pcbox1, log_x_min, log_x_max, log_y_min, log_y_max, dx, dy)'メモリの呼び出し
GraphBase(g1, x_min, x_max, y_min, y_max)'ダイアログを開いていない場合ここで終了
If ret = Windows.Forms.DialogResult.OK Then
'U と W をプロット
PlotPoint(g1, W_index, U_index, Plot)
End If
For i = 0 To line
A = consts(i, 0)c = consts(i, 1)
W = consts(i, 2)m = consts(i, 3)
curve(g1, A, c, W, m, x_min, x_max)
Next
End Sub'曲線の描画
PrivateSub

Button1_Click(ByVal

senderAs

System.Object,

ByVal

e

System.EventArgs) Handles Button1.Click
Dim A, c, W, m As Double
consts(line, 0) = Val(TextBox3.Text)
consts(line, 1) = Val(TextBox2.Text)
consts(line, 2) = Val(TextBox4.Text)
consts(line, 3) = Val(TextBox1.Text)
A = consts(line, 0)c = consts(line, 1)
W = consts(line, 2)m = consts(line, 3)

g1 = New Graph(pcbox1, log_x_min, log_x_max, log_y_min, log_y_max, dx, dy)
GraphBase(g1, x_min, x_max, y_min, y_max)

35

As

'ファイルを開いているときだけプロットと誤差の計算-------------------------If ret = Windows.Forms.DialogResult.OK Then
k1 = Val(TextBox5.Text)
k2 = Val(TextBox6.Text)'誤差を計算する範囲を設定する処理
If k1 = 0 Or k2 < k1 Or k1 > Plot Then
k1 = 1End IfIf k2 > Plot Or k2 = 0 Or k2 > Plot Then
k2 = Plot
End IfTextBox5.Text = k1
TextBox6.Text = k2
'U と W をプロット
PlotPoint(g1, W_index, U_index, Plot, k1, k2)'誤差の表示Label7.Text = Format(LSM(A, c, W, m, W_index, U_index, k2, k1), "0.0000")k(line) = k2
End If '--------------------------------------------------------------------'line の数だけ曲線を描く
For i = 0 To line
A = consts(i, 0)c = consts(i, 1)
W = consts(i, 2)m = consts(i, 3)
curve(g1, A, c, W, m, x_min, x_max)
g1.Drawtext(i + 1 & " 段 階 ", Math.Log10(x_min) + (i * 0.8) + 0.2,
Math.Log10(y_min) + 1.2, Color.Black, 9)
g1.Drawtext("m

="

&

m,

Math.Log10(x_min)

+

(i

*

0.8)

+

0.2,

Math.Log10(y_min) + 1.0, Color.Black, 9)
g1.Drawtext("c =" & c, Math.Log10(x_min) + (i * 0.8) + 0.2, Math.Log10(y_min)+ 0.8, Color.Black, 9)
g1.Drawtext("A =" & A, Math.Log10(x_min) + (i * 0.8) + 0.2, Math.Log10(y_min)+ 0.6, Color.Black, 9)36

g1.Drawtext("W

="

&

W,

Math.Log10(x_min)

+

(i

*

0.8)

+

0.2,

Math.Log10(y_min) + 0.4, Color.Black, 9)
Next
End Sub'曲線を自動的に引く、計算速度遅Private Sub 自動描写 ToolStripMenuItem_Click(ByVal sender As System.Object, ByVal e
As System.EventArgs) Handles 自動描写 ToolStripMenuItem.Click
Dim A, c, W, m, s1, s2 As Double
Dim min As Double = 1000
Dim A_sub, c_sub, W_sub As Double
Dim U As Double = 0
Dim U2 As Double = 0
m=1
s1 = Math.Log10(x_max) - Math.Log10(x_min)s2 = s1 / (s1 * 10)consts(line, 0) = Val(TextBox3.Text)
consts(line, 1) = Val(TextBox2.Text)
consts(line, 2) = Val(TextBox4.Text)
consts(line, 3) = Val(TextBox1.Text)
g1 = New Graph(pcbox1, log_x_min, log_x_max, log_y_min, log_y_max, dx, dy)
GraphBase(g1, x_min, x_max, y_min, y_max)'Excel ファイルを開いているときだけプロットの計算
If ret = Windows.Forms.DialogResult.OK Then
k1 = Val(TextBox5.Text)
k2 = Val(TextBox6.Text)
If k1 = 0 Or k2 < k1 Or k1 > Plot Then
k1 = 1End IfIf k2 > Plot Or k2 = 0 Or k2 > Plot Then
k2 = Plot
End IfTextBox5.Text = k1
TextBox6.Text = k2
'U と W をプロット
37

PlotPoint(g1, W_index, U_index, Plot, k1, k2)'一番誤差が低いときの定数を返す処理'定数を変化させるループ
For N As Integer = 0 To s1 * 10
W_sub = 10 ^ (Math.Log10(x_min) + (N * s2))
For A_sub = 0.01 To 4 Step 0.01
For c_sub = -1 To 2 Step 0.05
'最小二乗法の計算()
U = LSM(A_sub, c_sub, W_sub, m, W_index, U_index, k2, k1)' 最小値を返す()
If U < min Thenmin = UA = A_subc = c_sub
W = W_sub
End If
Next
Next
Next
W = keta(W, 3, 0)
consts(line, 0) = keta(A, 3, 0)
consts(line, 1) = keta(c, 3, 0)
consts(line, 2) = keta(W, 3, 0)
consts(line, 3) = keta(m, 3, 0)
TextBox3.Text = consts(line, 0)
TextBox2.Text = consts(line, 1)
TextBox4.Text = consts(line, 2)
TextBox1.Text = consts(line, 3)'誤差の表示Label7.Text = Format(min, "0.0000")k(line) = k2
End If'line の数だけ曲線を描く
For i = 0 To line38A = consts(i, 0)c = consts(i, 1)
W = consts(i, 2)m = consts(i, 3)
curve(g1, A, c, W, m, x_min, x_max)
g1.Drawtext(i + 1 & "段階", Math.Log10(x_min) + (i * 0.8) + 0.2,
Math.Log10(y_min) + 1.2, Color.Black, 9)
g1.Drawtext("m =" & m, Math.Log10(x_min) + (i * 0.8) + 0.2,
Math.Log10(y_min) + 1.0, Color.Black, 9)
g1.Drawtext("c =" & c, Math.Log10(x_min) + (i * 0.8) + 0.2,
Math.Log10(y_min) + 0.8, Color.Black, 9)
g1.Drawtext("A =" & A, Math.Log10(x_min) + (i * 0.8) + 0.2,
Math.Log10(y_min) + 0.6, Color.Black, 9)
g1.Drawtext("W =" & W, Math.Log10(x_min) + (i * 0.8) + 0.2,
Math.Log10(y_min) + 0.4, Color.Black, 9)
Next
End Sub'終了
Private Sub 終了 ToolStripMenuItem_Click(ByVal sender As System.Object, ByVal e As
System.EventArgs) Handles 終了 ToolStripMenuItem.Click
Me.Close()End
End Sub'曲線のやり直し
Private Sub 消去 CToolStripMenuItem_Click(ByVal sender As System.Object, ByVal e As
System.EventArgs) Handles 消去 CToolStripMenuItem.Click
pcbox1.Refresh()
g1 = New Graph(pcbox1, log_x_min, log_x_max, log_y_min, log_y_max, dx, dy)
GraphBase(g1, x_min, x_max, y_min, y_max)line = 0Label10.Text = 1
Erase kReDim k(10)If ret = Windows.Forms.DialogResult.OK Then
'U と W をプロット
39

PlotPoint(g1, W_index, U_index, Plot)
End IfEnd Sub'曲線を追加Private Sub 追加描写 AToolStripMenuItem_Click(ByVal sender As System.Object, ByVale As System.EventArgs) Handles 追加描写 AToolStripMenuItem.Click
Dim A, c, W, m As Double

g1 = New Graph(pcbox1, log_x_min, log_x_max, log_y_min, log_y_max, dx, dy)
GraphBase(g1, x_min, x_max, y_min, y_max)'Excel ファイルを開いているときだけプロットの計算
If ret = Windows.Forms.DialogResult.OK Then
k1 = Val(TextBox5.Text)
k2 = Val(TextBox6.Text)
If k1 = 0 Or k2 < k1 Or k2 > Plot Then
k1 = 1End IfIf k2 > Plot Or k2 = 0 Or k2 > Plot Then
k2 = Plot
End IfTextBox5.Text = k1
TextBox6.Text = k2
'U と W をプロット
PlotPoint(g1, W_index, U_index, Plot, k1, k2)'誤差の表示Label7.Text = Format(LSM(A, c, W, m, W_index, U_index, k2, k1), "0.0000")k(line) = k2
End If'曲線を一本ずつ増やす処理-----------------------------------------------------------consts(line, 0) = Val(TextBox3.Text)
consts(line, 1) = Val(TextBox2.Text)
consts(line, 2) = Val(TextBox4.Text) '上限値,これより増えるとエラー
consts(line, 3) = Val(TextBox1.Text)40

'line の数だけ曲線を描く
For i = 0 To line
A = consts(i, 0)c = consts(i, 1)
W = consts(i, 2)m = consts(i, 3)
curve(g1, A, c, W, m, x_min, x_max)
g1.Drawtext(i + 1 & " 段 階 ", Math.Log10(x_min) + (i * 0.8) + 0.2,
Math.Log10(y_min) + 1.2, Color.Black, 9)
g1.Drawtext("m

="

&

m,

Math.Log10(x_min)

+

(i

*

0.8)

+

0.2,

Math.Log10(y_min) + 1.0, Color.Black, 9)
g1.Drawtext("c =" & c, Math.Log10(x_min) + (i * 0.8) + 0.2, Math.Log10(y_min)+ 0.8, Color.Black, 9)
g1.Drawtext("A =" & A, Math.Log10(x_min) + (i * 0.8) + 0.2, Math.Log10(y_min)+ 0.6, Color.Black, 9)
g1.Drawtext("W

="

&

W,

Math.Log10(x_min)

+

(i

*

0.8)

+

0.2,

Math.Log10(y_min) + 0.4, Color.Black, 9)
Next'曲線の数を増やす
line += 1
Label10.Text = line + 1
End Sub'計算結果
Private Sub 計算結果 FToolStripMenuItem_Click(ByVal sender As System.Object, ByVale As System.EventArgs) Handles 計算結果 FToolStripMenuItem.Click
Dim xlsSheet As Object
Dim xls As Object
Dim xlsBook As Object
xls = CreateObject("Excel.Application")
xlsBook = xls.workbooks.add
xlsSheet = xlsBook.sheets("Sheet1")
xlsBook.Application.visible = True
xlsSheet.cells(1, 2) = "w"
xlsSheet.cells(1, 3) = "u"xlsSheet.cells(1, 4) = "A"41

xlsSheet.cells(1, 5) = "c"xlsSheet.cells(1, 6) = "W"
xlsSheet.cells(1, 7) = "m"'U と W の計算結果
For i = 1 To Plot
xlsSheet.cells(i + 1, 1) = i
xlsSheet.cells(i + 1, 2) = W_index(i - 1)
xlsSheet.cells(i + 1, 3) = U_index(i - 1)
Next
'定数の計算結果
For j = 0 To line
xlsSheet.cells(k(j) + 1, 4) = consts(j, 0)
xlsSheet.cells(k(j) + 1, 5) = consts(j, 1)
xlsSheet.cells(k(j) + 1, 6) = consts(j, 2)
xlsSheet.cells(k(j) + 1, 7) = consts(j, 3)
Next
xls.Quit()
xlsSheet = Nothingxls = Nothing
xlsBook = Nothing
End Sub'グラフを画像として保存
Private Sub グラフの保存 SToolStripMenuItem_Click(ByVal sender As System.Object,
ByVal e As System.EventArgs) Handles グラフの保存 SToolStripMenuItem.Click
Dim ext As String
SaveFileDialog1.FileName = "無題.gif"
SaveFileDialog1.Filter = "GIF ファイル(*.gif)|*.gif"
If SaveFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Thenext = System.IO.Path.GetExtension(SaveFileDialog1.FileName).ToUpper()
Select Case ext
Case ".BMP", ".bmp"
pcbox1.Image.Save(SaveFileDialog1.FileName,
System.Drawing.Imaging.ImageFormat.Bmp)
Case ".jpg", ".JPG"
pcbox1.Image.Save(SaveFileDialog1.FileName,
System.Drawing.Imaging.ImageFormat.Jpeg)

42

Case ".gif", ".GIF"
pcbox1.Image.Save(SaveFileDialog1.FileName,
System.Drawing.Imaging.ImageFormat.Gif)
End SelectEnd IfEnd Sub'拡大
PrivateSub

Button2_Click(ByVal

senderAs

System.Object,

ByVal

e

As

System.EventArgs) Handles Button2.Click
Dim x_max2, x_min2, y_max2, y_min2 As Double
Dim z1, z2 As Integer
If ret = Windows.Forms.DialogResult.OK Then'UW の計算z1 = Val(TextBox7.Text)
z2 = Val(TextBox8.Text)
x_max2 = W_index(z2 - 1)
x_min2 = W_index(z1)y_max2 = U_index(z1)y_min2 = U_index(z2 - 1)
'グラフメモリの再計算
max_min(x_min, x_max, y_min, y_max, x_min2, x_max2, y_min2, y_max2)'グラフの書き込み()
MainForm_Load(sender, e)'U と W をプロット
ZoomPoint(g1, W_index, U_index, Plot, z1, z2)
End IfEnd Sub'縮小
PrivateSub

Button3_Click(ByVal

senderAs

System.EventArgs) Handles Button3.Click
Dim x_max2, x_min2, y_max2, y_min2 As Double
If ret = Windows.Forms.DialogResult.OK Thenx_max2 = W_index(Plot - 1)
x_min2 = W_index(0)y_max2 = U_index(0)43System.Object,

ByVal

e

As

y_min2 = U_index(Plot - 1)
'グラフメモリの再計算
max_min(x_min, x_max, y_min, y_max, x_min2, x_max2, y_min2, y_max2)'グラフの書き込み()
MainForm_Load(sender, e)'U と W をプロット
PlotPoint(g1, W_index, U_index, Plot)
End IfEnd Sub
End Class'MainForm で使う関数
Module figure1'グラフメモリ
Public Sub GraphBase(ByVal g1 As Graph, ByVal x_min As Double, ByVal x_max As
Double, ByVal y_min As Double, ByVal y_max As Double)
Try'グラフ軸の設定
g1.Axis(Math.Log10(x_min), Math.Log10(x_max) + 0.1, Math.Log10(y_min),Math.Log10(y_max) + 0.1, Color.Black, 2)'X 軸のメモリ(対数)'主軸
For i = Math.Log10(x_min) + 1 To Math.Log10(x_max)
g1.Line(i, Math.Log10(y_min), i, Math.Log10(y_max) + 0.1, Color.Gray, 2)
Next
For j = 0 To Math.Log10(x_max) - Math.Log10(x_min) - 1
For i As Double = x_min * 10 ^ j To x_min * 10 ^ (j + 1) Step x_min * 10 ^ j
g1.Line(Math.Log10(i),Math.Log10(y_min),Math.Log10(i),Math.Log10(y_min) + 0.1, Color.Black, 1)
Next
Next
'Y 軸のメモリ(対数)'主軸
For i = Math.Log10(y_min) + 1 To Math.Log10(y_max)
g1.Line(Math.Log10(x_min), i, Math.Log10(x_max) + 0.05, i, Color.Gray, 2)
Next
For j = 0 To Math.Log10(y_max) - Math.Log10(y_min) - 1
44For i As Double = y_min * 10 ^ j To y_min * 10 ^ (j + 1) Step y_min * 10 ^ j
g1.Line(Math.Log10(x_min), Math.Log10(i), Math.Log10(x_min) + 0.05,
Math.Log10(i), Color.Black, 1)
Next
Next
'X 軸メモリの数値
For i = Math.Log10(x_min) To Math.Log10(x_max) Step 1
g1.Drawtext("10^" & i, Math.Log10(10 ^ i), Math.Log10(y_min) - 0.1,
Color.Black, 10)
Next'Y 軸メモリの数値
For i = Math.Log10(y_min) To Math.Log10(y_max) Step 1
g1.Drawtext_y_Axis1("10^" & i, Math.Log10(x_min), Math.Log10(10 ^ i),
Color.Black, 8)
Next
g1.Drawtext_x_Axis1("w", (Math.Log10(x_max) + Math.Log10(x_min)) / 2,
Math.Log10(y_min), Color.Black, 14)
g1.Drawtext_y_Axis2("U",Math.Log10(x_min),(Math.Log10(y_max)+

Math.Log10(y_min)) / 2, Color.Black, 14)
Catch ex As Exception
MsgBox(ex.Message)
Exit Sub
End Try
End Sub'Excel のセルの値から W と U の値を計算
Public Sub figure(ByVal Sheet As Object, ByRef W_index() As Double, ByRef U_index()
As Double, ByRef i As Integer)
Dim a, b, c, X As Double
Dim s, t, r As Integer
Dim N As Integer

Try'数値がある行を探索
For i = 1 To 300
If IsNumeric(Sheet.Cells(i, 1).value) = True Then
N=i
Exit ForEnd If45

Next'空セルまで繰り返す
For i = N To 300
If Sheet.Cells(i + 1, 1).value = Nothing Then Exit For'時間間隔
s = Sheet.Cells(i, 1).value 't1
t = Sheet.Cells(i + 1, 1).value 't2
r = t - s 't2 - t1
'幹樹材積の差
a = Sheet.Cells(i, 2).value 'w1
b = Sheet.Cells(i + 1, 2).value 'w2
c = b - a 'w2 -w1
X = Math.Log(b) - Math.Log(a) 'lnw2-lnw1'U-W の値
W_index(i - N) = c / X
U_index(i - N) = (X ^ 2) / (c * r)
Next
i=i-N
Catch
MsgBox("セルの値、または値の場所が正しくありません。")
End Try
End Sub'CSV ファイルの値から W と U の値を計算
Public Sub figure(ByVal sr As System.IO.StreamReader, ByRef W_index() As Double,
ByRef U_index() As Double, ByRef i As Integer)
Dim dat As String
Dim sbuf(1) As String
Dim delim() As Char = {","c}Dim data(30, 1) As Double
Dim a, b, c, X As Double
Dim s, t, r As IntegerTry
i=0'空セルまで繰り返す
46

Do While sr.Peek() >= 0
'1 行分
dat = sr.ReadLine()sbuf = dat.Split(delim)data(i, 0) = CDbl(sbuf(0))
data(i, 1) = CDbl(sbuf(1))
i += 1Loop
For j = 0 To i - 2
If data(j + 1, 0) = False Then Exit For'時間間隔
s = data(j, 0) 't1
t = data(j + 1, 0) 't2
r = t - s 't2 - t1
'幹樹材積の差
a = data(j, 1) 'w1
b = data(j + 1, 1) 'w2
c = b - a 'w2 -w1
X = Math.Log(b) - Math.Log(a) 'lnw2-lnw1'U-W の値
W_index(j) = c / X
U_index(j) = (X ^ 2) / (c * r)
Next
i=i-1
Catch
MsgBox("セルの値、または値の場所が正しくありません。")
End Try
End Sub'計算結果をプロット
Sub PlotPoint(ByVal g1 As Graph, ByVal W_index() As Double, ByVal U_index() As
Double, ByRef j As Integer, ByVal k1 As Integer, ByVal k2 As Integer)
Try
For i = 0 To j - 1
If k1 - 1 <= i And k2 - 1 >= i Then
g1.Plot(Math.Log10(W_index(i)), Math.Log10(U_index(i)), Color.Blue, 4,
47

3)
Else
g1.Plot(Math.Log10(W_index(i)), Math.Log10(U_index(i)), Color.Blue, 4,
2)
End IfNext
Catch ex As Exception
MsgBox(ex.Message)
End Try
End Sub'k1 と k2 を返さない PlotPoint
Sub PlotPoint(ByVal g1 As Graph, ByVal W_index() As Double, ByVal U_index() As
Double, ByRef j As Integer)
Try
For i = 0 To j - 1
g1.Plot(Math.Log10(W_index(i)), Math.Log10(U_index(i)), Color.Blue, 4, 2)
Next
Catch ex As Exception
MsgBox(ex.Message)
End Try
End Sub'拡大図用
Sub ZoomPoint(ByVal g1 As Graph, ByVal W_index() As Double, ByVal U_index() As
Double, ByRef j As Integer, ByVal z1 As Integer, ByVal z2 As Integer)
Try
For i = z1 To z2 - 1
If i + 1 > i And z2 <= j Then
g1.Plot(Math.Log10(W_index(i)), Math.Log10(U_index(i)), Color.Blue, 4,
2)
End IfNext
Catch ex As Exception
MsgBox(ex.Message)
End Try
End Sub'残差平方和48

Function LSM(ByVal A As Double, ByVal c As Double, ByVal W As Double, ByVal m As
Double, ByRef W_index() As Double, ByRef U_index() As Double, ByVal K2 As Integer, ByVal
K1 As Double) As Double
Dim X, fx, u As Double
Dim sum As Double = 0'誤差を計算する範囲を設定
For i = K1 - 1 To K2 - 1 Step 1
X = W_index(i)fx = (A * (X ^ -c)) * (1 - (X / W) ^ m)
u = U_index(i)sum += ((Math.Log10(fx) - Math.Log10(u)) ^ 2)
Next
'差の二乗の和を返す
Return sum
End Function'未使用
Function LSM_CS(ByVal A As Double, ByVal c As Double, ByVal W As Double, ByVal m
As Double, ByRef W_index() As Double, ByRef U_index() As Double, ByVal K2 As Integer,
ByVal K1 As Double) As Double
Dim X, fx, u, xx As Double
Dim sum As Double = 0
Dim sum2 As Double = 0
Dim xbar, sigma2, sigma As Double'誤差を計算する範囲を設定
For i = K1 - 1 To K2 - 1 Step 1
X = W_index(i)fx = (A * (X ^ -c)) * (1 - (X / W) ^ m)
u = U_index(i)'xx = Math.Abs((Math.Log10(fx) - Math.Log10(u)))
xx = (Math.Log10(fx) - Math.Log10(u)) ^ 2
sum += xx
sum2 += xx ^ 2

49

Next
xbar = sum / (K2 - K1 + 1)
sigma2 = (sum2 - (K2 - K1 + 1) * xbar ^ 2) / (K2 - K1 + 1)
sigma = Math.Sqrt(sigma2)'差の二乗の和を返す
Return sigma
End Function'曲線を引く処理
Sub curve(ByVal g1 As Graph, ByVal A As Double, ByVal c As Double, ByVal W As
Double, ByVal m As Double, ByVal x_min As Double, ByVal x_max As Double)
Dim X, Y As Double'曲線の一番最初の点
X = x_min'テキストボックスの値が正しくないときの処理
If W <= X Or Math.Abs(A) = 0 Or c = 0 Or m = 0 Then Exit Sub
Y = (A * (X ^ -c)) * (1 - (X / W) ^ m)
g1.PutPoint(Math.Log10(X), Math.Log10(Y), Color.Black)'曲線を続けて引く処理
For n = 0 To Math.Log10(x_max) - Math.Log10(x_min) - 1
For X = x_min * 10 ^ n To x_min * 10 ^ (n + 1) Step x_min * 10 ^ (n - 2)
Y = (A * (X ^ -c)) * (1 - (X / W) ^ m)
'w が上限値を超えるまで線を引く
If W > X Then
g1.LineTO(Math.Log10(X), Math.Log10(Y), Color.Red)
ElseIf W <= X Then
g1.LineTO(Math.Log10(X), -20, Color.Red)
Exit SubEnd IfNext
Next
End Sub'UW の値によってグラフを変える処理
Sub max_min(ByRef x_min As Double, ByRef x_max As Double, ByRef y_min As Double,
ByRef y_max As Double, ByVal x_min2 As Double, ByVal x_max2 As Double, ByVal y_min250

As Double, ByVal y_max2 As Double)
x_min = xmin(x_min, x_min2)
x_max = xmax(x_max, x_max2)
y_min = ymin(y_min, y_min2)
y_max = ymax(y_max, y_max2)
End Sub
Function keta(ByVal x As Double, ByVal N As Integer, ByVal c As Integer)
Dim t, m As Double'ｘをｎ桁かつ約（ｃ＝１：切りあげ，ｃ＝０：四捨五入：ｃ＝－１：切り捨て）
t = 0.5 * c + 0.5
m = Int(Math.Log(x) / Math.Log(10))
keta = Int(10 ^ (N - 1 - m) * x + t) * 10 ^ (m + 1 - N)
End Function'以下 4 つ初期値を変える処理
Function xmin(ByVal x_min As Double, ByVal min As Double) As Double
For n = -10 To 10
For i = 10 ^ n To 10 ^ (n + 1) Step 10 ^ n
If min < i Thenx_min = 10 ^ nReturn x_min
End IfNext
Next
End Function
Function xmax(ByVal x_max As Double, ByVal max As Double) As Double
For n = -10 To 10
For i = 10 ^ n To 10 ^ (n + 1) Step 10 ^ n
If max < i Thenx_max = 10 ^ (n + 2)
Return x_max
End IfNext
Next
End Function
Function ymin(ByVal y_min As Double, ByVal min As Double) As Double
For n = -10 To 10
For i = 10 ^ n To 10 ^ (n + 1) Step 10 ^ n
If min < i Then51

y_min = 10 ^ nReturn y_min
End IfNext
Next
End Function
Function ymax(ByVal y_max As Double, ByVal max As Double) As Double
For n = -10 To 10
For i = 10 ^ n To 10 ^ (n + 1) Step 10 ^ n
If max < i Theny_max = 10 ^ (n + 1)
Return y_max
End IfNext
Next
End Function
End Module52