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

#|

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

|#