(define R2 (make-manifold R^n 2))

(define U (patch 'origin R2))

; these are coordinate systems
(define R2-rect (coordinate-system 'rectangular U))
(define R2-polar (coordinate-system 'polar/cylindrical U))

; these are charts, and their inverses
(define R2-rect-chi (chart R2-rect))
(define R2-rect-chi-inverse (point R2-rect))
(define R2-polar-chi (chart R2-polar))
(define R2-polar-chi-inverse (point R2-polar))

((compose R2-polar-chi R2-rect-chi-inverse) (up 'x0 'y0))
; (up 
;   (sqrt (+ (expt x0 2) (expt y0 2)))  ; radius
;   (atan y0 x0))                       ; angle (theta)

((compose R2-rect-chi R2-polar-chi-inverse) (up 'r0 'theta0))
; (up
;   (* r0 (cos theta0))     ; x
;   (* r0 (sin theta0)))    ; y


((D (compose R2-rect-chi R2-polar-chi-inverse)) (up 'r0 'theta0))
; (down
;   (up
;       (cos theta0)
;       (sin theta0))
;   (up
;       (* -1 r0 (sin theta0))
;       (* r0 (cos theta0))))

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

(define R2-rect-point (R2-rect-chi-inverse (up 'x0 'y0)))

(define corresponding-polar-point
  (R2-polar-chi-inverse
    (up (sqrt (+ (square 'x0) (square 'y0)))
        (atan 'y0 'x0))))

(f R2-rect-point)
; (f-rect (up x0 y0))

(f corresponding-polar-point)
; (f-rect (up x0 y0))

; confirms that the CAS simplifies to the same point

(define-coordinates (up x y) R2-rect)
(define-coordinates (up r theta) R2-polar)

(x (R2-rect-chi-inverse (up 'x0 'y0)))
; x0

; expect r0 * cos(theta0)
(x (R2-polar-chi-inverse (up 'r0 'theta0)))
; (* r0 (cos theta0))

; expect r0 * sin(theta0)
(y (R2-polar-chi-inverse (up 'r0 'theta0)))
; (* r0 (sin theta0))


;(r (R2-polar-chi-inverse (up 'r0 'theta0)))
; r0

; expect sqrt(x0^2 + y0^2)
(r (R2-rect-chi-inverse (up 'x0 'y0)))
; (sqrt (+ (expt x0 2) (expt y0 2)))

; expect atan(y0, x0)
(theta (R2-rect-chi-inverse (up 'x0 'y0)))
; (atan y0 x0)

; h: x * r^2 + y^3
(define h (+ (* x (square r)) (cube y)))

(h R2-rect-point)
; (+ (expt x0 3) (* x0 (expt y0 2)) (expt y0 3))
; aka: x0^3 + x0 y0^2 + y0^3 
;      x0 (x0^2 + y0^2) + y0^3

(h (R2-polar-chi-inverse (up 'r0 'theta0)))
; (+ (* (expt r0 3) (expt (sin theta0) 3)) (* (expt r0 3) (cos theta0)))


;;; Exersize 2.1a

((- r (* 2 'a (+ 1 (cos theta))))
 ((point R2-rect) (up 'x 'y)))
; (/ 
;   (+ (* -2 a x)
;      (* -2 a (sqrt (+ (expt x 2) (expt y 2))))
;      (expt x 2) (expt y 2))
;   (sqrt (+ (expt x 2) (expt y 2))))


; "Lemniscate of Bernoulli"
; (x^2 + y^2)^2 = 2 a^2 (x^2 - y^2)

((-
   (square (+ (square x) (square y)))
   (* 2
      (square 'a)
      (- (square x) (square y))))
 ((point R2-polar) (up 'r0 'theta0)))
; (+ (* -4 (expt a 2) (expt r0 2) (expt (cos theta0) 2))
;    (* 2 (expt a 2) (expt r0 2))
;    (expt r0 4))
;
; r^2 / a^2 + 2 = 4 cos^2 (theta)
; => r^2 = 2 a^2 cos(2 theta)

; this matches Wikipedia, with substitution of a^2 = 2 c^2 (in the rectangular version)


;;; Exersize 2.1b

(define R3-rect-chi (chart R3-rect))
(define R3-rect-chi-inverse (point R3-rect))
(define R3-cyl-chi (chart R3-cyl))
(define R3-cyl-chi-inverse (point R3-cyl))

(define-coordinates (up x y z) R3-rect)
(define-coordinates (up r theta z-cyl) R3-cyl)

; R->R3
(define (helix-rect t)
  (up (* 'R (cos t))
      (* 'R (sin t))
      (* 's t)))

; R->R3
(define (helix-cyl t)
  (up 'R
      t
      (* 's t)))

((compose R3-rect-chi R3-cyl-chi-inverse helix-cyl) 't)
; (up
;   (* R (cos t))
;   (* R (sin t))
;   (* s t))


((compose R3-cyl-chi R3-rect-chi-inverse helix-rect) 't)
; (up R
;     t
;     (* s t))

((- helix-rect (compose R3-rect-chi R3-cyl-chi-inverse helix-cyl)) 't)
; (up 0 0 0)

;;; Exersize 2.2

; given polar coordinates of point on plane, get to spherical coordinates of point on sphere with:
((compose
   (chart S2-spherical)
   (point S2-Riemann)
   (chart R2-rect)
   (point R2-polar))
 (up 'rho 'theta))
; (up (acos (/ (+ -1 (expt rho 2))
;              (+ 1 (expt rho 2))))
;     theta)

; given spherical coordinates of point on sphere, what are the polar
; coordinates of corresponding point on the plane?
((compose
   (chart R2-polar)
   (point R2-rect)
   (chart S2-Riemann)
   (point S2-spherical))
 (up 'phi 'lambda))
; (up (/ (sin phi)
;        (+ -1 (cos phi)))
;     (atan (* -1 (sin lambda))
;           (* -1 (cos lambda))))

; this is just the inverse of the above, right?

; huh, I would expect the second term to just be lambda.
; atan (-sin(lambda), -cos(lambda)) -> tan^-1(tan(lambda)) -> lambda
; seems like an identity, ok