(load "incdec.sc")

(define first car)

(define rest cdr)

(define (map-n f n)
 (let loop ((i 0)) (if (= i n) '() (cons (f i) (loop (+ i 1))))))

(define (reduce f l i) (if (null? l) i (f (first l) (reduce f (rest l) i))))

(define (sqr x) (* x x))

(define (v+ u v) (map + u v))

(define (v- u v) (map - u v))

(define (k*v k v) (map (lambda (x) (* k x)) v))

(define (dot u v) (reduce + (map * u v) 0))

(define (distance u v) (let ((d (v- v u))) (sqrt (dot d d))))

(define (replace-ith x i xi)
 (if (zero? i)
     (cons xi (rest x))
     (cons (first x) (replace-ith (rest x) (- i 1) xi))))

(define (gradient f)
 (lambda (x)
  (map-n (lambda (i)
	  ((derivative (lambda (xi) (f (replace-ith x i xi)))) (list-ref x i)))
	 (length x))))

(define *print?* #f)
(define x-initial '(0 8))
(define xdot-initial '(0.75 0))
(define w0 0)
(define error-tolerance 1e-1)
(define delta-t 1e-1)

(define (naive-euler w)
 (let ((charges (list (list 10 (- 10 w)) (list 10 0))))
  (define (p x) (reduce + (map (lambda (c) (/ 1 (distance x c))) charges) 0))
  (let loop ((x x-initial) (xdot xdot-initial))
   (if *print?*
       (begin (write (list-ref x 0))
	      (display " ")
	      (write (list-ref x 1))
	      (newline)))
   (let* ((xddot (k*v -1 ((gradient p) x)))
	  (x-new (v+ x (k*v delta-t xdot))))
    (if (positive? (list-ref x-new 1))
	(loop x-new (v+ xdot (k*v delta-t xddot)))
	(let* ((delta-t-f (/ (- 0 (list-ref x 1)) (list-ref xdot 1)))
	       (x-t-f (v+ x (k*v delta-t-f xdot))))
	 (sqr (list-ref x-t-f 0))))))))

(define (path w)
 (let ((print? *print?*))
  (set! *print?* #t)
  (let ((e (naive-euler w))) (set! *print?* print?) e)))

(define (truncate-number x)
 (let ((s (number->string x))) (substring s 0 (min (string-length s) 6))))

(define (argmin-using-textbook-newtons-method f x)
 (call-with-output-file "particle.plt"
  (lambda (output-port)
   (display "set output 'particle.eps'" output-port)
   (newline output-port)
   (display "set size square" output-port)
   (newline output-port)
   (display "set terminal postscript eps enhanced monochrome blacktext \"Times-Roman\" 14"
	    output-port)
   (newline output-port)
   (display "set xrange [0:10]" output-port)
   (newline output-port)
   (display "set yrange [0:10]" output-port)
   (newline output-port)
   (display "set style data line" output-port)
   (newline output-port)
   (display "set key on right bottom Left reverse spacing 1" output-port)
   (newline output-port)
   (display "set title 'Path of Charged Particle'" output-port)
   (newline output-port)
   (display "plot " output-port)
   (let loop ((x x) (i 0))
    (let ((path-filename
	   (string-append "particle-" (number->string i) ".data")))
     (write path-filename output-port)
     (display " title \"w^{(" output-port)
     (write i output-port)
     (display ")}=" output-port)
     (display (truncate-number x) output-port)
     (display "\"" output-port)
     (let ((df-dx ((derivative f) x)))
      (if (< (abs df-dx) error-tolerance)
	  (newline output-port)
	  (display ", " output-port))
      (write (list x (f x) df-dx))
      (newline)
      (with-output-to-file path-filename (lambda () (path x)))
      (if (< (abs df-dx) error-tolerance)
	  x
	  (loop (- x (/ df-dx ((derivative (derivative f)) x))) (+ i 1)))))))))

(define (particle) (argmin-using-textbook-newtons-method naive-euler w0))