(use random-bsd)
(use srfi-1)

(define (mk-grid n m)
  (define (rand) (* 1.5 (- (random-real) 0.5)))
  (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 (rand) (rand)) lsi))))))))))

(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)
  (define (ease d) (- (* 3 (* d d)) (* 2 (* d d d))))
  (define (gridnode b a) (list-ref (list-ref grid a) b))
  (define (dot2 v w) (+ (* (car v) (car w)) (* (cadr v) (cadr w))))
  (let* ((x0 (- x (floor x)))
         (x1 (- x (ceiling x)))
         (y0 (- y (floor y)))
         (y1 (- y (ceiling y))))
    (let* ((s (dot2 (list x0 y0) (gridnode (floor x) (floor y))))
           (t (dot2 (list x0 y1) (gridnode (floor x) (ceiling y))))
           (u (dot2 (list x1 y0) (gridnode (ceiling x) (floor y))))
           (v (dot2 (list x1 y1) (gridnode (ceiling x) (ceiling y))))
           (sx (ease x0))
           (sy (ease y0))
           (a (+ s (* sy (- t s))))
           (b (+ u (* sy (- v u)))))
      (+ a (* sx (- b a))))))


(define (to-matlab gr num)
  (define (allbutlast ls) (take ls (- (length ls) 1)))
  (define (ls->string ls)
    (fold (lambda (x acc) (string-append acc " " (number->string x))) "" ls))
  (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 (- (length (car gr)) 1))
         (maxy (- (length gr) 1))
         (xs (allbutlast (linspace 0.0 maxx num)))
         (ys (allbutlast (linspace 0.0 maxy num))))
    (define (print-xmat)
      (print "Xrow = [" (ls->string xs) "];")
      (print "X = repmat(Xrow," (length ys) ",1);"))
    (define (print-ymat)
      (print "Ycol = [" (ls->string ys) "];")
      (print "Y = repmat(Ycol', 1, " (length xs) ");"))
    (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);")))

(set! gr (mk-grid 10 10.0))
(to-matlab gr 40)