(defun seheult (x y sd &optional (mu 1.5181932)
                (tau 0.0036737733))
  (let* ((l (/ sd tau))
         (ll (* l l))
         (u (/ (- x y) (* sd (sqrt 2))))
         (z (* .5 (+ x y)))
         (v (/ (- z mu) (* tau (expt (+ 1 (* .5 ll)) .5))))
         (factor (* (/ (+ 1 ll) (* l (expt (+ 2 ll) .5)))
                    (exp (* (/ -1 (* 2 (+ 1 ll)))
                            (- (* u u) (* v v)))))))
    factor))

(defun random-seheult (&optional (samples 10000) 
		       (min.x 1.4) (max.x 1.6)
		       (min.y 1.4) (max.y 1.6)
		       (min.sd 1) (max.sd 25)
		       (min.mu 1.5181932) (max.mu 15.181932)
		       (min.tau 0.0036737733) (max.tau 0.036737733))
  (labels ((nr (min max)
	     (+ min (random (- max min)))))
    (let ((listie '((x y sd mu tau lr)))
	  (x 0)
	  (y 0)
	  (sd 0)
	  (mu 0)
	  (tau 0)
	  (tempval 0))
      (dotimes (i samples)
	(setf  x (nr min.x max.x)
	       y (nr min.y max.y)
	       sd (* 1.e-5 (nr min.sd max.sd))
	       mu (nr min.mu max.mu) 
	       tau (nr min.tau max.tau))
	(handler-case
	    (progn
	      (setf tempval 
		    (seheult x y sd mu tau)))
	  
	  (floating-point-overflow ()
	    (setf tempval 10)))
	(if (> tempval 10)
	    (setf tempval 10))
	(setf listie (append listie (list (list x y sd mu tau tempval)))))
      listie)))

(defun write-seheult (datalist &optional (filename "randomsdata.lisp"))
  (with-open-file (out filename
		       :direction :output
		       :if-exists :supersede)
    (format out "(defun randomseheult () ~%'~A)" datalist))
  'done)