; deps
(define R2->R (-> (UP Real Real) Real))
(define R2-rect-chi (chart R2-rect))
(define R2-rect-chi-inverse (point R2-rect))
(define R2-rect-point (R2-rect-chi-inverse (up 'x0 'y0)))

(define R3->R (-> (UP Real Real Real) Real))
(define R3-rect-chi (chart R3-rect))
(define R3-rect-chi-inverse (point R3-rect))
(define R3-rect-point (R3-rect-chi-inverse (up 'x0 'y0 'z0)))

(define-coordinates (up x y z) R3-rect)

; two regular old vectors
(define u (+ (* 'u^0 d/dx) (* 'u^1 d/dy)))
(define v (+ (* 'v^0 d/dx) (* 'v^1 d/dy)))

(((wedge dx dy) u v) R3-rect-point)
; (+ (* u^0 v^1)
;    (* -1 u^1 v^0))

(define-coordinates (up r theta z) R3-cyl)

; two more vectors
(define a (+ (* 'a^0 d/dr) (* 'a^1 d/dtheta)))
(define b (+ (* 'b^0 d/dr) (* 'b^1 d/dtheta)))

(((wedge dr dtheta) a b) ((point R3-cyl) (up 'r0 'theta0 'z0)))
; (+ (* a^0 b^1) (* -1 a^1 b^0))

; the same thing, as expected

; fancy verification of determinant
; note that we are re-defining u,v

(define u (+ (* 'u^0 d/dx) (* 'u^1 d/dy) (* 'u^2 d/dz)))
(define v (+ (* 'v^0 d/dx) (* 'v^1 d/dy) (* 'v^2 d/dz)))
(define w (+ (* 'w^0 d/dx) (* 'w^1 d/dy) (* 'w^2 d/dz)))

(((wedge dx dy dz) u v w) R3-rect-point)
; (+ (* u^0 v^1 w^2)
;    (* -1 u^0 v^2 w^1)
;    (* -1 u^1 v^0 w^2)
;    (* u^1 v^2 w^0)
;    (* u^2 v^0 w^1)
;    (* -1 u^2 v^1 w^0))

(- (((wedge dx dy dz) u v w) R3-rect-point)
   (determinant
     (matrix-by-rows (list 'u^0 'u^1 'u^2)
                     (list 'v^0 'v^1 'v^2)
                     (list 'w^0 'w^1 'w^2))))
; 0

; continuing with bits after exersize 5.1...

(define a (literal-manifold-function 'alpha R3-rect))
(define b (literal-manifold-function 'beta R3-rect))
(define c (literal-manifold-function 'gamma R3-rect))

(define theta (+ (* a dx) (* b dy) (* c dz)))

(define X (literal-vector-field 'X-rect R3-rect))
(define Y (literal-vector-field 'Y-rect R3-rect))

(((- (d theta)
     (+ (wedge (d a) dx)
        (wedge (d b) dy)
        (wedge (d c) dz)))
  X Y)
 R3-rect-point)
; 0 

(define omega
  (+ (* a (wedge dy dz))
     (* b (wedge dz dx))
     (* c (wedge dx dy))))

(define Z (literal-vector-field 'Z-rect R3-rect))

(((- (d omega)
     (+ (wedge (d a) dy dz)
        (wedge (d b) dz dx)
        (wedge (d c) dx dy)))
  X Y Z)
 R3-rect-point)
; 0

(define v (literal-vector-field 'v-rect R2-rect))
(define w (literal-vector-field 'w-rect R2-rect))

(define R2-rect-basis (coordinate-system->basis R2-rect))

(define alpha (literal-function 'alpha R2->R))
(define beta (literal-function 'beta R2->R))

(let ((dx (ref (basis->1form-basis R2-rect-basis) 0))
      (dy (ref (basis->1form-basis R2-rect-basis) 1)))
  (((- (d (+ (* (compose alpha (chart R2-rect)) dx)
             (* (compose beta (chart R2-rect)) dy)))
       (* (compose (- ((partial 0) beta)
                      ((partial 1) alpha))
                   (chart R2-rect))
          (wedge dx dy)))
    v w)
   R2-rect-point))
; 0

(define a (literal-manifold-function 'a-rect R3-rect))
(define b (literal-manifold-function 'b-rect R3-rect))
(define c (literal-manifold-function 'c-rect R3-rect))

(define flux-through-boundary-element
  (+ (* a (wedge dy dz))
     (* b (wedge dz dx))
     (* c (wedge dx dy))))

(define production-in-volume-element
  (* (+ (d/dx a) (d/dy b) (d/dz c))
     (wedge dx dy dz)))

(define X (literal-vector-field 'X-rect R3-rect))
(define Y (literal-vector-field 'Y-rect R3-rect))
(define Z (literal-vector-field 'Z-rect R3-rect))

(((- production-in-volume-element
     (d flux-through-boundary-element))
  X Y Z)
 R3-rect-point)
; 0

;;; Exercise 5.1: Wedge Product

; let's do R2-polar! small and simple

(define-coordinates (up r theta) R2-polar)
(define R2-polar-chi (chart R2-polar))
(define R2-polar-chi-inverse (point R2-polar))
(define R2-polar-point (R2-polar-chi-inverse (up 'r0 'theta0)))

; a) verify that the wedge product is associative for forms in R2-polar

; A wedge ( B wedge C ) == ( A wedge B) wedge C

(define A (literal-1form-field 'a R2-polar))
(define B (literal-1form-field 'b R2-polar))
(define C (literal-1form-field 'c R2-polar))

(define m (+ (* 'm^0 d/dr) (* 'm^1 d/dtheta)))
(define n (+ (* 'n^0 d/dr) (* 'n^1 d/dtheta)))
(define p (+ (* 'p^0 d/dr) (* 'p^1 d/dtheta)))


(((- (wedge (wedge A B) C)
     (wedge A (wedge B C))) m n p) R2-polar-point)
; 0

; b) verify that formula (5.17) is true in R2-polar
; (dx wedge dy wedge ...) (d/dx, d/dy, ...) = 1

(((wedge dr dtheta) d/dr d/dtheta) R2-polar-point)
; 1

; exersize 5.2 and 5.3 would be in notebook, but not planning to do these right
; now