;;; 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->vecotr-field components coordinate-system)
  (procedure->vector-field
   (vector-field-procedure components coordinate-system)))

(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))))

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

(define ((form-field-over-map mu:N->M) w-on-M)
  (let ((k (get-rank w-on-M)))
    (procedure->nform-field
     (lambda vecotrs-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 sphere-basis-over-mu
  (basis->basis-over-map mu sphere-basis))