(defun grovenew (x y sd &optional (mu 1.5181932)
		 (tau 0.0036737733))
  (let* ((u (/ tau sd))
	 (v (/ (* (- x y)(- x y) ) (* -4 sd sd) ))
	 (w (/ (* (- y mu)(- y mu) ) (* 2 tau tau)))
	 (lr (* u  (exp ( + v w )))))
	 lr))

(defun grove-main (datarray filename)
  (let ((lm (make-array (list (array-dimension datarray 0)
                              (array-dimension datarray 0))
                        :initial-element nil))
        (tempval 0) (totalpairs 0) (totalmatches 0)
        (tempsd 1) (prev 1) (curr 0) (count 0))
    (with-open-file (out filename
                         :direction :output
                         :if-exists :supersede)
      (format out "### splot '~a' index 0" filename) 
      (dotimes (i (array-dimension lm 0))
        (setf curr (aref datarray i 2))
        (cond ((not (= prev curr))
               (setf count (+ count 1))
               (format out ", '~a' index ~a"
                       filename count)))
        (setf prev curr))
      (format out "~%####SD:1~%")
      (dotimes (i (array-dimension lm 0))
        (if (not (= tempsd (aref datarray i 2)))
            (format out "~%~%###SD:~a~%" (aref datarray i 2)))
        (setf tempsd (aref datarray i 2))
        (dotimes (j (array-dimension lm 0))
          (if (equalp (aref lm i j) nil)
              (cond ((= (aref datarray i 2) (aref datarray j 2))
                     (handler-case
			 (progn
			   (setf tempval
				 (grovenew (aref datarray i 1)
					  (aref datarray j 1)
					  (* (aref datarray i 2)
					     (expt 10 -5)))))
                       (floating-point-overflow ()
                           (setf tempval 10)))
		     (if (> tempval 10)
			 (setf tempval 10))
                     (setf (aref lm i j) tempval)
                     (setf (aref lm j i) tempval)
                     (if (not (= i j))
                         (setf totalpairs (+ totalpairs 1))) 
                     (cond ((and (>= tempval 10)
                                 (not (= i j)))
                            (setf totalmatches (+ totalmatches 1))))
                     (if (not (= i j))
                         (with-standard-io-syntax
                           (format out
                                   "~a ~a ~a~%"
                                   (aref datarray i 1)
                                   (aref datarray j 1)
                                   tempval))))))))
      (format t "Total Matches: ~a~%"  totalmatches)
      (format t "Total Pairs:   ~a~%"  totalpairs)
      'done)))

(defun grove-random (min max min.sd max.sd 
			   timespersd filename &optional
			   (mu 1.5181932) (tau 0.0036737733))
  (let ((x 0)
	(y 0)
	(sd 0)
	(tempval 0))
    (with-open-file (out filename
			 :direction :output
			 :if-exists :supersede)
      (format out "### splot '~a' index 0" filename)
      (dotimes (i (- max.sd min.sd))
	(format out ", '~a' index ~a" filename (+ i 1)))
      (dotimes (i (+ (- max.sd min.sd) 1))
	(format out "~%~%###SD: ~a~%" (+ i 1))
	(dotimes (j timespersd)
	  (setf x (+ (random (- max min)) min))
	  (setf y (+ (random (- max min)) min))
	  (handler-case
	      (progn
		(setf tempval 
		      (grovenew x y (* (+ i 1)  
				      (expt 10 -5)) mu tau)))
	    (floating-point-overflow ()
	      (setf tempval 10)))
	  (if (> tempval 10)
	      (setf tempval 10))
	  (format out "~a ~a ~a~%" x y tempval))))
    'done))