서윤의 이글루

http://jogamp.org/wiki/index.php/Downloading_and_installing_JOGL

위 페이지로 가서  Downloading a released version  섹션에서 Signed Release  링크를 따라가면

jogl 압축파일이 보인다.  자신의 시스템에 맞는 파일을 받은 다음.
(현재는 jogl-2.0-b23-20110303-macosx-universal.7z)  압축을 풀면,  

폴더들이 보이는데 이중 jar  과 lib 폴더를 사용자의 홈 디렉토리로 복사한다.

http://jogamp.org/wiki/index.php/Using_JOGL_in_AWT_SWT_and_Swing

에서 첫번째와 끝에서 두번째 자바파일을 복사, 저장한 후 , 아래의 명령으로 컴파일한다.

javac -classpath "jar/gluegen-rt.jar:jar/nativewindow.all.jar:jar/jogl.all.jar:jar/newt.all.jar" OneTriangle.java

javac -classpath "jar/gluegen-rt.jar:jar/nativewindow.all.jar:jar/jogl.all.jar:jar/newt.all.jar:." OneTriangleSwingGLCanvas.java

이제 실행시켜본다.

java -classpath "jar/gluegen-rt.jar:jar/nativewindow.all.jar:jar/jogl.all.jar:jar/newt.all.jar:." -Djava.library.path=lib OneTriangleSwingGLCanvas

책이름을 보고, 또 그저그런 잡담류책인가 했었다.

하지만 , 속은 아주 알차다.  고등학교과정까지의 수학들을 일목요연하게 . 빡빡하지 않은 편집으로 정리해놨다.

군데군데 들어있는 수학자들의 이야기가 책을 더 재미있게해준다.

(def mat1 '(3 4 6 7))
(def mat2 '(2 3 4 5))
(def mat3 '(6 2 7   5 8 5   7 4 5))
(def mat4 '(5 3 8 4    8 5 4 9    5 7 2 1    6 4 3 8))

(defn <-m [mat n i j]
    (nth mat (+ j (* (- i 1) n) -1)))

클로저에 있는 nth 함수를 이용하려고 하니 이 함수가 바뀌었다.
두번째 인자인 n 은  m x n 행렬임을 뜻한다.  즉 (2 x 3) 행렬일때 3 이 들어간다.

(defn arr1 [a b c]
    (range a (+ b 1) c))
    
(defn m+ [m1 m2]
    (map  + m1 m2))
   
(defn m- [m1 m2]
    (map  - m1 m2))
   
(defn m*c [m1 a]
    (map #(* a %) m1))
   
행렬의 더하기 , 빼기 , 상수와의 곱은 아주 쉽지만 , 행렬곱을 구하려면 좀 복잡해진다.

(defn mk-pair [arr1 arr2]
    (for [x arr1 y arr2]
        [x y]))

;(mk-pair '(1 2 3 4) '(1 2 3 4))

(defn ksum [k c sum i j m1 n1 m2 n2]
    (if (> c k)
        sum
        (recur k (+ c 1) (+ sum (* (<-m m1 n1 i c) (<-m m2 n2 c j))) i j m1 n1 m2 n2)))

행렬곱의 정의는 위와 같으므로 이를 그대로 옮긴것이 ksum 함수이다.

(defn m* [m1 m n m2 i j]
    (let [mat0 (mk-pair (arr1 1 m 1) (arr1 1 j 1))]
        (map (fn [x] (ksum n 1 0 (first x) (first (rest x)) m1 n m2 j)) mat0)))

; (m* mat1 2 2 mat2 2 2)

행렬곱을 할때는 그 행렬이 m x n 행렬임을 인자를 통해 알려줘야한다.

(def mat5 '(0 -1 2 4 11 2))
(def mat6 '(3 -1 1 2 6 1))
(m* mat5 2 3 mat6 3 2)

이제 이걸로 openGL 공부를 해보자. 행렬에 대한 이해없이도 opengl 이나 directx 프로그래밍을 할수있을지는 모른다.
하지만 대가는 될수없을것이다.

행렬 라이브러리:  지난 몇주간 아이패드와 Pixie SCHEME III 를 갖고 논 결과물이다.
버그가 버글댈것같지만, 한번 올려본다. ㅡ ㅡ

(define mat1 '(3 4 6 7))
(define mat2 '(2 3 4 5))
(define mat3 '(6 2 7   5 8 5   7 4 5))
(define mat4 '(5 3 8 4    8 5 4 9    5 7 2 1    6 4 3 8))

위와같이 행렬은 실제로는 하나의 리스트로 저장된다. opengl 에서도 쓰는 방법이고,
나 또한 3d파이프라인 구현에 사용할려고 만드는거라 그렇게 만들었다.

(define (nth arr i)
   (if (eq? i 1)
        (car arr)
   (nth (cdr arr) (- i 1))))

(define (<-m mat n i j)
   (nth mat (+ j (* (- i 1) n))))

행렬은 실제로 하나의 리스트이므로 (i,j)위치의 값을 읽기 위해서는 그 행렬이 몇 행으로 이루어져있는지를 알아야한다.

(define (arr1 a b c ret)
   (if (> a b) (reverse ret)
       (arr1 (+ a c) b c (cons a ret))))

함수프로그래밍을 위해서는 , 즉  중첩된 for 를 없애기 위해서는 먼저 (1 2 … n)  까지 만들어주는 도구가 필요

(define (m+ m1 m2 ret)
   (if (null? m1)
       (reverse ret)
       (m+ (cdr m1) (cdr m2) (cons (+ (car m1) (car m2)) ret))))

; (m+ mat1 mat2 '())

행렬의 덧셈은 정말 간단하다.

(define (m*const m1 x ret)
   (if (null? m1)
       (reverse ret)
       (m*const (cdr m1) x (cons (* (car m1) x) ret))))

;(m*const mat1 5 '())

(define (m- m1 m2)
   (m+ m1 (m*const m2 -1 '()) '()))

상수를 곱하고 행렬끼리 뺄수도 있다.  이제 행렬곱에 도전해보자.

(define (cat+ a arr ret)
   (if (null? arr)
        (reverse ret)
        (cat+ a (cdr arr) (cons (list a (car arr)) ret))))

(define (mk-pair i j ret)
   (if (null? i)
        (reverse ret)
        (mk-pair (cdr i) j (cons (cat+ (car i) j '()) ret))))

;(mk-pair '(1 2 3 4) '(1 2 3 4) '())

행렬을 표현해주는 도구도 필요하다.

(define (rm-curl src ret)
   (if (null? src)
       ret
       (rm-curl (cdr src) (append ret (car src)))))

컬리브레이스가 너무 많을때는 없애준다.

(define (ksum k c sum i j m1 n1 m2 n2)
   (if (> c k)
       sum
       (ksum k (+ c 1) (+ sum (* (<-m m1 n1 i c) (<-m m2 n2 c j))) i j m1 n1 m2 n2)))

행과 열을 서로 곱해서 모아준다.

(define (m* m1 m n m2 i j)
   (let ((mat0 (mk-pair (arr1 1 m 1 '()) (arr1 1 j 1 '()) '())))
       (let ((mat1 (rm-curl mat0 '())))
           (map (lambda (x) (ksum n 1 0 (car x) (cadr x) m1 m m2 i)) mat1))))

; (m* mat1 2 2 mat2 2 2)

이제 행렬곱을 하게 되었다.  ^ ^

(define (m-ptr m1 n i)
   (if (null? m1)
       (display "-------------------------------------------")
       (begin
           (display (car m1))
           (if (= (remainder i n) 0)
               (display "\n")
               (display "  "))

           (m-ptr (cdr m1) n (+ i 1)))))

(define (ptr1 matrix)
   (m-ptr matrix 1 1))

(define (ptr2 matrix)
   (m-ptr matrix 2 1))

(define (ptr3 matrix)
   (m-ptr matrix 3 1))

(define (ptr4 matrix)
   (m-ptr matrix 4 1))

결과를 보기 좋게(는 아니지만) 3개씩 또는 4개씩 출력해준다.

이제 행렬식 연산에 도전해본다. 먼저 여인수 부터..

; (define p1 (rm-curl (mk-pair (arr1 1 4 1 '()) (arr1 1 4 1 '()) '()) '()))

(define (cofactor m1 m n i j ret)
   (if (null? m1)
       (reverse ret)
       (cofactor (cdr m1) m n i j (if (not (or (= i (caar m1)) (= j (cadar m1))))
                                                     (cons (car m1) ret)
                                                     (cons '() ret)))))

(define p2 (cofactor p1 4 4 1 2 '()))

(4 , 4 ) 행렬에서 1 행과 2열을 제외한 여인수를 구할수있다.

(define (null-filter m1 ret)
   (if (null? m1)
       (reverse ret)
       (if (null? (car m1))
           (null-filter (cdr m1) ret)
           (null-filter (cdr m1) (cons (car m1) ret)))))

(null-filter (cofactor p1 4 4 1 2 '()) '())

널이 포함되어 있으므로 삭제해주고,

(define (peek-co m1 m n i j)
   (map (lambda (x) (<-m m1 m (car x) (cadr x))) (null-filter (cofactor
                        (rm-curl (mk-pair (arr1 1 m 1 '()) (arr1 1 n 1 '()) '()) '()) m n i j '()) '())))

(peek-co mat4 4 4 1 1)
(5 4 9 7 2 1 4 3 8)

사용하기 좋게 만들었다.

(define (det mat n)
   (if (= n 2)
       (- (* (nth mat 4) (nth mat 1)) (* (nth mat 3) (nth mat 2)))
       (apply + (map (lambda (x)
             (* (expt (- 1) (+ 1 x)) (<-m mat n 1 x) (det (peek-co mat n n 1 x) (- n 1))))
             (arr1 1 n 1 '())))))

만들고보면 별거 아닌데 , 시간을 아주 많이 잡아먹은 행렬식 연산 함수이다.  

(define mat5
'(3 -4 -5 1
  0 4 0 1
  3 1 7 1
  0 1 1 1))

(define mat6
'(3 2 -1 4
  2 1 5 7
  0 5 2 -6
 -1 2 1 0))

(det mat6 4)
-418.

음의 행렬까진 구해야하는데 , 다른 일정으로 중단 ..  ㅡ ㅡ

P402 Note: Numerical Solutions of Second-Order Differential Equations — The Runge-Kutta Method  참조
Elementary Differential Equations and Boundary Value Problems ,  8.6 Systems of First Order Equations  참조

무더운 여름, 컴퓨터를 켜지않고 아이패드에서 pixie scheme III 로 노는것이 너무 재밌다. ^ ^

그래서 전에 올렸던 너무나 Imperative 한  Runge-Kutta Method 소스를 함수프로그래밍 기법으로 다시 작성해봤다.

먼저 , P402 Note 를 보면, 단진동의 미분방정식을 어떻게 프로그래밍하는지를 가장 쉽게 알수있다.

f = ma 는  a = f/m  즉, 위의 두번째 식으로 쓸수있다. (첫번째 식은 후크의 법칙)

위의 이계도 미분방정식은 다시 두개의 first-order 방정식으로 변환된다.
미분의 정의만 알아도 위 식들을 쉽게 이해할수있다. 

P402 Note 의  C  소스를 그대로 scheme 으로 옮겨보면,

(define mass 1.0)
(define spring-constant 1.0)

(define (f time x v)
   (- (/(* spring-constant x) mass)))

(define (euler fn time x v tmax n i)
   (display time) (display " ")
   (display x) (display " ")
   (display v) (display " ")
   (newline)

   (if (> i n)
       '()
       (euler fn (+ time (/ tmax (- n 1)))
                     (+ x (* v (/ tmax (- n 1))))
                     (+ v (* (fn time x v) (/ tmax (- n 1))))
                     tmax n (+ i 1))))

;(euler f 0.0 1.0 1.0 20 100 0)

마지막의 i 는 카운터 역활이다.  런지 쿠타 방식을 적용하기 위해서는 함수의 인자수를 바꾸어야한다.
(이 부분이 현재의 난점이다)

(define (g time x)
   (- (/(* spring-constant x) mass)))

(define (rk fn time x v tmax n i)
 (begin
    (display time) (display " ")
    (display x) (display " ")
    (display v) (display " ")
    (newline)

    (if (> i n)
          '()
          (let* ((dt (/ tmax (- n 1)))
             (k1 (* dt (fn time x)))
             (k2 (* dt (fn (+ time (/ dt 2)) (+ x (/ k1 2)))))
             (k3 (* dt (fn (+ time (/ dt 2)) (+ x (/ k2 2)))))
             (k4 (* dt (fn time (+ x k3))))
             (v2 (/ (+ k1 (* 2  k2) (* 2  k3) k4) 6))
             (dx (* dt v2))
             (x2 (/ (+ dx (* 2  dx) (* 2  dx) dx) 6))
             (t2 (+ dt time)))

                (rk fn t2 x2 v2 tmax n (+ i 1))))))

(rk g 0.0 1.0 1.0 20 100 0)   와 같이 사용한다.


이제 보다 일반화 시켜서
Elementary Differential Equations and Boundary Value Problems ,  
8.6 Systems of First Order Equations 의 문제에 적용시켜보면


(define (f1 t x y)
   (- x (* 4 y)))

(define (g1 t x y)
   (- y x))

(define (euler2 f g t0 x0 y0 tmax n i)
 (begin
    (display x0) (display "    ") (display y0)
    (newline)

    (if (> i n)
          '()
          (let* ((h (/ tmax n))
             (k1 (f t0 x0 y0))
             (x1 (+ x0 (* h k1)))
             (p1 (g t0 x0 y0))
             (y1 (+ y0 (* h p1)))
             (t1 (+ h t0)))

                (euler2 f g t1 x1 y1 tmax n (+ i 1))))))

(euler2 f1 g1 0.0 1.0 0 1 10 0)

1.    0
1.1000000000000001    -0.10000000000000001
1.25    -0.22000000000000003
1.4630000000000001    -0.36699999999999999
1.7561    -0.55
2.15171    -0.78061
2.679125    -1.073842
3.3765743000000001    -1.4491387
4.2938872100000003    -1.93171
5.4959599310000007    -2.5542697209999998
7.0672638125000011    -3.3592926861999999
9.117707268230001    -4.4019483360700002

와 같이 사용할수있다.  다시 런지쿠타 를 위해 인자수를 두개로 조절한다. ㅡ ㅡ

(define (f2 x y)
   (- x (* 4 y)))

(define (g2 x y)
   (- y x))

(define (rk2 f g t0 x0 y0 tmax n i)
 (begin
    (display x0) (display "    ") (display y0)
    (newline)

    (if (> i n)
          '()
          (let* ((h (/ tmax n))
             (k1 (f x0 y0))
             (p1 (g x0 y0))
             (k2 (f (+ x0 (* h 0.5 k1)) (+ y0 (* h 0.5 p1))))
             (p2 (g (+ x0 (* h 0.5 k1)) (+ y0 (* h 0.5 p1))))
             (k3 (f (+ x0 (* h 0.5 k2)) (+ y0 (* h 0.5 p2))))
             (p3 (g (+ x0 (* h 0.5 k2)) (+ y0 (* h 0.5 p2))))
             (k4 (f (+ x0 (* h k3)) (+ y0 (* h p3))))
             (p4 (g (+ x0 (* h k3)) (+ y0 (* h p3))))
             (x1 (+ x0 (* (/ h 6) (+ k1 (* 2  k2) (* 2  k3) k4))))
             (y1 (+ y0 (* (/ h 6) (+ p1 (* 2  p2) (* 2  p3) p4))))
             (t1 (+ h t0)))

                (rk2 f g t1 x1 y1 tmax n (+ i 1))))))


(rk2 f2 g2 0 1.0 0 1 10 0)   를 실행하면,

1.    0
1.1273374999999999    -0.11125
1.3203960889062498    -0.25083259375
1.6001525300960995    -0.42966705404746086
1.9951137919483317    -0.6623967515154204
2.5439331488549146    -0.9687011072157672
3.2989432289082439    -1.3750656472659641
4.330926625352701    -1.917170503340736
5.7355568684951788    -2.6431137893803749
7.6420939775114505    -3.6177619931556868
10.225123206327096    -4.9286217159572976
13.720401496213769    -6.6937650404168991

(rk2 f2 g2 0 1.0 0 1 5 0)

1.    0
1.3200666666666667    -0.2506666666666667
1.9939111155555558    -0.6617934222222223
3.2956538046225186    -1.3734184898607409
5.727563671522403    -2.63911452154088
10.206918044252664    -4.91951640306368
18.406447393088314    -9.05262374289693

등과 같은 결과를 얻는다.  

원래 shoot 까지 구현하려고 했는데,  함수인자 갯수 부분이 완전히 이해 되지않아 다음으로 미룬다.