直列多重剛体振り子

(問題) $n$個の剛体が直列に連結された多重剛体振り子(Series-Multiple Rigid Pendulum)がある。 これらは重力と逆向きの方向を$y$軸とする$xy$平面上で回転運動するものとする。 剛体1は質量$m_1$で、原点$O$に固定された回転軸$O_1$をもち、 $O_1$から剛体1の重心$G_1$までの距離を$c_1$とする。 また、剛体1の$G_1$周りの慣性モーメントを$I_{G_1}$とする。 その他の剛体 $i(i=2,\ldots\,n)$は質量$m_i$で、 剛体$i-1$上の$O_{i-1}$から距離$b_{i-1}$、 角度 $\angle G_{i-1}O_{i-1}O_i=\beta_{i-1}$の回転軸$O_i$をもち、 $O_i$から剛体$i$の重心$G_i$までの距離を$c_i$とする。 また、剛体$i$の$G_i$周りの慣性モーメントを$I_{G_i}$とする。 この多重剛体振り子のラグランジュ運動方程式を記述せよ。 また、計算機シミュレーションなどをしやすいように、 求められた2階常微分方程式を連立の1階常微分方程式に変換せよ。 また、この振り子が各々の角速度に比例する抵抗 や外力 を受けている場合においても同様に求めよ。

(解答) 剛体 $i(i=1,\ldots,n)$のy軸に対する角度を$\theta_i$とする。 剛体 $i(i=1,\ldots,n)$の重心$G_i$の座標 $(x_{G_i},y_{G_i})$およびその速度$V_{G_i}$の2乗は、

$$x_{G_i} = \sum_{j=1}^{i-1} b_j\sin(\theta_j+\beta_j) + c_i\sin\theta_i$$ (187)

$$y_{G_i} = -\sum_{j=1}^{i-1} b_j\cos(\theta_j+\beta_j) - c_i\cos\theta_i$$ (188)

$$V_{G_i}^2 = \dot x_{G_i}^2 + \dot y_{G_i}^2$$

$$= \left\{ \sum_{j=1}^{i-1} b_j\dot\theta _j\cos(\theta_j+\beta_j) +... ..._j\dot\theta _j\sin(\theta_j+\beta_j) + c_i\dot\theta _i\sin\theta_i \right\}^2$$

$$= \sum_{j=1}^{i-1} b_j^2\dot\theta _j^2 + c_i^2\dot\theta _i^2 + 2\... ...sum_{j=1}^{i-1} b_jc_i\dot\theta _j\dot\theta _i\cos(\theta_j+\beta_j-\theta_i)$$

となる。 ここで、剛体 $i(i=1,\ldots,n)$の運動エネルギー$T_i$および位置エネルギー$U_i$は、

$$T_i = \frac12 I_{G_i} \dot\theta_i^2 + \frac12 m_i V_{G_i}^2$$

$$= \frac12 I_{G_i} \dot\theta _i^2 + \frac12 m_i \left( \sum_{j=1}^{i-1} b_j^2\dot\theta _j^2 + c_i^2\dot\theta _i^2 \right)$$

$$+ m_i \sum_{j=1}^{i-1} \sum_{k=j+1}^{i-1} b_jb_k\dot\theta _j\dot... ...sum_{j=1}^{i-1} b_jc_i\dot\theta _j\dot\theta _i\cos(\theta_j+\beta_j-\theta_i)$$ (189)

$$U_i = m_i g y_{G_i}$$

$$= -m_i g \left\{ \sum_{j=1}^{i-1} b_j\cos(\theta_j+\beta_j) + c_i\cos\theta_i \right\}$$ (190)

となる。 この多重剛体振り子のラグランジュ関数$L$は、

$$L = T - U$$

$$= \sum_{i=1}^nT_i - \sum_{i=1}^nU_i$$

$$= \sum_{i=1}^n \left\{ \frac12 I_{G_i} \dot\theta _i^2 + \frac12 m_... ...( \sum_{j=1}^{i-1} b_j^2\dot\theta _j^2 + c_i^2\dot\theta _i^2 \right) \right\}$$

$$+ \sum_{i=1}^n \left\{ m_i \sum_{j=1}^{i-1} \sum_{k=j+1}^{i-1} b_... ...^{i-1} b_jc_i\dot\theta _j\dot\theta _i\cos(\theta_j+\beta_j-\theta_i) \right\}$$

$$+ \sum_{i=1}^n m_i g \left\{ \sum_{j=1}^{i-1} b_j\cos(\theta_j+\beta_j) + c_i\cos\theta_i \right\}$$

$$= \frac12 \sum_{i=1}^n \left( I_{G_i} + m_i c_i^2 \right) \dot\theta _i^2 + \frac12 \sum_{i=1}^n m_i \sum_{j=1}^{i-1} b_j^2\dot\theta _j^2$$

$$+ \sum_{i=1}^n m_i \sum_{j=1}^{i-1} \sum_{k=j+1}^{i-1} b_jb_k\dot... ...sum_{j=1}^{i-1} b_jc_i\dot\theta _j\dot\theta _i\cos(\theta_j+\beta_j-\theta_i)$$

$$+ \sum_{i=1}^n m_i g \sum_{j=1}^{i-1} b_j\cos(\theta_j+\beta_j) + \sum_{i=1}^n m_i g c_i\cos\theta_i$$

$$= \frac12 \sum_{i=1}^n \left( I_{G_i} + m_i c_i^2 + \sum_{k=i+1}^n m_k b_i^2 \right) \dot\theta _i^2$$

$$+ \frac12 \sum_{i=1}^n \sum_{j=1}^{i-1} \left\{ \sum_{k=i+1}^n m_... ... + m_i b_jc_i\dot\theta _i\dot\theta _j\cos(\theta_j+\beta_j-\theta_i) \right\}$$

$$+ \frac12 \sum_{i=1}^n \sum_{j=i+1}^{n} \left\{ \sum_{k=j+1}^n m_... ... + m_j b_ic_j\dot\theta _i\dot\theta _j\cos(\theta_i+\beta_i-\theta_j) \right\}$$

$$+ \sum_{i=1}^n \left\{ \sum_{k=i+1}^n m_k g b_i\cos(\theta_i+\beta_i) + m_i g c_i\cos\theta_i \right\}$$ (191)

となる。 このラグランジュ関数$L$から $\theta_d(d=1,\ldots,n)$ に関するラグランジュ運動方程式をたてるには、 以下の関係式を用いればよい。

$$\begin{eqnarray*} \frac{d}{dt} \frac{\partial L}{\partial\dot\theta _d} - \frac{\partial L}{\partial\theta_d} &=& 0 \end{eqnarray*}$$

これを具体的に求めると、

$$\frac{\partial L}{\partial\dot\theta_d} = (I_{G_d} + m_d c_d^2 + m_{d+1..n} b_d^2)\dot\theta _d$$

$$+ \frac12 \sum_{\stackrel{i=1}{d=i}}^n \sum_{\stackrel{j=1}{d=j}}... ...artial}{\partial\dot\theta _d}\dot\theta _i\dot\theta _j{\cal C}_{j,i} \right\}$$

$$+ \frac12 \sum_{\stackrel{i=1}{d=i}}^n \sum_{\stackrel{j=i+1}{d=j... ...artial}{\partial\dot\theta _d}\dot\theta _i\dot\theta _j{\cal C}_{i,j} \right\}$$

$$= (I_{G_d} + m_d c_d^2 + m_{d+1..n} b_d^2)\dot\theta _d$$

$$+ \sum_{i=1}^{d-1} m_{d+1..n} b_db_i\dot\theta _iC_{d,i} + \sum_{... ...theta _i{\cal C}_{i,d} + \sum_{i=d+1}^{n} m_i b_dc_i\dot\theta _i{\cal C}_{d,i}$$ (192)

$$\Rightarrow \frac{d}{dt} \frac{\partial L}{\partial\dot\theta_d} = (I_{G_d} + m_d c_d^2 + m_{d+1..n} b_d^2)\ddot\theta _d$$

$$+ \sum_{i=1}^{d-1} m_{d+1..n} b_db_i \left\{\ddot\theta _iC_{d,i}... ...dot\theta _iC_{d,i} - \dot\theta _i(\dot\theta _d-\dot\theta _i)S_{d,i}\right\}$$

$$+ \sum_{i=1}^{d-1} m_d b_ic_d \left\{\ddot\theta _i{\cal C}_{i,d}... ...cal C}_{d,i} - \dot\theta _i(\dot\theta _d-\dot\theta _i){\cal S}_{d,i}\right\}$$

$$= (I_{G_d} + m_d c_d^2 + m_{d+1..n} b_d^2)\ddot\theta _d$$

$$+ \sum_{i=1}^{d-1} m_{d+1..n} b_db_i \ddot\theta _iC_{d,i} - \sum_{i=1}^{d-1} m_{d+1..n} b_db_i \dot\theta _i(\dot\theta _d-\dot\theta _i)S_{d,i}$$

$$+ \sum_{i=d+1}^{n} m_{i+1..n} b_db_i \ddot\theta _iC_{d,i} - \sum_{i=d+1}^{n} m_{i+1..n} b_db_i \dot\theta _i(\dot\theta _d-\dot\theta _i)S_{d,i}$$

$$+ \sum_{i=1}^{d-1} m_d b_ic_d \ddot\theta _i{\cal C}_{i,d} {+}\sum_{i=1}^{d-1} m_d b_ic_d \dot\theta _i(\dot\theta _d-\dot\theta _i){\cal S}_{i,d}$$

$$+ \sum_{i=d+1}^{n} m_i b_dc_i \ddot\theta _i{\cal C}_{d,i} - \sum_{i=d+1}^{n} m_i b_dc_i \dot\theta _i(\dot\theta _d-\dot\theta _i){\cal S}_{d,i}$$

$$= (I_{G_d} + m_d c_d^2 + m_{d+1..n} b_d^2)\ddot\theta _d$$

$$+ \sum_{i=1}^{d-1} m_{d+1..n} b_db_i \ddot\theta _iC_{d,i} - \sum... ...dot\theta _iS_{d,i} + \sum_{i=1}^{d-1} m_{d+1..n} b_db_i \dot\theta _i^2S_{d,i}$$

$$+ \sum_{i=d+1}^{n} m_{i+1..n} b_db_i \ddot\theta _iC_{d,i} - \sum... ...dot\theta _iS_{d,i} + \sum_{i=d+1}^{n} m_{i+1..n} b_db_i \dot\theta _i^2S_{d,i}$$

$$+ \sum_{i=1}^{d-1} m_d b_ic_d \ddot\theta _i{\cal C}_{i,d} {+}\su... ...a _i{\cal S}_{i,d} {-}\sum_{i=1}^{d-1} m_d b_ic_d \dot\theta _i^2{\cal S}_{i,d}$$

$$+ \sum_{i=d+1}^{n} m_i b_dc_i \ddot\theta _i{\cal C}_{d,i} - \sum... ...ta _i{\cal S}_{d,i} + \sum_{i=d+1}^{n} m_i b_dc_i \dot\theta _i^2{\cal S}_{d,i}$$ (193)

$$- \frac{\partial L}{\partial\theta_d} = \frac12 \sum_{\stackrel{i=1}{d=i}}^n \sum_{\stackrel{j=1}{d=j}}^{... ...t\theta _i\dot\theta _j\frac{\partial}{\partial\theta_d}{\cal C}_{j,i} \right\}$$

$$+ \frac12 \sum_{\stackrel{i=1}{d=i}}^n \sum_{\stackrel{j=i+1}{d=j... ...{i,j} \right\} + m_{d+1..n} g b_d\sin(\theta_d+\beta_d) + m_d g c_d\sin\theta_d$$

$$= \sum_{i=1}^{d-1} m_{d+1..n} b_db_i \dot\theta _d\dot\theta _iS_{d,i} + \sum_{i=d+1}^{n} m_{i+1..n} b_db_i \dot\theta _d\dot\theta _iS_{d,i}$$

$${-}\sum_{i=1}^{d-1} m_d b_ic_d \dot\theta _d\dot\theta _i{\cal S}... ...{\cal S}_{d,i} + m_{d+1..n} g b_d\sin(\theta_d+\beta_d) + m_d g c_d\sin\theta_d$$ (194)

但しここで、

$$m_{k..l} = \sum_{j=k}^l m_j$$

$$C_{i,j} = \cos(\theta_i+\beta_i-\theta_j-\beta_j) = \cos(\beta_i-\beta_j)\cos(\theta_i-\theta_j)-\sin(\beta_i-\beta_j)\sin(\theta_i-\theta_j)$$

$$S_{i,j} = \sin(\theta_i+\beta_i-\theta_j-\beta_j) = \sin(\beta_i-\beta_j)\cos(\theta_i-\theta_j)+\cos(\beta_i-\beta_j)\sin(\theta_i-\theta_j)$$

$${\cal C}_{i,j} = \cos(\theta_i+\beta_i-\theta_j) = \cos\beta_i\cos(\theta_i-\theta_j)-\sin\beta_i\sin(\theta_i-\theta_j)$$

$${\cal S}_{i,j} = \sin(\theta_i+\beta_i-\theta_j) = \sin\beta_i\cos(\theta_i-\theta_j)+\cos\beta_i\sin(\theta_i-\theta_j)$$

ゆえに、 $\theta_d(d=1,\ldots,n)$に関するラグランジュ運動方程式は、

$$(I_{G_d} + m_d c_d^2 + m_{d+1..n} b_d^2)\ddot\theta _d$$

$$+ \sum_{i=1}^{d-1} m_{d+1..n} b_db_i \ddot\theta _iC_{d,i} + \sum... ...eta _i{\cal C}_{i,d} + \sum_{i=d+1}^{n} m_i b_dc_i \ddot\theta _i{\cal C}_{d,i}$$

$$+ \sum_{i=1}^{d-1} m_{d+1..n} b_db_i \dot\theta _i^2S_{d,i} + \su... ... _i^2{\cal S}_{i,d} + \sum_{i=d+1}^{n} m_i b_dc_i \dot\theta _i^2{\cal S}_{d,i}$$

$$+ m_{d+1..n} g b_d\sin(\theta_d+\beta_d) + m_d g c_d\sin\theta_d = 0$$ (195)

となる¹¹。 $\Box$

これを連立1階微分方程式にするために、この2階微分方程式を以下のように表す。

$${R}\bm{\ddot\theta } = {\bm{v}} \qquad\mbox{(${R}\in\mathbb{R}^{n\times n},\quad \bm{\ddot\theta },{\bm{v}}\in\mathbb{R}^n$)}$$ (196)

さらに、

$$\vartheta_d = \dot\theta _d$$

とおくことにより以下の連立1階微分方程式を得る。

$$\dot\theta_d = \vartheta_d$$ (197)

$$\dot\vartheta_d = [{R}^{-1}{\bm{v}}]_d$$ (198)

但しここで、

$${R} = \left(\begin{array}{ccccc} I_{1..n}& & & & \\ &\ddots & &m_i b_d... ...d} + m_{d+1..n} b_db_i C_{d,i}& &\ddots & \\ & & & &I_{n..n}\end{array}\right)$$ (199)

$${\bm{v}} = \left[\begin{array}{c} \vdots \\ \multicolumn{1}{l}{\displaystyl... ...i} - m_{d+1..n} g b_d\sin(\theta_d+\beta_d) }}\\ \vdots \\ \end{array}\right]$$ (200)

$$I_i = I_{G_i} + m_i c_i^2$$

$$I_{i..n} = I_i + m_{i+1..n} b_i^2$$

である。 $\Box$ また、各々の角速度に比例する抵抗などの減衰 や強制振動のような外力 がある場合は、 $\dot\theta _d(d=1,\ldots,n)$それぞれに対する減衰定数を $\lambda_d$、外力を$\sigma_d$ とすると、 $\theta_d(d=1,\ldots,n)$に関するラグランジュ運動方程式は式 (199)から、

$$(I_{G_d} + m_d c_d^2 + m_{d+1..n} b_d^2)\ddot\theta _d$$

$$+ \sum_{i=1}^{d-1} m_{d+1..n} b_db_i \ddot\theta _iC_{d,i} + \sum... ...eta _i{\cal C}_{i,d} + \sum_{i=d+1}^{n} m_i b_dc_i \ddot\theta _i{\cal C}_{d,i}$$

$$+ \sum_{i=1}^{d-1} m_{d+1..n} b_db_i \dot\theta _i^2S_{d,i} + \su... ... _i^2{\cal S}_{i,d} + \sum_{i=d+1}^{n} m_i b_dc_i \dot\theta _i^2{\cal S}_{d,i}$$

$$+ m_{d+1..n} g b_d\sin(\theta_d+\beta_d) + m_d g c_d\sin\theta_d = \sigma_d - \lambda_d\dot\theta _d$$ (201)

となる¹²。 $\Box$

よって、式(201),(202)に対応する減衰 や外力 がある場合の連立一階微分方程式は以下のようになる。

$$\dot\theta_d = \vartheta_d$$ (202)

$$\dot\vartheta_d = [{R}^{-1}{\bm{w}}]_d$$ (203)

但しここで、

$${R} = \left(\begin{array}{ccccc} I_{1..n}& & & & \\ &\ddots & &m_i b_d... ...d} + m_{d+1..n} b_db_i C_{d,i}& &\ddots & \\ & & & &I_{n..n}\end{array}\right)$$ (204)

$${\bm{w}} = \left[\begin{array}{c} \vdots \\ \multicolumn{1}{l}{\displaystyl... ...aystyle{ + \sigma_d - \lambda_d \vartheta_d }}\\ \vdots \\ \end{array}\right]$$ (205)

である。 $\Box$ ［補足］ ハミルトンの正準方程式の導出

ここで、運動エネルギー$T$を以下のように表す。

$$T = \frac12 \dot{\bm\theta}^\top {R} \dot{\bm\theta}$$ (206)

但しここで、 $\bm\theta\in\mathbb{R}^n$である。すると、角運動量$\bm{p}$ は、 と表される。ハミルトニアン$H$は

$$H = \frac12 \bm{p}^\top {R}^{-1} \bm{p} + U$$ (207)

であり、ハミルトニアン$H$からハミルトンの正準方程式を得るには以下

$$\dot{\bm\theta} = \frac{\partial H}{\partial \bm{p}} = \frac{\partial}{\partial \bm{p}} \left( \frac12 \bm{p}^\top {R}^{-1} \bm{p} + U \right)$$

$$= {R}^{-1} \bm{p}$$ (208)

$$\dot{\bm{p}} = - \frac{\partial H}{\partial \bm\theta} + \bm{F} = - \frac{\parti... ...tial \bm\theta} \left( \frac12 \bm{p}^\top {R}^{-1} \bm{p} + U \right) + \bm{F}$$

$$= - \frac12 \frac{\partial}{\partial \bm\theta} \left( \bm{p}^\top {R}^{-1} \bm{p} \right) - \frac{\partial U}{\partial \bm\theta} + \bm{F}$$

$$= \frac12 \frac{\partial}{\partial \bm\theta} \left( \dot{\bm\theta... ...op {R} \dot{\bm\theta} \right) - \frac{\partial U}{\partial \bm\theta} + \bm{F}$$ (209)

の関係式を用いればよい。 ちなみに$\bm{F}$は、前述の抵抗や外力がある場合は、

$$\bm{F} = \left[\begin{array}{c} \vdots\\ \sigma_d - \lambda_d\dot\theta_d \\ \vdots \end{array} \right]$$ (210)

さもなくば${F}_d=0$($d=1,\ldots,n$)である。

具体的に求めると、ハミルトンの正準方程式は以下のようになる。

$$\dot{\bm\theta} = {R}^{-1} \bm{p}$$ (211)

$$\dot{\bm{p}} = \frac12 \frac{\partial}{\partial \bm\theta} \left( \dot{\bm\theta... ...bm\theta}, \ldots \right)^\top - \frac{\partial U}{\partial \bm\theta} + \bm{F}$$

$$= \left[\begin{array}{c} \vdots\\ \multicolumn{1}{l}{\displaystyle... ..._{d+1..n} g b_d\sin(\theta_d+\beta_d) \\ \vdots\\ \end{array}\right] + \bm{F}$$

$\Box$ ［発展］ 動座標系におけるラグランジュ運動方程式及びハミルトンの正準方程式の導出

$O_1$が時変の位置 $\bm{\Gamma}=(\Gamma_x, \Gamma_y)$によって動くことに より、振り子が強制振動される場合を考える。 $i$番目の質点の位置ベクトルを $\bm{r}_i=(x_{G_i},y_{G_i})$として、この 動座標系における運動エネルギー$T$は、

$$T = \frac12\sum_{i=1}^n m_i\vert\dot{\bm{r}}_i\vert^2 + ... ...bm{r}}_i + \frac12\sum_{i=1}^n m_i\vert\dot{\bm\Gamma}\vert^2$$ (212)

となり、位置エネルギー$U$は、

$$U = \sum_{i=1}^n m_i g y_{G_i} + \sum_{i=1}^n m_i g \Gamma_y$$ (213)

となる。ここで、

$$\sum_{i=1}^n m_i\dot{\bm\Gamma}^\top\dot{\bm{r}}_i = \sum_{i=1}^n m_i (\dot\Gamma_x, \dot\Gamma_y)\left[\begin{array}{... ...theta_j\sin(\theta_j+\beta_j) + c_i\dot\theta_i\sin\theta_i} \end{array}\right]$$

$$= \left[\begin{array}{c} \vdots\\ \displaystyle{m_{i+1..n} b_i \le... ...ots\\ \dot\theta_i\\ \vdots \end{array}\right] = \bm{p}_O^\top\dot{\bm\theta}$$

とおくと、ラグランジュ関数$L$は$L=T-U$より、

$$L = \frac12 \dot{\bm\theta}^\top {R} \dot{\bm\theta} + \bm{p}_O^\top\dot{\bm\theta} + \frac12\sum_{i=1}^n m_i\vert\dot{\bm\Gamma}\vert^2 - U$$ (214)

となり、ラグランジュ運動方程式は、

$${R} \ddot{\bm\theta} + \dot{R} \dot{\bm\theta} + \frac{\p... ...heta}\right) + \frac{\partial U}{\partial\bm\theta} = \bm{F}$$ (215)

となる。但しここで、

$$\frac{\partial \bm{p}_O}{\partial t} = \left[\begin{array... ...dot\Gamma_y\sin\theta_i\right)}\\ \vdots \end{array}\right]$$ (216)

である。 この場合の角運動量$\bm{p}$は、 であるので、 より、ハミルトニアン$H$は、

$$H = \frac12 \left( \bm{p} - \bm{p}_O \right)^\top {R}^{-1} ... ...ight) + \frac12\sum_{i=1}^n m_i\vert\dot{\bm\Gamma}\vert^2 + U$$ (217)

であり、ハミルトニアン$H$からハミルトンの正準方程式を得るには以下

$$\dot{\bm\theta} = \frac{\partial H}{\partial \bm{p}} = \frac{\partial}{\partial \bm... ...{p}_O \right) + \frac12\sum_{i=1}^n m_i\vert\dot{\bm\Gamma}\vert^2 + U \right\}$$

$$= {R}^{-1} \left( \bm{p} - \bm{p}_O \right)$$ (218)

$$\dot{\bm{p}} = - \frac{\partial H}{\partial \bm\theta} + \bm{F} = - \frac{\parti... ...ght) + \frac12\sum_{i=1}^n m_i\vert\dot{\bm\Gamma}\vert^2 + U \right\} + \bm{F}$$

$$= - \frac12 \frac{\partial}{\partial \bm\theta} \left\{ \left( \bm{... ...p} - \bm{p}_O \right) \right\} - \frac{\partial U}{\partial \bm\theta} + \bm{F}$$

$$= - \left( \frac12 \left( \bm{p} - \bm{p}_O \right)^\top \frac{\par... ..._O \right) ,\ldots\right)^\top - \frac{\partial U}{\partial \bm\theta} + \bm{F}$$

$$= - \left( - \frac12 \dot{\bm\theta}^\top \frac{\partial {R}}{\part... ...\bm\theta} ,\ldots\right)^\top - \frac{\partial U}{\partial \bm\theta} + \bm{F}$$

$$= \frac12 \frac{\partial}{\partial \bm\theta} \left( \dot{\bm\theta... ...O^\top \dot{\bm\theta} \right) - \frac{\partial U}{\partial \bm\theta} + \bm{F}$$ (219)

の関係式を用いればよい。但しここで、

$$\frac{\partial}{\partial \bm\theta}\left( \bm{p}_O^\top \dot... ...os\theta_i\right) \dot\theta_i}\\ \vdots \end{array}\right]$$ (220)

である。

... となる¹¹

確認のために$n=2$として書き下すと、確かに式 (16),(17)に対応する。

... となる¹²

確認のために$n=2$として書き下すと、確かに式 (26),(27)に対応する。

このサイトについてのご意見ご要望は taiji@aihara.co.jp までお願いします。 Copyright (C) 1996-1998,2003-2004 Taiji Yamada, AIHARA Electrical Engineering Co., Ltd.
All rights reserved.