旋轉矩陣程式記錄

三度空間中一個點， 分別繞 X Y Z 軸旋轉 Ω Φ Κ 弧度， 結果可以表示成旋轉矩陣乘上該點的座標向量； 實作中需要用到矩陣的相乘， 將多次旋轉的相乘，用一個矩陣代表三次轉旋。

低階的 C 語言

一開始用 C 寫，結果覺得宣告矩陣超麻煩。 我是用二維陣列模擬矩陣， 但 C 又只能在初始化陣列時對多個陣列賦值。 覺得不管怎麼作都很麻煩， 在副函數宣告的矩陣又不能回傳， 試了很多寫法都不滿意。

C 的問題是對底層過於嚴格， 我認為需要親自管理記憶體的語言都太低階， 不想用。

JS 物向導向

後來改用我唯一熟練的語言 JavaScript 寫， 就她不尋常的 物件導向 風格來寫。 宣告一個 matrix class ， 有 public 方法 multiply ， 可以接受另一個 matrix 作輸入， 返回一個新的 matrix ； 所以可以像 jQuery 鏈式呼叫 ： m4 = m1.multiply(m2).multiply(m3) 。

用 matlab 寫，還可以用匿名函數， 一行搞定一種矩陣宣告。 mw = @(w) [1 0 0; 0 cos(w) sin(w); 0 -sin(w) cos(w)] 可惜老師規定不能用 matlab QQ

其實我是用 JS 的方言 CoffeeScript 寫的， 然後再編譯成 JS ， 最原始是 wfk.coffee ， 編譯好的 JS 是 wfk.js 。

可以直接引用原始碼， 看要用 nodejs require 或直接在網頁中引用。 或用我包好的版本： 計算旋轉矩陣 。 （媽的超麻煩，寫這個浪費了我 寫核心程式 2 到 3 倍的時間。）

class Matrix3
# Matrix3 is 2 demition array,
# access m(i,j) by m[i][j].

    constructor: (a = 0) ->
    # a = number, [1,2,3,4,5,6,7,8,9], [[1,2,3],[4,5,6],[7,8,9]].
    # if no arguments, a = 0.
        a = [1..9].map( -> a ) if typeof a == 'number'

        if a.length >= 9
            for i in [0..2]
                this[i] = a.slice i*3,(i+1)*3
        else if a.length >= 3
            for i in [0..2]
                this[i] = a[i].slice 0
        else
            throw new Error 'arguments not array!'

    ':': (n) ->
    # m[':'](n) to access nth column of m, 
    # just like syntax of matlab. 
        this.reduce(
            (s,v,i) ->
                s.push v if i[1] == n
                return s
            []
        )

    size: [3,3] # Matrix3 is [3,3] matrix.

    forEach: (callback) ->
    # callback will get 3 argument v,i,m, 
    # v = current value, i = [row, col], m = matrix it self.
    # row first order.
        for i in [0..2]
            for j in [0..2]
                callback this[i][j], [i,j], this
        return undefined

    map: (callback) ->
    # map will return a new Matrix3 instance.
        m = new Matrix3()
        this.forEach (v,i,a) -> m[i[0]][i[1]] = callback v,i,a
        return m

    reduce: (callback, sum) ->
        this.forEach (v,i,a) -> sum = callback sum, v, i, a
        return sum

    multiply: (m) ->
    # matrix to matrix multiply, return a new matrix.
        @map (v,i,a) ->
            a[i[0]].reduce(
                (s,v,j) -> s + v*m[j][i[1]]
                0
            )

    add: (m) -> @map (v,i) -> v + m[i[0]][i[1]]
    # add by matrix, return a new matrix.

    # convert Matrix3 into 1 demantion array, 
    # and slice.
    slice: (start, end) ->
        (this.reduce ((s,v) -> s.push v), []).slice start, end


sin = Math.sin
cos = Math.cos

# define omega, phi, kappa rotate matrix.
rotateMatrix =
    omega: (w) -> new Matrix3 [
        [1, 0, 0]
        [0, cos(w), sin(w)]
        [0, -sin(w), cos(w)]
    ]
    phi: (f) -> new Matrix3 [
        [cos(f), 0, -sin(f)]
        [0, 1, 0]
        [sin(f), 0, cos(f)]
    ]
    kappa: (k) -> new Matrix3 [
        [cos(k), sin(k), 0]
        [-sin(k), cos(k), 0]
        [0, 0, 1]
    ]

# define a function directly return 
# a matrix rotate by omega, phi, kappa. 
rotateWFKMatrix = (w, f, k) ->
    rotateMatrix.kappa k
        .multiply rotateMatrix.phi f
        .multiply rotateMatrix.omega w


# check is excute in nodejs or browser.
if  typeof exports == 'object'
    md = exports
else if typeof window == 'object'
    window.rotateWFK = {}
    md = window.rotateWFK
else
    md = {}

# export modules.
md.Matrix3 = Matrix3
md.rotateMatrix = rotateMatrix
md.rotateWFKMatrix = rotateWFKMatrix

轉換成果

成果和 australis 的計算結果相比， 誤差大約在 1E-9 ， 而 australis 的輸出小數到 1E-8 而已， 我的計算結果則是到 1E-16 ， 可以視為資料型態轉換造成的誤差。 （比較上是把矩陣平展，再逐項比較。 反正都差很少，也看不出系統誤差。）