N88-BASICで行列(matrix) (2回目)

2021/6/24(木)

N88-BASICで行列(matrix) (2回目)

 

(by ULproject for N88-BASIC, NL-BASIC)

今回は、前回の併進移動に回転移動を追加します

(キー操作はプログラムの最終行コメント参照)

 

前回同様、行列は大文字、ベクトルとスカラーは

小文字で、行列は4×4(4行4列)とし、ベクトルは

3次元(4次元目は常に1)とします

またベクトル上部の矢印は省いています

 

まず、位置ベクトルp=(x,y,z)を単位ベクトル

n=(n₁,n₂,n₃)まわりにθ度回転した

点p'=(x',y',z')を求める式は

 

図1 位置ベクトルの回転

 

|p'- q|=|p - q|, |a| = |p - q|cosθ

a  = |a|(p - q)/|p - q|

   = (p - q)cosθ = pcosθ - qcosθ

 

|n×p|はnとpでできる平行四辺形の

面積なので底辺|n|=1と高さ|p - q|の積

|n×p| = |n||p - q| = |p - q|より、

b = |p - q|sinθ(n×p)/|n×p|

  = (n×p)sinθ

q = (n･p)n

p'= q + a + b

  = q+pcosθ-qcosθ+(n×p)sinθ

  = pcosθ+(1-cosθ)q+(n×p)sinθ

 

p'= pcosθ+(1-cosθ)(n･p)n+(n×p)sinθ

を、p' =Rpの形に変形していきます

 

外積a×bと内積a･bの成分表記について、

x,y,z軸方向の単位ベクトルをi,j,k、

a=(x,y,z),b=(x',y',z')とすると、

 

i×i = j×j = k×k = 0ベクトル

i×j =  k , j×k =  i , k×i =  j

j×i = -k , k×j = -i , i×k = -j

より

a×b

= (ix + jy + kz)×(ix'+ jy'+ kz')

= kxy'-jxz'-kyx'+iyz'+jzx'-izy'

= (yz'-zy')i+(zx'-xz')j+(xy'-yx')k

= (yz'- zy', zx'- xz', xy'- yx')

 

i･i = j･j = k･k = 1、これ以外0より、

a･b

= (ix + jy + kz)･(ix' + jy' + kz')

= xx' + yy' + zz'

 

a×bは、大きさ|a×b|がaとbで出来る

平行四辺形の面積で、aとbに垂直でaから

bに右手の親指以外を向けた時、親指の方向

 

a･b = |a||b|cosθで、|b|をaに垂直に

投影した長さと|a|の積となります

(θはaとbのなす角)

 

p = (x,y,z)より、

pcosθ= (xcosθ, ycosθ, zcosθ)

 |xcosθ| |cosθ   0     0  ||x|

=|ycosθ|=|  0   cosθ   0  ||y|

 |zcosθ| |  0     0   cosθ||z|

 

n = (n₁,n₂,n₃)より、

(n･p)n

 |(n₁x + n₂y + n₃z)n₁|

=|(n₁x + n₂y + n₃z)n₂|

 |(n₁x + n₂y + n₃z)n₃|

 |n₁n₁ n₁n₂ n₁n₃||x|

=|n₂n₁ n₂n₂ n₂n₃||y|

 |n₃n₁ n₃n₂ n₃n₃||z|

 

n×p

 |n₂z - n₃y| | 0  -n₃  n₂||x|

=|n₃x - n₁z|=| n₃  0  -n₁||y|

 |n₁y - n₂x| |-n₂  n₁  0 ||z|

 

cosθ+n₁n₁(1-cosθ) = n₁n₁+(1-n₁n₁)cosθ

 

p' = pcosθ+(1-cosθ)(n･p)n+(n×p)sinθ

は、p'=Rpとして、

  |r₁₁ r₁₂ r₁₃|

R=|r₂₁ r₂₂ r₂₃|

  |r₃₁ r₃₂ r₃₃|

r₁₁ = n₁n₁+(1-n₁n₁)cosθ

r₂₂ = n₂n₂+(1-n₂n₂)cosθ

r₃₃ = n₃n₃+(1-n₃n₃)cosθ

r₁₂ = n₁n₂(1-cosθ)-n₃sinθ

r₂₁ = n₂n₁(1-cosθ)+n₃sinθ

r₁₃ = n₁n₃(1-cosθ)+n₂sinθ

r₃₁ = n₃n₁(1-cosθ)-n₂sinθ

r₂₃ = n₂n₃(1-cosθ)-n₁sinθ

r₃₂ = n₃n₂(1-cosθ)+n₁sinθ

となり、回転行列Aは

  |r₁₁ r₁₂ r₁₃ 0|

A=|r₂₁ r₂₂ r₂₃ 0|

  |r₃₁ r₃₂ r₃₃ 0|

  |0   0   0   1|

です

 

ベクトル(x,y,z)軸回りにd°回転するプログラム

(角度は度→ラジアン、(x,y,z)→単位ベクトルn)

 

*ROTATE

PI = 3.141592653589793

RR = SQR(X*X + Y*Y + Z*Z)

N1 = X / RR

N2 = Y / RR

N3 = Z / RR

TT = D * PI / 180

SS = COS(TT)

RR = 1.0 - SS

A( 0) = N1 * N1 * RR + SS

A( 5) = N2 * N2 * RR + SS

A(10) = N3 * N3 * RR + SS

SS = SIN(TT)

TT = N1 * N2 * RR: NN = N3 * SS: A(4) = TT - NN: A(1) = TT + NN

TT = N1 * N3 * RR: NN = N2 * SS: A(2) = TT - NN: A(8) = TT + NN

TT = N2 * N3 * RR: NN = N1 * SS: A(9) = TT - NN: A(6) = TT + NN

A(12) = 0: A(13) = 0: A(14) = 0: A(15) = 1

A( 3) = 0: A( 7) = 0: A(11) = 0

GOSUB *MULT.MA

RETURN

となります

NL-BASICとblg～.zip(mat002.bas)は

以下のリンクからダウンロードできます

NL-BASIC(N88-BASIC互換?)ホームページ

Readme.txtを読んで遊んで下さい

次回のおすすめ

N88-BASICでクォータニオン(quaternion) (1回目)