;;; Jan 29th

(load "~/thesis/scmutils/src/calculus/load.scm")

;(define R2 (make-manifold R^n-type 2)) ; doesn't work, so...
(define R2 (rectangular 2))
;(define U (patch 'origin R2)) ; also undefined so going to AIM-2005-003.pdf

(define R2 (rectangular 2))
(define P2 (polar/cylindrical 2))


(define R2-chi-inverse (R2 '->point))
(define R2-chi (R2 '->coords))
(define P2-chi (P2 '->coords))
(define P2-chi-inverse (P2 '->point))

;;; Feb 1st

(print-expression
 ((compose R2-chi-inverse P2-chi)
  (up 'x0 'y0)))
;= (up (sqrt (+ (expt x0 2) (expt y0 2))) (atan y0 x0))

(print-expression
 ((compose R2-chi P2-chi-inverse)
  (up 'r0 'theta0)))

; (up (* r0 (cos theta0)) (* r0 (sin theta0)))

(define R2->R (-> (UP Real Real) Real))

(define f
  (compose (literal-function 'f-rect R2->R)
	   R2-chi))
; there is a simpler syntax...

(define R2-point (R2-chi-inverse (up 'x0 'y0)))
(define P2-point (P2-chi-inverse (up 'r0 'theta0)))

(print-expression (f R2-point))
; (f-rect (up x0 y0))
(print-expression (f P2-point))
; (f-rect (up (* r0 (cos theta0)) (* r0 (sin theta0))))

(define g (literal-manifold-function 'g-polar P2))

(print-expression (g R2-point))
; (g-polar (up (sqrt (+ (expt x0 2) (expt y0 2))) (atan y0 x0)))

(instantiate-coordinates R2 '(x y))
(instantiate-coordinates P2 '(r theta))

(print-expression (x R2-point))
;x0

(print-expression (theta R2-point))
;(atan y0 x0)

(define h (+ (* x (square r)) (cube y)))

(print-expression (h P2-point))
; (+ (* (expt r0 3) (expt (sin theta0) 3)) (* (expt r0 3) (cos theta0)))

(print-expression (h R2-point))
; (+ (expt x0 3) (* x0 (expt y0 2)) (expt y0 3))

(print-expression
 (D h))
; a-euclidean-derivative

; pe is print-expression

(pe ((D h) R2-point))
; (down (+ (* 3 (expt x0 2)) (expt y0 2)) (+ (* 2 x0 y0) (* 3 (expt y0 2))))

;(pe (D (h R2-point)))
; ERROR

(pe ((D (compose h R2-chi-inverse)) R2-point))
; (down (+ (* 3 (expt x0 2)) (expt y0 2)) (+ (* 2 x0 y0) (* 3 (expt y0 2))))

; TODO: formal definition of operator vs. function?

;;; Feb 2

#|
(define (vector-field-procedure components coordinate-system)
  (define (the-procedure f)
    (compose (* (D (compose f (coordinate-system '->point)))
		components)
	     (coordinate-system '->coords)))
  the-procedure)

(define (components->vector-field components coordinate-system)
  (procedure->vector-field
   (vector-field-procedure components coordinate-system)
   'foo))

(define v
  (components->vector-field
   (up (literal-function 'vx (-> (UP Real Real) Real))
       (literal-function 'vy (-> (UP Real Real) Real)))
   R2))

; shorthand:
; (define w (literal-vector-field 'v R2))

;(pe ((v (literal-manifold-function 'f R2)) R2-point))

(define (coordinatize vector-field coordsys)
  (define ((v f) x)
    (let ((b (compose (vector-field (coordsys '->coords))
		      (coordsys '->point))))
      (* ((D f) x) (b x))))
  (make-operator v))
|#

#|
(pe (((coordinatize v R2)
      (literal-function 'f (-> (UP Real Real) Real)))
     (up 'x0 'y0)))
|#

; Note: nice summary of vector field properties on p12

;(pe ((d/dx (square r)) R2-point))
;(pe ((d/dx (square r)) P2-point))

(define J (- (* x d/dy) (* y d/dx)))

(series:for-each pe 
		 (((exp (* 'a J)) R2-chi)
		  ((R2 '->point) (up 1 0)))
		 6)
;(up 1 0)
;(up 0 a)
;(up (* -1/2 (expt a 2)) 0)
;(up 0 (* -1/6 (expt a 3)))
;(up (* 1/24 (expt a 4)) 0)
;(up 0 (* 1/120 (expt a 5)))

; now do evolution on coordinates

(define ((((evolution order)
	   delta-t vector-field)
	  manifold-function)
	 manifold-point)
  (series:sum
   (((exp (* delta-t vector-field))
     manifold-function)
    manifold-point)
   order))

(pe ((((evolution 6) 'a J) R2-chi)
     ((R2 '->point) (up 1 0))))

(pe ((((evolution 6) 2. J) R2-chi)
     ((R2 '->point) (up 1 0))))

#|
; super chunky, i would rewrite these for sure
(define mywindow (frame -4 4 -4 4))

(define (plot-evolution win order initial a step)
  (letrec 
      ((dostep 
	(lambda (val)
	  (cond
	   ((< val a) (let ((this-point ((((evolution order) val J) R2-chi)
					 ((R2 '->point) initial))))
			(plot-point win
				    (time this-point)
				    (coordinate this-point))
			(dostep (+ val step))))))))
    (dostep 0.)))

(plot-evolution mywindow 6 (up 1. 0.) 6. .01)



(define (explore-evolution window order length)
  (define (iterate-mything i x y)
    (if (< i length)
	(let ((this-point ((((evolution order) i J) R2-chi)
			   ((R2 '->point) (up x y)))))
	  (plot-point window (time this-point) (coordinate this-point))
	  (iterate-mything (+ .01 i) x y))
	(button-loop x y)))
  (define (button-loop ox oy)
    (pointer-coordinates
     window
     (lambda (x y button)
       (let ((temp button))
         (cond ((eq? temp 0) (write-line (cons x (cons y (quote ()))))
                             (display " started.")
                             (iterate-mything 0 x y))
               ((eq? temp 1) (write-line (cons ox (cons oy (quote ()))))
                             (display " continued.")
                             (iterate-mything 0 ox oy))
               ((eq? temp 2) (write-line (cons x (cons y (quote ()))))
                             (display " hit.")
                             (button-loop ox oy)))))))
  (newline)
  (display "Left button starts a trajectory.")
  (newline)
  (display "Middle button continues a trajectory.")
  (newline)
  (display "Right button interrogates coordinates.")
  (button-loop 0. 0.))

(explore-evolution mywindow 5 .4)

|#

(define R3 (rectangular 3))
(instantiate-coordinates R3 '(x y z))
(define R3-chi (R3 '->coords))
(define R3-chi-inverse (R3 '->point))
(define R3->R (-> (UP Real Real Real) Real))
(define R3-point (R3-chi-inverse (up 'x0 'y0 'z0)))

(define omega (literal-1form-field 'omega R3))

(define v (literal-vector-field 'v R3))
(define w (literal-vector-field 'w R3))
(define c (literal-manifold-function 'c R3))

(pe ((- (omega (+ v w)) (+ (omega v) (omega w))) R3-point))
;0

(pe ((- (omega (* c v)) (* c (omega v))) R3-point))
;0

(define omega
  (components->1form-field
   (down (literal-manifold-function 'a_0 R3)
	 (literal-manifold-function 'a_1 R3)
	 (literal-manifold-function 'a_2 R3))
   R3))

(pe ((dx d/dx) R3-point))
; 1

(pe ((omega (literal-vector-field 'v R3)) R3-point))

; back to R2?
(instantiate-coordinates R2 '(x y))
(pe ((dy J) R2-point))

(define e0
  (+ (* (literal-manifold-function 'e0x R2) d/dx)
     (* (literal-manifold-function 'e0y R2) d/dy)))

(define e1
  (+ (* (literal-manifold-function 'e1x R2) d/dx)
     (* (literal-manifold-function 'e2x R2) d/dy)))

(define e-vector-basis (down e0 e1))
(define e-dual-basis (vector-basis->dual e-vector-basis P2))

;(pe ((e-dual-basis e-vector-basis) m))
; TODO: m undefined

(define l-m-f literal-manifold-function) ; type this a lot!

(define v
  (* e-vector-basis
     (up (l-m-f 'bx R2)
	 (l-m-f 'by R2))))

#|
(pe ((e-dual-basis v) R2-point))
; takes a super long time... still have R3 definitions?

(let ((polar-vector-basis (basis->vector-basis polar-basis))
      (polar-dual-basis (basis->1form-basis polar-basis)))
  (pe ((- (commutator e0 e1 f)
	  (* (- (e0 (polar-dual-basis e1))
		(e1 (polar-dual-basis e0)))
	     (polar-vector-basis f)))
       R2-point)))

; TODO: missing defs... goddamn it's cold

|#

; back to R3
(define R3 (rectangular 3))
(instantiate-coordinates R3 '(x y z))
(define R3->R (-> (UP Real Real Real) Real))
(define g 
  (compose (literal-function 'g R3->R)
	   (R3 '->coords)))
(define R3-point ((R3 '->point) (up 'x 'y 'z)))

(define Jz (- (* x d/dy) (* y d/dx)))
(define Jx (- (* y d/dz) (* z d/dy)))
(define Jy (- (* z d/dx) (* x d/dz)))

(pe (((+ (commutator Jx Jy) Jz) g) R3-point))
; 0
; so [Jx,Jy] = -Jz, nice

;;; Feb 4

(define setup-R2
  (lambda ()
    (let ()
      (define R2 (rectangular 2))
      (instantiate-coordinates R2 '(x y))
      (define R2-point ((R2 '->point) (up 'x0 'y0))))))

(define setup-R3
  (lambda ()
    (let ()
      (define R3 (rectangular 3))
      (instantiate-coordinates R3 '(x y z))
      (define R3-point ((R3 '->point) (up 'x0 'y0 'z0))))))

(setup-R2)

(define v (+ (* 'a d/dx) (* 'b d/dy)))
(define w (+ (* 'c d/dx) (* 'd d/dy)))

;(pe (((wedge dx dy) v w) R2-point))

(setup-R3)

(define u (+ (* 'a d/dx) (* 'b d/dy) (* 'c d/dz)))
(define v (+ (* 'd d/dx) (* 'e d/dy) (* 'f d/dz)))
(define w (+ (* 'g d/dx) (* 'h d/dy) (* 'i d/dz)))

;(pe (((wedge dx dy dz) u v w) R3-point))
; determinant

(define a (l-m-f 'alpha R3))
(define b (l-m-f 'beta R3))
(define c (l-m-f 'gamma R3))

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

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

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

|#

(define ((vector-field-over-map mu:N->M) v-on-M)
  (procedure->vector-field
   (lambda (f-on-M)
     (compose (v-on-M f-on-M) mu:N->M))
   'foo))

(define (vector-field-over-map->vector-field V-over-mu n)
  (procedure->vector-field
   (lambda (f)
     (lambda (m) ((V-over-mu f) n)))
   'foo))

(define ((form-field-over-map mu:N->M) w-on-M)
  (let ((k (get-rank w-on-M)))
    (procedure->nform-field
     (lambda vectors-over-map
       (lambda (n)
	 ((apply w-on-M
		 (map (lambda (V-over-mu)
			(vector-field-over-map->vector-field
			 V-over-mu n))
		      vectors-over-map))
	  (mu:N->M n))))
     'athing
     k)))

(pp procedure->nform-field)

(define sphere (S2 1))
(instantiate-coordinates sphere '(theta phi))
(define sphere-basis (coordinate-system->basis sphere))
(instantiate-coordinates the-real-line 't)

(define mu
  (compose (sphere '->point)
	   (up (literal-function 'alpha)
	       (literal-function 'beta))
	   (the-real-line '->coords)))

(define (((differential mu) v) f)
  (v (compose f mu)))

;;; Feb 5
;;; skipped a bit here...

(define setup-P2
  (lambda ()
    (define P2 (polar/cylindrical 2))
    (define P2-chi (P2 '->coords))
    (define P2-chi-inverse (P2 '->point))
    (define polar-basis (coordinate-system->basis P2))
    (instantiate-coordinates P2 '(r theta))))

(setup-R2)
(define rectangular-basis (coordinate-system->basis R2))
(setup-P2)
(define polar-basis (coordinate-system->basis P2))

(define rectangular-Christoffel
  (make-Christoffel
   (let ((zero (lambda (m) 0)))
     (down (down (up zero zero)
		 (up zero zero))
	   (down (up zero zero)
		 (up zero zero))))
   rectangular-basis))

(define polar-Christoffel
  (make-Christoffel
   (let ((zero (lambda (m) 0)))
     (down (down (up zero zero)
		 (up zero (/ 1 r)))
	   (down (up zero (/ 1 r))
		 (up (* -1 r) zero))))
   polar-basis))

(define rectangular-Cartan
  (Christoffel->Cartan rectangular-Christoffel))

(define polar-Cartan
  (Christoffel->Cartan polar-Christoffel))

(define f
  (compose (literal-function 'f-rect R2->R) R2-chi))

#|
(pe (((((covariant-derivative rectangular-Cartan) d/dx)
       J)
      f)
     R2-point))

(pe ((d/dy f) R2-point))
|#

(instantiate-coordinates the-real-line 't)
(define M (rectangular 2))
(instantiate-coordinates M '(theta phi))
(define M-basis (coordinate-system->basis M))

(define G-S2-1
  (make-Christoffel
   (let ((zero (lambda (point) 0)))
     (down (down (up zero zero)
		 (up zero (/ 1 (tan theta))))
	   (down (up zero (/ 1 (tan theta)))
		 (up (- (+ (sin theta) (cos theta))) zero))))
   M-basis))

(define gamma:N->M
  (compose (M '->point)
	   (up (literal-function 'alpha)
	       (literal-function 'beta))
	   (the-real-line '->coords)))

(define basis-over-gamma
  (basis->basis-over-map gamma:N->M M-basis))
; this gave me an error; i tried from a blank run and didn't get the error so 
; I must have clobbered something in this file

; ah, it was easy, i /was/ passing the wrong arguments to 
; procedure->vector-field

(define w
  (basis-components->vector-field
   (up (compose (literal-function 'w0)
		(the-real-line '->coords))
       (compose (literal-function 'w1)
		(the-real-line '->coords)))
   (basis->vector-basis basis-over-gamma)))

(define sphere-Cartan-over-gamma
  (Christoffel->Cartan-over-map G-S2-1 gamma:N->M))

#|
(pe
 (s:map/r
  (lambda (omega)
    ((omega
      (((covariant-derivative sphere-Cartan-over-gamma)
	d/dt)
       w))
     ((the-real-line '->point) 'tau)))
  (basis->1form-basis basis-over-gamma)))

|#

(define gamma:N->M
  (compose (M '->point)
	   (up (literal-function 'alpha)
	       (literal-function 'beta))
	   (the-real-line '->coords)))

#|
(se
 (((((covariant-derivative sphere-Cartan-over-gamma) d/dt)
    ((differential gamma:N->M) d/dt))
   (M '->coords))
  ((the-real-line '->point) 't)))
|#

(define (Lfree s)
  (* 1/2 (square (velocity s))))

(define (sphere->R3 s)
  (let ((q (coordinate s)))
    (let ((theta (ref q 0)) (phi (ref q 1)))
      (up (* (sin theta) (cos phi))
	  (* (sin theta) (sin phi))
	  (cos theta)))))

(define Lsphere
  (compose Lfree (F->C sphere->R3)))

#|
(pe (((Lagrange-equations Lsphere)
      (up (literal-function 'alpha)
	  (literal-function 'beta)))
     't))
|#

(define ((Riemann-curvature Cartan) u v)
  (let ((nabla (covariant-derivative Cartan)))
    (- (commutator (nabla u) (nabla v))
       (nabla (commutator u v)))))

(define sphere-Cartan (Christoffel->Cartan G-S2-1))

(pe (((Riemann sphere-Cartan) dphi d/dtheta d/dphi d/dtheta)
     ((M '->point) (up 'theta0 'phi0))))
;1

(define SR (rectangular 4))
(instantiate-coordinates SR '(t x y z))
(define an-event ((SR '->point) (up 't0 'x0 'y0 'z0)))
(define c 'c)

(define (g-Lorentz u v)
  (+ (* (dx u) (dx v))
     (* (dy u) (dy v))
     (* (dz u) (dz v))
     (* -1 (square c) (dt u) (dt v))))

(define SR-vector-basis
  (down (* (/ 1 c) d/dt)
	d/dx
	d/dy
	d/dz))

(define SR-1form-basis
  (up (* c dt) dx dy dz))

(define SR-basis
  (make-basis SR-vector-basis SR-1form-basis))

(define (Faraday Ex Ey Ez Bx By Bz)
  (+ (* Ex c (wedge dx dt))
     (* Ey c (wedge dy dt))
     (* Ez c (wedge dz dt))
     (* Bx (wedge dy dz))
     (* By (wedge dz dx))
     (* Bz (wedge dx dy))))

(define (Maxwell Ex Ey Ez Bx By Bz)
  (+ (* -1 Bx c (wedge dx dt))
     (* -1 By c (wedge dy dt))
     (* -1 Bz c (wedge dz dt))
     (* Ex (wedge dy dz))
     (* Ey (wedge dz dx))
     (* Ez (wedge dx dy))))

; blah doesn't interpret

(define SR-star
  (Hodge-star g-Lorentz SR-basis))

(SR-star (Faraday 'Ex 'Ey 'Ez 'Bx 'By 'Bz))