; based off the development in ch6.2 of Reflections on Relativity at
; mathpages.com
;

;========================================================================
; first prepare a schwartzchild
; definitions are in calculus/manifold.scm

#| UNUSED see orbital-plane below

(install-coordinates spacetime-rect (up 't 'x 'y 'z))
(pp (spacetime-rect 'coordinate-prototype))
; #(t x y z)

(install-coordinates spacetime-spher (up 't 'r 'theta 'phi))
(pp (spacetime-spher 'coordinate-prototype))
; #(t r phi theta)

(spacetime-spher 'check-coordinates)

(define dT (1form-field->tensor-field dt))

(define dX (1form-field->tensor-field dx))
(define dY (1form-field->tensor-field dy))
(define dZ (1form-field->tensor-field dz))

(define dR (1form-field->tensor-field dr))
(define dTheta (1form-field->tensor-field dtheta))
(define dPhi (1form-field->tensor-field dphi))

(define rect-event ((spacetime-rect '->point) (up 't0 'x0 'y0 'z0)))
(define spher-event ((spacetime-rect '->point) (up 't0 'r0 'theta0 'phi0)))
|#

; for orbital plane: (tau,r,phi)
(attach-coordinate-system 'spacetime-polar 'origin R^n-type
  (lambda (manifold)
    (define (me m)
      (case m
    ((check-coords)
     (lambda (coords)
       (and (up? coords)
        (fix:= (s:dimension coords) 3))))
    ((coords->point)
     (lambda (coords)
       (if (and (up? coords)
		(fix:= (s:dimension coords) 3))
           (let ((t (ref coords 0))
		 (r (ref coords 1))
		 (phi   (ref coords 2)))
	     (make-manifold-point
	      (up t
		  (* r (cos phi))
		  (* r (sin phi)))
	      manifold
	      me
	      coords))
           (error "Bad coordinates: spacetime-polar"
		  coords))))
    ((check-point)
     (lambda (point)
       (my-manifold-point? point manifold)))
    ((point->coords)
     (lambda (point)
       (assert ((me 'check-point) point) "Bad point assertion: spacetime-polar")
       (get-coordinates point me
         (lambda ()
           (let ((prep (manifold-point-representation point)))
         (if (and (up? prep) (fix:= (s:dimension prep) 3))
             (let ((t (ref prep 0))
		   (x (ref prep 1))
		   (y (ref prep 2)))
               (let ((r (sqrt (+ (square x) (square y)))))
		 (up t
		     r
		     (atan y x))))
             (error "Bad point: spacetime-polar" point)))))))
    (else (error "Bad message: coordinate-transformer" m))))
    me))

(define orbital-plane (make-manifold R^n-type 3))
(define orbital-plane-rect (coordinate-system-at 'rectangular 'origin R3))
(define orbital-plane-polar (coordinate-system-at 'spacetime-polar 'origin R3))

(install-coordinates orbital-plane-rect (up 'tau 'x 'y))
(install-coordinates orbital-plane-polar (up 'tau 'r 'phi))

(install-coordinates the-real-line (up 'k))

(define dT (1form-field->tensor-field dtau))

(define dX (1form-field->tensor-field dx))
(define dY (1form-field->tensor-field dy))

(define dR (1form-field->tensor-field dr))
(define dPhi (1form-field->tensor-field dphi))

(define rect-event ((orbital-plane-rect '->point) (up 'tau0 'x0 'y0)))
(define polar-event ((orbital-plane-polar '->point) (up 'tau0 'r0 'phi0)))

(define (Schwartzchild-plane-metric M)
  (+ (* (/ (- r (* 2 M)) r) (square dT))
     (* -1 (/ r (- r (* 2 M))) (square dR))
     (* -1 (square r) (square dPhi))))

(define S-p-metric (Schwartzchild-plane-metric 'M))

;(pe (coordinate-system->basis orbital-plane-polar))

(define Gamma
  (Christoffel S-p-metric orbital-plane-polar))

#|

(show-time 
 (lambda ()
   (pec ((tensor-field->coefficient-structure Gamma orbital-plane-polar)
	 polar-event))))
|#

(begin
  (define orbit-Cs-2
    (metric->Christoffel-2 S-p-metric 
			   (coordinate-system->basis orbital-plane-polar)))
  (define orbit-Cs-1
    (metric->Christoffel-1 S-p-metric 
			   (coordinate-system->basis orbital-plane-polar)))
  (se ((Christoffel->symbols orbit-Cs-2)
       ((orbital-plane-polar '->point) (up 'tau0 'r0 'phi0)))))

(install-coordinates the-real-line (up 'k))

(define test-path
  (compose (orbital-plane-polar '->point)
	   (up (literal-function 'alpha)
	       (literal-function 'beta)
	       (literal-function 'gamma))
	   (the-real-line '->coords)))

;(se ((orbital-plane-polar '->coords) (test-path ((the-real-line '->point) 'k))))

;((orbital-plane-polar '->point)  ((orbital-plane-polar '->coords) (test-path ((the-real-line '->point) 'k))))

(define polar-Cartan-over-path
  (Christoffel->Cartan-over-map orbit-Cs-1 test-path))

(((((covariant-derivative polar-Cartan-over-path) d/dk)
   ((differential test-path) d/dk))
  (orbital-plane-polar '->coords))
 ((the-real-line '->point) 'k))


(((((covariant-derivative polar-Cartan-over-path) d/dk)
   ((differential test-path) d/dk))
  (orbital-plane-polar '->coords))
 ((the-real-line '->point) 'k))


(se ((orbital-plane-polar '->coords) (test-path ((the-real-line '->point) 't))))

(define polar-Cartan (Christoffel->Cartan orbit-Cs-1))




(define basis-over-path
  (basis->basis-over-map test-path 
			 (coordinate-system->basis orbital-plane-polar)))

;(basis-over-path (test-path ((the-real-line '->point) 'k)))

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

(pe
 (s:map/r
  (lambda (omega)
    ((omega
      (((covariant-derivative polar-Cartan-over-path)
	d/dt)
       w))
     ((the-real-line '->point) 'k)))
  (basis->1form-basis basis-over-path)))

(((basis->1form-basis basis-over-path)
  (((covariant-derivative polar-Cartan-over-path) d/dk)
   ((differential test-path) d/dk))
  (orbital-plane-polar '->coords))
 ((the-real-line '->point) 'k))







;===========================================================
; for mercury, [TODO: cite these]
; precesion:                ~43.92 arcsec/century
;                           ~0.103 arcsec/rev
; eccentricity k:           0.2060
; semimajor axis:           57.9 x 10^6 km