(use random-bsd)

(use srfi-1)

(define (mk-grid n m)
  (let loop ((j 0) (ls '()))
    (if (> j m)
      ls
      (loop (+ j 1)
            (append
              ls
              (list (let loopi ((i 0) (lsi '()))
                      (if (> i n)
                        lsi
                        (loopi (+ i 1)
                               (cons (list (- (random-real) 0.5)
                                           (- (random-real) 0.5))
                                     lsi))))))))))

(define (dot2 v w) (+ (* (car v) (car w)) (* (cadr v) (cadr w))))

(define old-floor floor)

(define old-ceiling ceiling)

(define (floor x) (inexact->exact (old-floor x)))

(define (ceiling x) (inexact->exact (old-floor (+ x 1.0))))


(define (ingrid x y grid)
  (assert (> (length grid) (ceiling x)))
  (assert (> (length (car grid)) (ceiling y)))
  (assert (not (= (floor x) (ceiling x))))
  (assert (not (= (floor y) (ceiling y))))
  (define (ease d) (- (* 3 (* d d)) (* 2 (* d d))))
  (define (gridnode a b) (list-ref (list-ref grid a) b))
  (let* ((x0 (- x (floor x)))
         (x1 (- (+ 1.0 (ceiling x)) x))
         (y0 (- y (floor y)))
         (y1 (- (+ 1.0 (ceiling y)) y)))
    (let* ((d0 (dot2 (list x0 y0) (gridnode (floor x) (floor y))))
           (d1 (dot2 (list x0 y1) (gridnode (floor x) (ceiling y))))
           (d2 (dot2 (list x1 y0) (gridnode (ceiling x) (floor y))))
           (d3 (dot2 (list x1 y1) (gridnode (ceiling x) (ceiling y)))))
      (+ (/ (+ d0 (ease (- d0 d1))) 2.0) (/ (+ d2 (ease (- d2 d3))) 2.0)))))

(define (to-matlab gr num) 
  (define (allbutlast ls)
	(take ls (- (length ls) 1)))
  (define (linspace start stop n)
    (let ((stepsize (/ (- stop start) (exact->inexact n))))
      (let loop ((i start) (ls '()))
        (if (>= i stop) ls (loop (+ i stepsize) (append ls (list i)))))))
  (let ((maxx (exact->inexact (- (length (car gr)) 1)))
        (maxy (exact->inexact (- (length gr) 1))))
    (let ((xs (allbutlast (linspace 0.0 maxx num))) (ys (allbutlast (linspace 0.0 maxy num))))
      (define (print-xmat)
        (begin
          (print "X = [")
          (map (lambda (y)
                 (begin
                   (print "[")
                   (print (fold (lambda (x acc)
                                  (string-append acc " " (number->string x)))
                                ""
                                xs))
                   (print "],")))
               ys)
          (print "];")))
      (define (print-ymat)
        (begin
          (print "Y = [")
          (map (lambda (y)
                 (begin
                   (print "[")
                   (print (fold (lambda (x acc)
                                  (string-append acc " " (number->string y)))
                                ""
                                xs))
                   (print "],")))
               ys)
          (print "];")))
      (define (print-zmat)
        (begin
          (print "Z = [")
          (map (lambda (y)
                 (begin
                   (print "[")
                   (print (fold (lambda (x acc)
                                  (begin
                                    (string-append
                                      acc
                                      " "
                                      (number->string (ingrid x y gr)))))
                                ""
                                xs))
                   (print "],")))
               ys)
          (print "];")))
      (print-xmat)
      (print-ymat)
      (print-zmat)
      (print "surfc(X,Y,Z); sleep(200);"))))



(define dim 20.0)

(define gr (mk-grid dim dim))

(to-matlab gr 30)