;;; -*- Mode: LISP; Syntax: Common-lisp; Base: 10 -*-
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Copyright (c) 92, 93, 94, 98 Fernando Lopez Lezcano.  All rights reserved.
;;; Use and copying of this software and preparation of derivative works
;;; based upon this software are permitted and may be copied as long as 
;;; no fees or compensation are charged for use, copying, or accessing
;;; this software and all copies of this software include this copyright
;;; notice.  Suggestions, comments and bug reports are welcome.  Please 
;;; address email to: nando@ccrma.stanford.edu
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Dynamic quadraphonic signal locator
;;;
;;; by Fernando Lopez Lezcano
;;;     CCRMA, Stanford University
;;;     nando@ccrma.stanford.edu
;;;

;;;;;;;;;;
;;; To do:
;;;
;;; * gain and rgain envelopes are not really necessary. Would be better to replace them
;;;   with an interpolated array representing gain at different distance which could be
;;;   accesed with distance of delay envelope.
;;; v compensate for non square distribution of speakers (done 5/4/93)
;;; * use phase information, in particular use three speakers at a time to locate sound. 
;;; v put gain normalization as optional parameter
;;; * add optional additional delay on back speakers (6/21/93)
;;; * change supersonic movement check to check only radial velocity (6/24/93)
;;;
;;;;;;;;;;;;
;;; Changes:
;;;
;;; - Changed ,in-sig in the run macro to prevent multiple evaluation with side effects
;;; - Setup gain normalization as optional with a controlling parameter
;;; - Compensation for non square loudspeaker distribution
;;; - fixed bug due to precision of clm envelopes
;;; - changed pan table law to constant power

(defvar speed-of-sound 344)		;in air in meters per second (under normal conditions)

(defvar sound-mag 0)

(defun set-sound-mag ()			;conversion factor for meters/sec to samples/meter
  (setf sound-mag (/ sampling-rate 
		     speed-of-sound)))

(defun in-samples (x)			;x = distance in meters, converts to samples
  (* x sound-mag))

(defun in-seconds (x)                   ;x = distance in meters, convert to seconds of delay
  (/ x speed-of-sound))

(defun distance (x y)			;pythagorean theorem
  (sqrt (+ (* x x) (* y y))))

(defun nearest-point (x0 y0 x1 y1 sx sy)
  ;; find the point of intersection with the line between (x0 y0) and (x1 y1) and a 
  ;; perpendicular to it dropped from (sx sy) this gives the first part of the calculation
  ;; of the virtual source (i.e. reflection) given any current source (sx sy) and any line
  ;; the line equation is y = (d-b)/(c-a) x + (bc-ad)/(c-a) 
  ;; the perpendicular to it that passes through (e f) is y = f + (c-a)/(d-b) e - (c-a)/(d-b) x
  ;; we can solve these two equations for x, then plug in that value to get y (to get the point
  ;; of intersection) x = (f + (c-a)/(d-b) e - (bc-ad)/(c-a)) / ((d-b)/(c-a) + (c-a)/(d-b))
  ;; all this algebra is just an explicit solution of the vector equations given in Borish's article.
  (if (= x0 x1)
      (values x0 sy)
    (if (= y0 y1)
	(values sx y0)
      (let* ((c-a (- x1 x0))
	     (d-b (- y1 y0))
	     (s1 (/ c-a d-b))
	     (s2 (/ d-b c-a))
	     (bc (* y0 x1))
	     (ad (* x0 y1))
	     (offset (/ (- bc ad) c-a))
	     (inX (/ (+ (* s1 sx) sy (- offset)) (+ s1 s2)))
	     (inY (+ (* inX s2) offset)))
	(values inX inY)))))

(defun make-gpolar (R freq 
		    &key (gain (* (- 1 R)(- 1 R))))
  (make-two-pole gain 
		 (- (* 2.0 R (cos (* freq frequency-mag)))) 
		 (* R R)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Define the loudspeaker control structure

;;; pan-total-length: length in table bins of the 360 degree circle
;;; pan-offset: offset in radians of the first speaker with respect to the y axis
;;; source-table-xy: pan interpolation table for space between x and y speakers
;;; source-lower-xy: lower bound of the xy table (in table bins)
;;; source-upper-xy: upper bound of the xy table (in table bins)

(defstruct pan
  total-length
  offset
  table-ab
  table-bc
  table-cd
  table-da
  lower-ab
  lower-bc
  lower-cd
  lower-da
  upper-ab
  upper-bc
  upper-cd
  upper-da)
  
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Create a new loudspeaker control structure

(defparameter dlocsig-pan nil)

;;; total-table-length: number of table bins the circle is split into
;;;                     (/ 360 total-table-length) gives resolution of 
;;;                     loudspeaker placement in degrees.
;;; aspect-ratio: desired ratio of loudspeaker rectangle (depth / width)
;;; angle-x: angle of speaker x (a/b/c/d) with respect to y axis

(defun make-pan-tables (&key
			(total-table-length 120) 
			(aspect-ratio nil)
			(angle-a nil)
			(angle-b nil)
			(angle-c nil)
			(angle-d nil))
  (flet ((norm(angle)
	   (if angle 
	       (mod angle (* 2 pi)) 
	     nil)))
    (if (not (or aspect-ratio angle-a angle-b angle-c angle-d))
	(error "~%Error: aspect ratio or angles must be specified!"))
    (let* (ang-ab ang-bc ang-cd ang-da offset)
      (setf angle-a (norm angle-a)
	    angle-b (norm angle-b)
	    angle-c (norm angle-c)
	    angle-d (norm angle-d))
      (if aspect-ratio
	  (let ((beta (atan 1 aspect-ratio)))
	    (if (or angle-a angle-b angle-c angle-d)
		(error "~%Error: both aspect ratio and angles specified!")
	      (setf offset beta
		    ang-ab (* 2 beta)
		    ang-bc (- pi beta beta)
		    ang-cd (* 2 beta)
		    ang-da (- pi beta beta))))
	(if (not (and angle-a angle-b angle-c angle-d))
	    (error "~%Error: unspecified speaker angle! a:~s b:~s c:~s d:~s"
		   angle-a angle-b angle-c angle-d)
	  (if (not (< angle-b angle-c angle-d angle-a ))
	      (error "~%Error: unordered speaker angle! b:~s c:~s d:~s a:~s"
		     angle-b angle-c angle-d angle-a)
	    (setf offset angle-a 
		  ang-ab (norm (- angle-b angle-a))
		  ang-bc (norm (- angle-c angle-b))
		  ang-cd (norm (- angle-d angle-c))
		  ang-da (norm (- angle-a angle-d))))))
      (let* ((len-ab (round (* (/ ang-ab (* 2 pi)) total-table-length)))
	     (len-bc (round (* (/ ang-bc (* 2 pi)) total-table-length)))
	     (len-cd (round (* (/ ang-cd (* 2 pi)) total-table-length)))
	     (len-da (round (* (/ ang-da (* 2 pi)) total-table-length))))
	(make-pan :offset offset
		  :lower-ab 0
		  :upper-ab len-ab
		  :table-ab (let ((tab (make-table (+ len-ab 1))))
			      (loop  for index from 0 to len-ab and angle from (- (/ pi 4)) by (/ pi 2 len-ab) do 
				(setf (aref tab index) (make-short-float (- (cos angle)(sin angle)))))
			      tab)
		  :lower-bc len-ab
		  :upper-bc (+ len-ab len-bc)
		  :table-bc (let ((tab (make-table (+ len-bc 1))))
			      (loop  for index from 0 to len-bc and angle from (- (/ pi 4)) by (/ pi 2 len-bc) do 
				(setf (aref tab index) (make-short-float (- (cos angle)(sin angle)))))
			      tab)
		  :lower-cd (+ len-ab len-bc)
		  :upper-cd (+ len-ab len-bc len-cd)
		  :table-cd (let ((tab (make-table (+ len-cd 1))))
			      (loop  for index from 0 to len-cd and angle from (- (/ pi 4)) by (/ pi 2 len-cd) do 
				(setf (aref tab index) (make-short-float (- (cos angle)(sin angle)))))
			      tab)
		  :lower-da (+ len-ab len-bc len-cd)
		  :upper-da (+ len-ab len-bc len-cd len-da)
		  :table-da (let ((tab (make-table (+ len-da 1))))
			      (loop  for index from 0 to len-da and angle from (- (/ pi 4)) by (/ pi 2 len-da) do 
				(setf (aref tab index) (make-short-float (- (cos angle)(sin angle)))))
			      tab)
		  :total-length (+ len-ab len-bc len-cd len-da))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Define the dlocsig structure

;;; dlocs-sources: number of sources that will be modeled (source 0 is the real source)
;;; dlocs-start: array of starting samples for each source (sample at which signal arrives to the listener position)
;;; dlocs-angle: array of angle envelopes 
;;; dlocs-gain: array of attenuation envelopes 
;;;     gain>1 --> virtual source inside speaker perimeter (only happens to real source!)
;;;     gain<=1 --> virtual source outside speaker perimeter (normal case)
;;; dlocs-rgain: array of envelopes containing the gain for the reverberator channel
;;; dlocs-delay: array of envelopes that describe the axial movement of the virtual source
;;; dlocs-path: delay line representing the path between sources and listener
;;; dlocs-pan-xy: pan interpolation table from speaker x to speaker y
;;; dlocs-upper-xy: max value of angle parameter for zone from speaker x to y

(def-clm-struct dlocs
  (start integer)
  (end integer)
  (angle envelope)
  (gain envelope)
  (rgain envelope)
  (delay envelope)
  (pan-ab array real)
  (pan-bc array real)
  (pan-cd array real)
  (pan-da array real)
  (upper-ab integer)
  (upper-bc integer)
  (upper-cd integer)
  (upper-da integer)
  (scaler real)
  (path zdly)

  rev)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Include virtual image handling

(def-clm-struct dlocs-virtual
  (sources integer)
  (start array integer)
  (end integer)
  (angle array envelope)
  (gain array envelope)
  (rgain array envelope)
  (delay array envelope)
  (filter array smpflt)
  (pan-ab array real)
  (pan-bc array real)
  (pan-cd array real)
  (pan-da array real)
  (upper-ab integer)
  (upper-bc integer)
  (upper-cd integer)
  (upper-da integer)
  (scaler real)
  (path zdly)

  rev
  only-to)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Define the dlocsig default parameters

(defparameter dlocsig-path '())
(defparameter dlocsig-scaler 1.0)
(defparameter dlocsig-power nil)
(defparameter dlocsig-reverb-power nil)
(defparameter dlocsig-initial-delay nil)

(defparameter dlocsig-x-listener 5)
(defparameter dlocsig-y-listener 5)
(defparameter dlocsig-W 10)
(defparameter dlocsig-L 20)
(defparameter dlocsig-max-time 0.1)

(defparameter dlocsig-reflect-coeff 1.0)
(defparameter dlocsig-front-reflect-coeff 0.9)
(defparameter dlocsig-left-reflect-coeff 0.9)
(defparameter dlocsig-right-reflect-coeff 0.9)
(defparameter dlocsig-back-reflect-coeff 0.9)

(defparameter dlocsig-filter-coeff 1.0)
(defparameter dlocsig-front-filter-coeff 0.9)
(defparameter dlocsig-left-filter-coeff 0.9)
(defparameter dlocsig-right-filter-coeff 0.9)
(defparameter dlocsig-back-filter-coeff 0.9)

(defparameter dlocsig-filter-freq 3000)
(defparameter dlocsig-filter-R 0.5)

(defparameter dlocsig-reverb-amount 0.04)

(defparameter dlocsig-normalize-gain nil)
(defparameter dlocsig-max-rooms 1)
(defparameter dlocsig-only-to 0)
(defparameter dlocsig-verbose t)

(defparameter dlocsig-speaker-aspect-ratio 1.5)

(defparameter dlocsig-debug '())
(defmacro IFDEBUG (lev x)
  #+mcl (declare (ignore x)) 
  #+(or excl gcl) `(if (if (listp dlocsig-debug)
			   (member ,lev dlocsig-debug)
			 (<= ,lev dlocsig-debug))
		       ,x))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Create a new dlocsig structure

(defun make-dlocsig (&key 
		    (start-time 0)
		    (duration 1)
		    (path dlocsig-path)
		    (scaler dlocsig-scaler)
		    (power dlocsig-power)
		    (rev-power dlocsig-reverb-power)
		    (reverb-amount dlocsig-reverb-amount)
		    (initial-delay dlocsig-initial-delay)
		    (normalize-gain dlocsig-normalize-gain)
		    (aspect-ratio dlocsig-speaker-aspect-ratio)
		    (verbose dlocsig-verbose))
  (let* ((start 0)
	 (dist '())
	 (angle '())
	 (gain '())
	 (rgain '())
	 (delay '())
	 (sound-beg 0)
	 (sound-dur 0)
	 (prev-angle 0.0)
	 (prev-time -1.0)
	 (max-dist 0.0)
	 (max-delay 0.0)
	 (min-dist 999999)
	 (unity-dist 1.0)
	 (xpoints '())
	 (ypoints '())
	 (tpoints '())
	 (xrpoints '())
	 (yrpoints '())
	 (trpoints '())
	 (run-beg start-time)
	 (run-end 0.0))

    ;; Set value of individual component lists from the path
    (setf tpoints (get-t-coords path)
	  xpoints (get-x-coords path)
	  ypoints (get-y-coords path))

    ;; Find the minimum distance of the real sound source
    (setf min-dist (loop 
		       for x in xpoints
		       for y in ypoints
		       minimize (distance x y)))
    (if normalize-gain (setf unity-dist min-dist))

    ;; Create reverse order path lists for pushing in main loop
    (setf xrpoints (reverse xpoints)
	  yrpoints (reverse ypoints)
	  trpoints (reverse tpoints))

    ;; Check the envelope for supersonic movement of the source,
    ;; we only check radial movement which is the one that affects
    ;; resampling in the delay line
    (loop 
	for ptime in tpoints
	for px in xpoints
	for py in ypoints
	for time in (cdr tpoints)
	for x in (cdr xpoints)
	for y in (cdr ypoints)
	for speed = 0 
	for ticket = nil 
	for span = (- (car (last tpoints))(car tpoints)) do
	  (if (/= 0 (distance (- x px)(- y py)))
	      (if (setf ticket 
		    (> (setf speed 
			 (/ (abs (- y py))
			    (* (/ (- time ptime) span) duration)))
		       speed-of-sound))
		  (format t "supersonic radial movement from [~f,~f,~f] to [~f,~f,~f], speed=~f~%" 
			  ptime px py time x y speed)))
	finally (if ticket (error "ticket due to overspeed!~%")))

    ;; Create the loudspeaker pan tables if necessary
    (if (not dlocsig-pan)
	(setf dlocsig-pan (make-pan-tables :aspect-ratio aspect-ratio)))
    
    ;; Loop through all virtual rooms for one breakpoint in the trajectory of the source.
    ;; Looping is done for only one quadrant. In the inner loop the function add-break-to-room
    ;; is called four times for the four rooms defined in each quadrant

    (flet ((walk-all-rooms (xs ys time)
	     (let ((room 0))
	       (IFDEBUG 1 (format t "source[x=~f y=~f t=~f]~%" xs ys time))
	       
	       (flet ((add-break-to-room (x y att q)
			(if (/= 0 att)
			    (let* ((d (if initial-delay 
					  (distance x y)
					(- (distance x y) min-dist)))
				   (previous-angle prev-angle)
				   (previous-time prev-time)
				   (arrival-time (+ time (/ d speed-of-sound)))
				   (current-angle (coerce (* (pan-upper-da dlocsig-pan)
							     (/ (mod (+ (atan x y)
									(pan-offset dlocsig-pan))
								     (* pi 2))
								(* pi 2)))
							  'short-float)))
			      (IFDEBUG 1 (format t "x=~f, y=~f, d=~f, a=~f~%" x y d att))
			      (push d dist)
			      (push (in-samples d) delay)
			      (push time delay)
			      (setf max-dist (max d max-dist))
			      (if (/= previous-time -1.0)
				  (if (<= 0.0 current-angle (pan-upper-ab dlocsig-pan)
					  (pan-lower-da dlocsig-pan) previous-angle (pan-upper-da dlocsig-pan))
				      (let ((new-time (/ (+ (* current-angle previous-time)
							    (* (- (pan-upper-da dlocsig-pan) previous-angle) arrival-time))
							 (+ current-angle (- (pan-upper-da dlocsig-pan) previous-angle)))))
					(push (pan-upper-da dlocsig-pan) angle)
					(push new-time angle)
					(push 0.0 angle)
					(push new-time angle))
				    (if (<= 0.0 previous-angle (pan-upper-ab dlocsig-pan)
					    (pan-lower-da dlocsig-pan) current-angle (pan-upper-da dlocsig-pan))
					(let ((new-time (/ (+ (* (- (pan-upper-da dlocsig-pan) current-angle) previous-time)
							      (* previous-angle arrival-time))
							   (+ (- (pan-upper-da dlocsig-pan) current-angle) previous-angle))))
					  (push 0.0 angle)
					  (push new-time angle)
					  (push (pan-upper-da dlocsig-pan) angle)
					  (push new-time angle)))))
			      (push current-angle angle)
			      (push arrival-time angle)
			      (let ((base (/ unity-dist (if initial-delay d (+ d min-dist)))))
				(push (* att (if power 
						 (expt base power)
					       base))
				      gain)
				(push arrival-time gain)
				(push (* att (if rev-power 
						 (expt base rev-power)
					       (sqrt base)))
				      rgain)
				(push arrival-time rgain))
			      (setf prev-angle current-angle
				    prev-time arrival-time)
			      (incf room)))))
		       
		       ;; x>0 y>0
		       (add-break-to-room xs ys scaler 1))
	       ;; return number of rooms processed
	       room)))
      
      ;; Loop for each point in the position envelope <pairs are reversed so push is all right>
      ;; this-x, this-y, this-t = position of the source
      ;; next-x, next-y, next-t = next position of the source
      ;; new-x, new-y = position of intermediate breakpoint if segment changes direction
      ;; room = sequential room number

      (loop 
	  for this-x in xrpoints
	  for this-y in yrpoints
	  for this-t in trpoints 
	  for next-x in (cdr xrpoints)
	  for next-y in (cdr yrpoints)
	  for next-t in (cdr trpoints) do
	    (walk-all-rooms this-x this-y this-t)
	    (multiple-value-bind (new-x new-y)(nearest-point this-x this-y next-x next-y 0 0)
	      (IFDEBUG 8 (format t "this[~f,~f] new[~f,~f] next[~f,~f]~%" this-x this-y new-x new-y next-x next-y))
	      (if (and (if (< this-x next-x)
			   (<= this-x new-x next-x)
			 (<= next-x new-x this-x))
		       (if (< this-y next-y)
			   (<= this-y new-y next-y)
			 (<= next-y new-y this-y)))
		  (walk-all-rooms new-x new-y 
				  (+ next-t (* (- this-t next-t)
					       (/ (distance (- next-x new-x)(- next-y new-y))
						  (distance (- next-x this-x)(- next-y this-y))))))))
	  finally 
	    (walk-all-rooms next-x next-y next-t)))
    
    ;; Calculate starting times and samples
    (setf max-delay (in-samples max-dist)
	  sound-beg (+ start-time (in-seconds (car dist)))
	  sound-dur (+ duration (in-seconds (- (car (last dist))(car dist))))
	  start (in-samples (* sound-beg speed-of-sound))
	  run-end (+ sound-beg sound-dur))
    
    ;; Create the runtime structure with all the information
    (values (make-dlocs :start start
			:end (+ (floor (* start-time sampling-rate))
				(floor (* duration sampling-rate)))
			:angle (make-env :envelope angle
					 :start-time sound-beg
					 :duration sound-dur)
			:gain (make-env :envelope gain
					:start-time sound-beg
					:duration sound-dur)
			:rgain (make-env :envelope rgain
					 :start-time sound-beg
					 :duration sound-dur)
			:delay (make-env :envelope delay
					 :start-time sound-beg
					 :duration sound-dur)
			:path (make-zdelay 0 :true-length (max 1 (+ (ceiling max-delay) 1)))
			:pan-ab (pan-table-ab dlocsig-pan)
			:pan-bc (pan-table-bc dlocsig-pan)
			:pan-cd (pan-table-cd dlocsig-pan)
			:pan-da (pan-table-da dlocsig-pan)
			:upper-ab (pan-upper-ab dlocsig-pan)
			:upper-bc (pan-upper-bc dlocsig-pan)
			:upper-cd (pan-upper-cd dlocsig-pan)
			:upper-da (pan-upper-da dlocsig-pan)
			:scaler scaler
			:rev reverb-amount)
	    (floor (* run-beg sampling-rate))
	    (floor (* run-end sampling-rate)))))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Run macro to localize samples
;;;

(defmacro dlocsig (l i in-sig)
  
  ;; Generate names for internal variables so there is no collision
  (let ((input (gensym))
	(a-sample (gensym))
	(b-sample (gensym))
	(c-sample (gensym))
	(d-sample (gensym))
	(r-sample (gensym))
	(angle (gensym))
	(delay (gensym))
	(sample (gensym))
	(scaled-sample (gensym))) 
    
    ;; And now, let's locate the signal
    `(let ((,input ,in-sig)
	   (,a-sample 0.0)
	   (,b-sample 0.0)
	   (,c-sample 0.0)
	   (,d-sample 0.0)
	   (,r-sample 0.0))
       (if (< ,i (dlocs-start ,l))
	   ;; Feed the delay while waiting for the signal to appear at the output
	   (zdelay (dlocs-path ,l) (if (> ,i (dlocs-end ,l)) 0.0 ,input) 0.0)
	 ;; Signal coming out of the delay line... process it
	 (let* ((,angle (env (dlocs-angle ,l)))
		(,delay (env (dlocs-delay ,l)))
		(,sample (zdelay (dlocs-path ,l) (if (> ,i (dlocs-end ,l)) 0.0 ,input) ,delay))
		(,scaled-sample (* ,sample (env (dlocs-gain ,l)))))
	   (if (or (stereo)(quad))
	       (cond ((< ,angle 0.0)
		      (setf ,d-sample (* ,scaled-sample (array-interp (dlocs-pan-da ,l)(- (+ (dlocs-upper-da ,l) ,angle)
											  (dlocs-upper-cd ,l)))))
		      (setf ,a-sample (* ,scaled-sample (array-interp (dlocs-pan-da ,l)(- ,angle)))))
		     ((<= ,angle (dlocs-upper-ab ,l))
		      (setf ,a-sample (* ,scaled-sample (array-interp (dlocs-pan-ab ,l) ,angle)))
		      (setf ,b-sample (* ,scaled-sample (array-interp (dlocs-pan-ab ,l)(- (dlocs-upper-ab ,l) ,angle)))))
		     ((<= ,angle (dlocs-upper-bc ,l))
		      (setf ,b-sample (* ,scaled-sample (array-interp (dlocs-pan-bc ,l)(- ,angle (dlocs-upper-ab ,l)))))
		      (setf ,c-sample (* ,scaled-sample (array-interp (dlocs-pan-bc ,l)(- (dlocs-upper-bc ,l) ,angle)))))
		     ((<= ,angle (dlocs-upper-cd ,l))
		      (setf ,c-sample (* ,scaled-sample (array-interp (dlocs-pan-cd ,l)(- ,angle (dlocs-upper-bc ,l)))))
		      (setf ,d-sample (* ,scaled-sample (array-interp (dlocs-pan-cd ,l)(- (dlocs-upper-cd ,l) ,angle)))))
		     ((<= ,angle (dlocs-upper-da ,l))
		      (setf ,d-sample (* ,scaled-sample (array-interp (dlocs-pan-da ,l)(- ,angle (dlocs-upper-cd ,l)))))
		      (setf ,a-sample (* ,scaled-sample (array-interp (dlocs-pan-da ,l)(- (dlocs-upper-da ,l) ,angle)))))
		     ((> ,angle (dlocs-upper-da ,l))
		      (setf ,a-sample (* ,scaled-sample (array-interp (dlocs-pan-ab ,l) (- ,angle (dlocs-upper-da ,l)))))
		      (setf ,b-sample (* ,scaled-sample (array-interp (dlocs-pan-ab ,l)(- (dlocs-upper-ab ,l) 
											  (- ,angle (dlocs-upper-da ,l)))))))
		     (t (clm-print "error: angle ~f out of bounds~%" ,angle)))
	     (setf ,a-sample ,scaled-sample))
	   ;; compute coefficient for the reverberator channel 
	   (setf ,r-sample (* ,sample (dlocs-rev ,l)(env (dlocs-rgain ,l))))))
       ;; send sample to output streams, fold back channels if stereo
       (outa ,i (if (stereo) 
		    (/ (+ ,a-sample ,d-sample) 2)
		  ,a-sample))
       (if (or (stereo)(quad))
	   (outb ,i (if (stereo)
			(/ (+ ,b-sample ,c-sample) 2)
		      ,b-sample)))
       (if (quad)
	   (progn (outc ,i ,c-sample)
		  (outd ,i ,d-sample)))
       (if *reverb* (outa ,i ,r-sample *reverb*)))))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Create a new dlocsig structure with virtual image handling
;;;

(defun make-dlocsig-virtual (&key 
			     (start-time 0)
			     (duration 1)
			     (path dlocsig-path)
			     (x-listener dlocsig-x-listener)
			     (y-listener dlocsig-y-listener)
			     (W dlocsig-W)
			     (L dlocsig-L)
			     (max-time dlocsig-max-time)
			     (max-rooms dlocsig-max-rooms)
			     (only-to dlocsig-only-to)
			     (scaler dlocsig-scaler)
			     (power dlocsig-power)
			     (rev-power dlocsig-reverb-power)
			     (reflect-coeff dlocsig-reflect-coeff)
			     (front-reflect-coeff (* reflect-coeff dlocsig-front-reflect-coeff))
			     (left-reflect-coeff (* reflect-coeff dlocsig-left-reflect-coeff))
			     (right-reflect-coeff (* reflect-coeff dlocsig-right-reflect-coeff))
			     (back-reflect-coeff (* reflect-coeff dlocsig-back-reflect-coeff))
			     (filter-coeff dlocsig-filter-coeff)
			     (front-filter-coeff (* filter-coeff dlocsig-front-filter-coeff))
			     (left-filter-coeff (* filter-coeff dlocsig-left-filter-coeff))
			     (right-filter-coeff (* filter-coeff dlocsig-right-filter-coeff))
			     (back-filter-coeff (* filter-coeff dlocsig-back-filter-coeff))
			     (filter-freq dlocsig-filter-freq)
			     (filter-R dlocsig-filter-R)
			     (reverb-amount dlocsig-reverb-amount)
			     (initial-delay dlocsig-initial-delay)
			     (normalize-gain dlocsig-normalize-gain)
			     (aspect-ratio dlocsig-speaker-aspect-ratio)
			     (verbose dlocsig-verbose))

  (let* ((radius (* max-time speed-of-sound))
	 (rooms (loop 
		    repeat (ceiling radius L) 
		    for i from 0 for n = 0 do
		      (loop 
			  repeat (ceiling (sqrt (/ (+ (* radius radius)(* i i W W))(* L L))))
			  for j from 0 do
			    (incf n (if (and (= 0 i)(= 0 j)) 1
				      (if (or (= 0 i)(= 0 j)) 2 4))))
		    finally (return n)))
	 (start (make-array rooms :element-type 'integer :adjustable t))
	 (dist (make-array rooms :adjustable t))
	 (angle (make-array rooms :adjustable t))
	 (gain (make-array rooms :adjustable t))
	 (rgain (make-array rooms :adjustable t))
	 (delay (make-array rooms :adjustable t))
	 (filter (make-array rooms :adjustable t))
	 (sound-beg (make-array rooms))
	 (sound-dur (make-array rooms))
	 (prev-angle (make-array rooms :element-type 'float :initial-element 0.0))
	 (prev-time (make-array rooms :element-type 'float :initial-element -1.0))
	 (max-dist (make-array rooms :element-type 'float :initial-element 0.0))
	 (max-delay 0.0)
	 (min-dist (+ radius (distance W L)))
	 (unity-dist 1.0)
	 (xpoints '())
	 (ypoints '())
	 (tpoints '())
	 (xrpoints '())
	 (yrpoints '())
	 (trpoints '())
	 (run-beg start-time)
	 (run-end 0.0))

    ;; Set value of individual component lists from the path
    (setf tpoints (get-t-coords path)
	  xpoints (get-x-coords path)
	  ypoints (get-y-coords path))

    ;; Find the minimum distance of the real sound source
    (setf min-dist (loop 
		       for x in xpoints
		       for y in ypoints
		       minimize (distance x y)))
    (if normalize-gain (setf unity-dist min-dist))

    ;; Create reverse order path lists for pushing in main loop
    (setf xrpoints (reverse xpoints)
	  yrpoints (reverse ypoints)
	  trpoints (reverse tpoints))

    ;; Check the envelope for supersonic movement of the source
    ;; we only check radial movement which is the one that affects
    ;; resampling in the delay line
    (loop 
	for ptime in tpoints
	for px in xpoints
	for py in ypoints
	for time in (cdr tpoints)
	for x in (cdr xpoints)
	for y in (cdr ypoints)
	for speed = 0 
	for ticket = nil 
	for span = (- (car (last tpoints))(car tpoints)) do
	  (if (/= 0 (distance (- x px)(- y py)))
	      (if (setf ticket 
		    (> (setf speed 
			 (/ (abs (- y py))
			    (* (/ (- time ptime) span) duration)))
		       speed-of-sound))
		  (format t "supersonic movement from [~f,~f,~f] to [~f,~f,~f], speed=~f~%" 
			  ptime px py time x y speed)))
	finally (if ticket (error "ticket due to overspeed!~%")))

    ;; Check the envelope for out of the room movements and redimension room if necessary
    (if (> max-rooms 1)
	(loop 
	    for time in tpoints
	    for x in xpoints
	    for y in ypoints
	    for low-x = (- x-listener)
	    for high-x = (- W x-listener)
	    for low-y = (- y-listener)
	    for high-y = (- L y-listener) 
	    for new-W = W 
	    for new-L = L 
	    for changed = nil do
	      (if (not (and (< low-x x high-x)
			    (< low-y y high-y)))
		  (setf low-x (min low-x x)
			low-y (min low-y y)
			high-x (max high-x x)
			high-y (max high-y y)
			new-W (- high-x low-x)
			new-L (- high-y low-y)
			x-listener (- low-x)
			y-listener (- low-y)
			changed t))
	    finally (progn
		      (if (and changed verbose)
			  (warn "Source outside room: size changed from [~s:~s] to [~s:~s]" W L new-W new-L))
		      (setf W new-W
			    dlocsig-W new-W
			    L new-L
			    dlocsig-L new-L))))

    ;; Create the loudspeaker pan tables if necessary
    (if (not dlocsig-pan)
	(setf dlocsig-pan (make-pan-tables :aspect-ratio aspect-ratio)))
    
    ;; Loop through all virtual rooms for one breakpoint in the trajectory of the source.
    ;; Looping is done for only one quadrant. In the inner loop the function add-break-to-room
    ;; is called four times for the four rooms defined in each quadrant

    (flet ((walk-all-rooms (xs ys time)
	     (let ((room 0))
	       (IFDEBUG 1 (format t "source[x=~f y=~f t=~f]~%" xs ys time))
	       ;; Loop for all rooms in the x <width> direction
	       ;; xe?/xo? x coordinate of virtual source for even/odd rooms
	       ;; x?p/x?m x coordinate for rooms to the right/left of the real room
	       ;; rx? reflection coefficient of walls
	       ;; fx? filtering coefficient of walls
	       (loop 
		   repeat (ceiling radius L)
		   for i from 0 
		   for xep = xs 
		   for xop = (- (* 2 (- W x-listener)) xs)
		   for xem = (- xs (* 2 W))
		   for xom = (- (+ (* 2 x-listener) xs))
		   for rxp = 1
		   for rxm = 1 do
		     ;; Loop for all rooms in the y <length> direction for each x <width> value
		     ;; ye?/yo? coordinate of virtual source for even/odd rooms
		     ;; y?p/y?m coordinate for rooms on from/back of the real room
		     (loop 
			 repeat (ceiling (sqrt (/ (+ (* radius radius)(* i i W W))(* L L))))
			 for j from 0 
			 for yep = ys 
			 for yop = (- (* 2 (- L y-listener)) ys)
			 for yem = (- ys (* 2 L))
			 for yom = (- (+ (* 2 y-listener) ys))
			 for ryp = 1
			 for rym = 1 do
			   (IFDEBUG 3 (format t "  x: xep=~f xop=~f xem=~f xom=~f ~%" xep xop xem xom))
			   (IFDEBUG 3 (format t "  y: yep=~f yop=~f yem=~f yom=~f ~%" yep yop yem yom))
			   (let* (xp yp xm ym)
			     ;; decide which coordinates to use depending on room we are in <even/odd>
			     (if (evenp i)
				 (setf xp xep xm xem)
			       (setf xp xop xm xom))
			     (if (evenp j)
				 (setf yp yep ym yem)
			       (setf yp yop ym yom))

			     ;; Add a breakpoint to a room. It pushes distance, angle, gain and delay
			     ;; information into the list corresponding to the room. Except on the 
			     ;; delay list the time is corrected for the arrival delay to the listener. 

			     (flet ((add-break-to-room (x y att q)
				      (if (and (/= 0 att)(< room max-rooms))
					  (let* ((d (if initial-delay 
							(distance x y)
						      (- (distance x y) min-dist)))
						 (previous-angle (aref prev-angle room))
						 (previous-time (aref prev-time room))
						 (arrival-time (+ time (/ d speed-of-sound)))
						 (current-angle (coerce (* (pan-upper-da dlocsig-pan)
									   (/ (mod (+ (atan x y)
										      (pan-offset dlocsig-pan))
										   (* pi 2))
									      (* pi 2)))
									'short-float)))
					    (IFDEBUG 1 (format t "room[~d.~d]<~f:~f>, x=~f, y=~f, d=~f, a=~f~%" q room i j x y d att))
					    (push d (aref dist room))
					    (push (in-samples d)(aref delay room))
					    (push time (aref delay room))
					    (setf (aref max-dist room)(max d (aref max-dist room)))
					    (if (/= previous-time -1.0)
						(if (<= 0.0 current-angle (pan-upper-ab dlocsig-pan)
							(pan-lower-da dlocsig-pan) previous-angle (pan-upper-da dlocsig-pan))
						    (let ((new-time (/ (+ (* current-angle previous-time)
									  (* (- (pan-upper-da dlocsig-pan) previous-angle) arrival-time))
								       (+ current-angle (- (pan-upper-da dlocsig-pan) previous-angle)))))
						      (push (pan-upper-da dlocsig-pan) (aref angle room))
						      (push new-time (aref angle room))
						      (push 0.0 (aref angle room))
						      (push new-time (aref angle room)))
						  (if (<= 0.0 previous-angle (pan-upper-ab dlocsig-pan)
							  (pan-lower-da dlocsig-pan) current-angle (pan-upper-da dlocsig-pan))
						      (let ((new-time (/ (+ (* (- (pan-upper-da dlocsig-pan) current-angle) previous-time)
									    (* previous-angle arrival-time))
									 (+ (- (pan-upper-da dlocsig-pan) current-angle) previous-angle))))
							(push 0.0 (aref angle room))
							(push new-time (aref angle room))
							(push (pan-upper-da dlocsig-pan) (aref angle room))
							(push new-time (aref angle room))))))
					    (push current-angle (aref angle room))
					    (push arrival-time (aref angle room))
					    (let ((base (/ unity-dist (if initial-delay d (+ d min-dist)))))
					      (push (* att (if power 
							       (expt base power)
							     base))
						    (aref gain room))
					      (push arrival-time (aref gain room))
					      (push (* att (if rev-power 
							       (expt base rev-power)
							     (sqrt base)))
						    (aref rgain room))
					      (push arrival-time (aref rgain room)))
					    (setf (aref prev-angle room) current-angle
						  (aref prev-time room) arrival-time)
					    (incf room)))))
			       
			       (IFDEBUG 2 (format t "  xp=~f xm=~f yp=~f ym=~f~%" xp xm yp ym))
			       ;; x>0 y>0
			       (add-break-to-room xp yp (* rxp ryp scaler) 1)
			       ;; x>0 y<0
			       (if (/= 0 j)
				   (add-break-to-room xp ym (* rxp rym scaler) 2))
			       ;; x<0 y>0
			       (if (/= 0 i)
				   (add-break-to-room xm yp (* rxm ryp scaler) 3))
			       ;; x<0 y<0
			       (if (and (/= 0 i)(/= 0 j))
				   (add-break-to-room xm ym (* rxm rym scaler) 4))))
			   
			   ;; we are going to update distances and reflection coeffs for the next loop
			   ;; in the <y> direction...
			   (if (evenp j)
			       (progn (incf yep (* 2 L))
				      (decf yem (* 2 L))
				      ;; next room will be odd!
				      (setf ryp (* ryp front-reflect-coeff)
					    rym (* rym back-reflect-coeff)))
			     (progn (incf yop (* 2 L))
				    (decf yom (* 2 L))
				    ;; next room will be even!
				    (setf ryp (* ryp back-reflect-coeff)
					  rym (* rym front-reflect-coeff)))))
		     ;; in the <x> direction...
		     (if (evenp i)
			 (progn (incf xep (* 2 W))
				(decf xem (* 2 W))
				;; next room will be odd!
				(setf rxp (* rxp right-reflect-coeff)
				      rxm (* rxm left-reflect-coeff)))
		       (progn (incf xop (* 2 W))
			      (decf xom (* 2 W))
			      ;; next room will be even!
			      (setf rxp (* rxp left-reflect-coeff)
				    rxm (* rxm right-reflect-coeff)))))
	       ;; return number of rooms processed
	       room)))
      
      (IFDEBUG 1 (format t "radius=~f~%estimated rooms=~f~%" radius rooms))

      ;; Loop for each point in the position envelope <pairs are reversed so push is all right>
      ;; this-x, this-y, this-t = position of the source
      ;; next-x, next-y, next-t = next position of the source
      ;; new-x, new-y = position of intermediate breakpoint if segment changes direction
      ;; room = sequential room number

      (loop 
	  for this-x in xrpoints
	  for this-y in yrpoints
	  for this-t in trpoints 
	  for next-x in (cdr xrpoints)
	  for next-y in (cdr yrpoints)
	  for next-t in (cdr trpoints) do
	    (walk-all-rooms this-x this-y this-t)
	    (multiple-value-bind (new-x new-y)(nearest-point this-x this-y next-x next-y 0 0)
	      (IFDEBUG 8 (format t "this[~f,~f] new[~f,~f] next[~f,~f]~%" this-x this-y new-x new-y next-x next-y))
	      (if (and (if (< this-x next-x)
			   (<= this-x new-x next-x)
			 (<= next-x new-x this-x))
		       (if (< this-y next-y)
			   (<= this-y new-y next-y)
			 (<= next-y new-y this-y)))
		  (walk-all-rooms new-x new-y 
				  (+ next-t (* (- this-t next-t)
					       (/ (distance (- next-x new-x)(- next-y new-y))
						  (distance (- next-x this-x)(- next-y this-y))))))))
	  finally 
	    (setf rooms (walk-all-rooms next-x next-y next-t))))
    
    (adjust-array start rooms)
    (adjust-array dist rooms)
    (adjust-array gain rooms)
    (adjust-array rgain rooms)
    (adjust-array angle rooms)
    (adjust-array delay rooms)
    (adjust-array filter rooms)

    (IFDEBUG 20 (format t "~%angle=~s" (aref angle 0)))
    (IFDEBUG 20 (format t "~%gain=~s" (aref gain 0)))

    ;; Calculate starting times and samples for all the virtual sources
    (loop 
	for room from 0 repeat rooms do
	  (setf max-delay (max max-delay (in-samples (aref max-dist room)))
		(aref sound-beg room)(+ start-time (in-seconds (car (aref dist room))))
		(aref sound-dur room)(+ duration (in-seconds (- (car (last (aref dist room)))
								(car (aref dist room)))))
		(aref start room)(in-samples (* (aref sound-beg room) speed-of-sound)))
	  (IFDEBUG 12 (format t "start[~s]=~s; ~%dist=~s~%delay=~s~%" room (aref start room)(aref dist room)(aref delay room)))
	maximize (+ (aref sound-beg room)(aref sound-dur room)) into sound-end
	finally (setf run-end sound-end))

    (IFDEBUG 10 (format t "max delay = ~f~%" max-delay))
    (IFDEBUG 11 (format t "rooms:~f~%begin:~f~%end:~f~%" rooms run-beg run-end))

    ;; Create the runtime structure with all the information
    (values (make-dlocs :sources rooms
			:start start
			:end (+ (floor (* start-time sampling-rate))
				(floor (* duration sampling-rate)))
			:angle (loop for room from 0 repeat rooms do
				(setf (aref angle room)
				  (make-env :envelope (aref angle room)
					    :start-time (aref sound-beg room)
					    :duration (aref sound-dur room)))
				 finally (return angle))
			:gain (loop for room from 0 repeat rooms do
				(setf (aref gain room)
				  (make-env :envelope (aref gain room)
					    :start-time (aref sound-beg room)
					    :duration (aref sound-dur room)))
				finally (return gain))
			:rgain (loop for room from 0 repeat rooms do
				 (setf (aref rgain room)
				   (make-env :envelope (aref rgain room)
					     :start-time (aref sound-beg room)
					     :duration (aref sound-dur room)))
				finally (return rgain))
			:delay (loop for room from 0 repeat rooms do
				 (setf (aref delay room)
				   (make-env :envelope (aref delay room)
					     :start-time (aref sound-beg room)
					     :duration (aref sound-dur room)))
				finally (return delay))
			:filter (loop for room from 0 repeat rooms do
				  (setf (aref filter room)
				    (if (/= room 0)
					(make-gpolar filter-R filter-freq)
				      (make-gpolar filter-R (/ sampling-rate 2.2))))
				 finally (return filter))
			:path (make-zdelay 0 :true-length (max 1 (+ (ceiling max-delay) 1)))
			:pan-ab (pan-table-ab dlocsig-pan)
			:pan-bc (pan-table-bc dlocsig-pan)
			:pan-cd (pan-table-cd dlocsig-pan)
			:pan-da (pan-table-da dlocsig-pan)
			:upper-ab (pan-upper-ab dlocsig-pan)
			:upper-bc (pan-upper-bc dlocsig-pan)
			:upper-cd (pan-upper-cd dlocsig-pan)
			:upper-da (pan-upper-da dlocsig-pan)
			:rev reverb-amount
			:scaler scaler
			:only-to only-to)
	    (floor (* run-beg sampling-rate))
	    (floor (* run-end sampling-rate)))))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Run macro to localize samples with reflections on the walls

(defmacro dlocsig-virtual (l i in-sig)
  
  ;; Generate names for internal variables so there is no colission
  (let ((last-source (gensym))
	(input (gensym))
	(a-sample (gensym))
	(b-sample (gensym))
	(c-sample (gensym))
	(d-sample (gensym))
	(r-sample (gensym))
	(angle (gensym))
	(gain (gensym))
	(delay (gensym))
	(sample (gensym))
	(scaled-sample (gensym))) 
    
    ;; And now, let's locate the signal
    `(let ((,last-source (- (dlocs-sources ,l) 1))
	   (,input ,in-sig)
	   (,a-sample 0.0)
	   (,b-sample 0.0)
	   (,c-sample 0.0)
	   (,d-sample 0.0)
	   (,r-sample 0.0))
       (loop for source from 0 to ,last-source do
	 (if (< ,i (aref (dlocs-start ,l) source))
	     ;; Feed the delay while waiting for the signal to appear at the output
	     (if (= source ,last-source)(zdelay (dlocs-path ,l) (if (> ,i (dlocs-end ,l)) 0.0 ,input) 0.0))
	   ;; Signal coming out of the delay line... process it
	   (let* ((,angle (env (aref (dlocs-angle ,l) source)))
		  (,gain (env (aref (dlocs-gain ,l) source)))
		  (,delay (env (aref (dlocs-delay ,l) source)))
		  (,sample (if (= source ,last-source)
			      (zdelay (dlocs-path ,l) (if (> ,i (dlocs-end ,l)) 0.0 ,input) ,delay)
			    (ztap (dlocs-path ,l) ,delay))))
	     (if (/= source 0)
		 (setf ,sample (ppolar (aref (dlocs-filter ,l) source) ,sample)))
	     ;; compute scaling coefficients for the four channels <source outside speakers>
	     (if (or (= 0 (dlocs-only-to ,l))(= (+ source 1)(dlocs-only-to ,l)))
		 (let* ((,scaled-sample (* ,sample ,gain)))
		   (if (or (stereo)(quad))
		       (cond ((< ,angle 0.0)
			      (incf ,d-sample (* ,scaled-sample (array-interp (dlocs-pan-da ,l)
									      (- (+ (dlocs-upper-da ,l) ,angle)(dlocs-upper-cd ,l)))))
			      (incf ,a-sample (* ,scaled-sample (array-interp (dlocs-pan-da ,l)(- ,angle)))))
			     ((<= ,angle (dlocs-upper-ab ,l))
			      (incf ,a-sample (* ,scaled-sample (array-interp (dlocs-pan-ab ,l) ,angle)))
			      (incf ,b-sample (* ,scaled-sample (array-interp (dlocs-pan-ab ,l)(- (dlocs-upper-ab ,l) ,angle)))))
			     ((<= ,angle (dlocs-upper-bc ,l))
			      (incf ,b-sample (* ,scaled-sample (array-interp (dlocs-pan-bc ,l)(- ,angle (dlocs-upper-ab ,l)))))
			      (incf ,c-sample (* ,scaled-sample (array-interp (dlocs-pan-bc ,l)(- (dlocs-upper-bc ,l) ,angle)))))
			     ((<= ,angle (dlocs-upper-cd ,l))
			      (incf ,c-sample (* ,scaled-sample (array-interp (dlocs-pan-cd ,l)(- ,angle (dlocs-upper-bc ,l)))))
			      (incf ,d-sample (* ,scaled-sample (array-interp (dlocs-pan-cd ,l)(- (dlocs-upper-cd ,l) ,angle)))))
			     ((<= ,angle (dlocs-upper-da ,l))
			      (incf ,d-sample (* ,scaled-sample (array-interp (dlocs-pan-da ,l)(- ,angle (dlocs-upper-cd ,l)))))
			      (incf ,a-sample (* ,scaled-sample (array-interp (dlocs-pan-da ,l)(- (dlocs-upper-da ,l) ,angle)))))
			     ((> ,angle (dlocs-upper-da ,l))
			      (incf ,a-sample (* ,scaled-sample (array-interp (dlocs-pan-ab ,l) (- ,angle (dlocs-upper-da ,l)))))
			      (incf ,b-sample (* ,scaled-sample (array-interp (dlocs-pan-ab ,l)
									      (- (dlocs-upper-ab ,l) (- ,angle (dlocs-upper-da ,l)))))))
			     (t (clm-print "error: angle ~f out of bounds~%" ,angle)))
		     (incf ,a-sample ,scaled-sample))
		   ;; compute coefficient for the reverberator channel <source outside speakers>
		   (incf ,r-sample (* ,sample (dlocs-rev ,l)
				      (env (aref (dlocs-rgain ,l) source)))))))))
       ;; send sample to output streams, fold back channels if stereo
       (outa ,i (if (stereo) 
		    (/ (+ ,a-sample ,d-sample) 2)
		  ,a-sample))
       (if (or (stereo)(quad))
	   (outb ,i (if (stereo)
			(/ (+ ,b-sample ,c-sample) 2)
		      ,b-sample)))
       (if (quad)
	   (progn (outc ,i ,c-sample)
		  (outd ,i ,d-sample)))
       (if *reverb* (outa ,i ,r-sample *reverb*)))))


;;;;;;;;;
;;; Paths

;; path in cartesian coordinates
;;   format "a": '((x0 y0 t0)(x1 y1 t1)...(xn yn tn))
;;     t = time of breakpoint (tn is last because it is options)
;;     x = x coordinate of source
;;     y = y coordinate of source
;;   format "b": '(x0 y0 x1 y1...xn yn)
;;     x = x coordinate of source
;;     y = y coordinate of source
;; path in polar coordinates
;;   format "a": '((d0 a0 t0)(d1 a1 t1)...(dn an tn))
;;     t = time of breakpoint (tn is last because it is optional)
;;     d = distance from source to listener
;;     a = angle in degres
;;   format "b": '(d0 a0 d1 a1...dn an)
;;     d = distance from source to listener
;;     a = angle in degres

(defstruct path
  path					; path description
  closed				; t-> closed path; nil->open path
  initial-tangent			; initial tangent for an open path
  final-tangent				; final tangent for an open path
  x y time				; cartesian components of the path (split-path)
  curvature				; curvature of bezier segments
  xc yc					; bezier control points
  (error 0.1)				; error bound for rendering path
  xr yr tr				; rendered cartesian components (render-path)
  xt yt tt)				; transformed cartesian components (transform-path)

(defstruct (polar-path			; makes it possible to differentiate polar paths
	    (:include path)))

(defstruct (literal-path
	    (:include path)))		; makes it possible to differenciate literal paths
					; literal paths are those that are rendered by
					; definition, includes all kinds of geometric
					; paths (for example, spirals)

;;; Define a closed cartesian path

(defmacro make-closed-path (&rest args &key path &allow-other-keys)
  `(if (is-closed ,path)
       (make-path :path ,path ,@args :closed t)
     (error "First and last x/y coordinates of a closed path must be the same.")))

;;; Define a closed polar path

(defmacro make-closed-polar-path (&rest args &key path &allow-other-keys)
  `(if (is-closed ,path)
       (make-polar-path :path ,path ,@args :closed t)
     (error "First and last x/y coordinates of a closed path must be the same.")))

;;; Check to see if a path is closed or not

(defun is-closed (path)
  (if (listp (first path))
      (let* ((start (first path))
	     (end (car (last path))))
	(and (= (first start)(first end))
	     (= (second start)(second end))))
    (and (= (first path)(clm::lastx path))
	 (= (second path)(last path)))))

;; Define geometric paths as types of literal paths
;; a) Spirals

(defun make-spiral-path (&key
			 (start-angle 0)
			 (total-angle 360)
			 (turns nil)
			 (step-angle 10)
			 (initial-distance 20)
			 (final-distance 10)
			 (mod-distance 0)
			 (cycles-turn 1.5))
  (if (< turns 0)
      (setf step-angle (* step-angle -1)))
  (multiple-value-bind (time x y)
      (loop 
	  with segments = (/ (+ start-angle 
				(if turns (* (abs turns) 360) total-angle))
			     (abs step-angle))
	  repeat (+ segments 1)
	  for angle from start-angle by step-angle
	  for radians from (/ (* start-angle 2 pi) 360) by (/ (* step-angle 2 pi) 360)
	  for distance from initial-distance by (/ (- final-distance initial-distance) segments)
	  for time from 0 by (/ segments)
	  for vect = (cis (* (/ angle 360) 2 pi))
	  for dist = (+ distance (* mod-distance (sin (* radians cycles-turn))))
	  collect time into time-r
	  collect (coerce (* dist (imagpart vect)) 'short-float) into x
	  collect (coerce (* dist (realpart vect)) 'short-float) into y
	  finally (return (values time-r x y)))
    (make-literal-path :tr time
		       :xr x
		       :yr y)))

;;; Inquiries into the state of the path

(defun not-split (path)
  (if (literal-path-p path)
      nil
    (null (path-x path))))

(defun not-fitted (path)
  (if (literal-path-p path)
      nil
    (null (path-xc path))))

(defun not-rendered (path)
  (if (literal-path-p path)
      nil
    (null (path-xr path))))

(defun not-transformed (path)
  (null (path-xt path)))

;;; Return the best possible set of x coordinates

(defun get-x-coords (path)
  (if (path-xt path)
      (path-xt path)
    (if (path-xr path)
	(path-xr path)
      (progn
	(render-path path)
	(path-xr path)))))

(defun get-y-coords (path)
  (if (path-yt path)
      (path-yt path)
    (if (path-yr path)
	(path-yr path)
      (progn
	(render-path path)
	(path-yr path)))))

(defun get-t-coords (path)
  (if (path-tt path)
      (path-tt path)
    (if (path-tr path)
	(path-tr path)
      (progn
	(render-path path)
	(path-tr path)))))

;;; Reset different stages of the rendering process

;;; Reset any transformations on the originally rendered path

(defun reset-transformation (path)
  (setf (path-tt path) nil
	(path-xt path) nil
	(path-yt path) nil)
  path)

;;; Reset the rendered path (and any transformations also)

(defun reset-rendering (path)
  (when (not (literal-path-p path))
    (setf (path-tr path) nil
	  (path-xr path) nil
	  (path-yr path) nil))
  (reset-transformation path))

;;; Reset the path split into separate x and y lists

(defun reset-split (path)
  (when (not (literal-path-p path))
    (setf (path-tt path) nil
	  (path-xt path) nil
	  (path-yt path) nil))
  (reset-rendering path))

;;; Set components of a path and reset appropriate part of the rendering process
;;; Rendering steps:
;;;   (split-path path)      split the path into x, y and t lists
;;;   (fit-path path)        calculate the bezier curve control points
;;;   (render-path path)     derive a linear approximation to the bezier segments

;;; Set a new path

(defun set-path (in path)
  (setf (path-path in) path)
  (reset-split in))

;;; Set a new polar path

(defun set-polar-path (in path)
  (when (polar-path-p in)
    (setf (path-path in) path)
    (reset-split in)))

;;; Set a new path curvature

(defun set-path-curvature (in curvature)
  (if (literal-path-p in)
      (error "a literal path does not have a rendering curvature"))
  (when curvature
    (setf (path-curvature in) curvature)
    (reset-rendering in)))

;;; Set a new path rendering error bound

(defun set-path-error (in error)
  (if (literal-path-p in)
      (error "a literal path does not have a rendering error"))
  (when error
    (setf (path-error in) error)
    (reset-rendering in)))

;;; Split a path into its time, x and y components

(defun split-path (path)
  (if (literal-path-p path)
      (return-from split-path path)
    (if (polar-path-p path)
	;; a path specified in polar coordinates
	(if (listp (first (path-path path)))
	    ;; decode a list of lists into x:y:t components
	    (loop for number from 0 
		for point in (polar-path-path path)
		for dist = (first point)
		for angle = (second point)
		for vect = (cis (* (/ angle 360) 2 pi))
		collect (third point) into time
		collect (coerce (* dist (imagpart vect)) 'short-float) into x
		collect (coerce (* dist (realpart vect)) 'short-float) into y
		finally (setf (path-x path) x
			      (path-y path) y 
			      (path-time path) time))
	  ;; decode a list of x:y components
	  (loop for number from 0 
	      for dist in (polar-path-path path) by #'cddr
	      for angle in (cdr (polar-path-path path)) by #'cddr
	      for vect = (cis (* (/ angle 360) 2 pi))
	      collect (coerce (* dist (imagpart vect)) 'short-float) into x
	      collect (coerce (* dist (realpart vect)) 'short-float) into y
	      finally (setf (polar-path-x path) x
			    (polar-path-y path) y
			    (polar-path-time path)(loop repeat (length x) collect nil))))
      (if (path-p path)
	  ;; has to be a cartesian path
	  (if (listp (first (path-path path)))
	      ;; decode a list of lists into x:y:t components
	      (loop for number from 0
		  for point in (path-path path) 
		  collect (first point) into x
		  collect (second point) into y
		  collect (third point) into time
		  finally (setf (path-x path) x
				(path-y path) y
				(path-time path) time))
	    ;; decode a list of x:y components
	    (loop for number from 0
		for xp in (path-path path) by #'cddr
		for yp in (cdr (path-path path)) by #'cddr
		collect xp into x
		collect yp into y
		finally (setf (path-x path) x
			      (path-y path) y
			      (polar-path-time path)(loop repeat (length x) collect nil))))
	;; path it was not a path!
	(format t "~s not a path or polar-path structure~%" path))))
  (loop 
      with length = (length (path-time path))
      with count-nil = 0
      for time in (path-time path) do
	(if (null time)
	    (incf count-nil))
      finally
	(if (and (/= count-nil 0)
		 (/= count-nil length))
	    (if (not (and (first (path-time path))
			  (car (last (path-time path)))))
		(error "path with nil times need to have first and last times specified"))))
  path)

;;; Bezier curve fitting auxiliary functions

(defparameter gtab nil)
(defparameter ftab nil)
(defparameter ak-even nil)
(defparameter ak-odd nil)
(defconstant maxcoeff 8)

(defun make-a-even()
  (flet ((g(m)
	   (if (null gtab)
	       (progn
		 (setf gtab (make-array maxcoeff)
		       (aref gtab 0) 1
		       (aref gtab 1) -4)
		 (loop for i from 2 to (- maxcoeff 1) do
		   (setf (aref gtab i)(- (* -4 (aref gtab (- i 1)))
					 (aref gtab (- i 2)))))))
	   (aref gtab m)))
    (setf ak-even (make-array (- maxcoeff 1)))
    (loop for m from 1 to (- maxcoeff 1) do
      (setf (aref ak-even (- m 1))(make-array m))
      (loop for k from 1 to m do
	(setf (aref (aref ak-even (- m 1))(- k 1))
	  (float (/ (- (g (- m k)))(g m))))))))

(defun make-a-odd()
  (flet ((f(m)
	   (if (null ftab)
	       (progn
		 (setf ftab (make-array maxcoeff)
		       (aref ftab 0) 1
		       (aref ftab 1) -3)
		 (loop for i from 2 to (- maxcoeff 1) do
		   (setf (aref ftab i)(- (* -4 (aref ftab (- i 1)))
					 (aref ftab (- i 2)))))))
	   (aref ftab m)))
    (setf ak-odd (make-array (- maxcoeff 1)))
    (loop for m from 1 to (- maxcoeff 1) do
      (setf (aref ak-odd (- m 1))(make-array m))
      (loop for k from 1 to m do
	(setf (aref (aref ak-odd (- m 1))(- k 1))
	  (float (/ (- (f (- m k)))(f m))))))))

(defun a(k n)
  (if (oddp (min (+ (* maxcoeff 2) 1) n))
      (progn
	(if (null ak-odd)(make-a-odd))
	(aref (aref ak-odd (/ (- n 3) 2))(- k 1)))
    (progn
      (if (null ak-even)(make-a-even))
      (aref (aref ak-even (/ (- n 4) 2))(- k 1)))))

(defun ac(k n)
  (setf n (min n maxcoeff))
  (if (null ak-even)(make-a-even))
  (aref (aref ak-even (- n 2))(- k 1)))

;;; Calculate bezier difference vectors for the given path

(defun calculate-fit-closed (path) 
  (let* ((n (- (length (path-x path)) 1))
	 (m (/ (- n (if (oddp n) 3 4)) 2))
	 (p (make-array `( 2 ,(+ n 1)) 
			:initial-contents (list (path-x path)
						(path-y path))))
	 (d (make-array `( 2 ,n) :initial-element 0.0)))
    (flet ((ref(z j i)
	     (cond ((> i (- n 1))(aref z j (- i n)))
		   ((< i 0)(aref z j (+ i n)))
		   (t (aref z j i)))))
      (loop for i from 0 below n do
	    (loop for k from 1 to m do 
		  (incf (aref d 0 i)
			(* (a k n)
			   (- (ref p 0 (+ i k))
			      (ref p 0 (- i k)))))
		  (incf (aref d 1 i)
			(* (a k n)
			   (- (ref p 1 (+ i k))
			      (ref p 1 (- i k)))))))
      (if (path-curvature path)
	  (loop for i from 0 to (- n 1)
	      for curve = (path-curvature path) do
		(setf (aref d 0 i)(* (aref d 0 i) curve)
		      (aref d 1 i)(* (aref d 1 i) curve))))
      (values (- n 1) p d))))


(defun calculate-fit-open (path) 
  (let* ((n (- (length (path-x path)) 1))
	 (m (- n 1))
	 (p (make-array `( 2 ,(+ n 1)) 
			:initial-contents (list (path-x path)
						(path-y path))))
	 (d (make-array `( 2 ,(+ n 1)) 
			:initial-element 0.0)))
    (flet ((ref(z j i)
	     (cond ((> i n)(aref z j (- i n)))
		   ((< i 0)(aref z j (+ i n)))
		   ((= i n)(- (aref z j n)(aref d j n)))
		   ((= i 0)(+ (aref z j 0)(aref d j 0)))
		   (t (aref z j i)))))
      (if (path-initial-tangent path)
	  (setf (aref d 0 0)(first (path-initial-tangent path))
		(aref d 1 0)(second (path-initial-tangent path)))
	(setf (aref d 0 0) 0.0
	      (aref d 1 0) 0.0))
      (if (path-final-tangent path)
	  (setf (aref d 0 n)(first (path-final-tangent path))
		(aref d 1 n)(second (path-final-tangent path)))
	(setf (aref d 0 n) 0.0
	      (aref d 1 n) 0.0))
      (loop for i from 1 below n do
	    (loop for k from 1 to (min (- maxcoeff 1) m) do 
		  (incf (aref d 0 i)
			(* (ac k n)
			   (- (ref p 0 (+ i k))
			      (ref p 0 (- i k)))))
		  (incf (aref d 1 i)
			(* (ac k n)
			   (- (ref p 1 (+ i k))
			      (ref p 1 (- i k)))))))
      (if (path-curvature path)
	  (loop for i from 0 to n 
	      for curve = (path-curvature path) do
		(setf (aref d 0 i)(* (aref d 0 i) curve)
		      (aref d 1 i)(* (aref d 1 i) curve))))
      (values n p d))))

;;; Calculate Bezier control points for the given path

(defun fit-path (path)
  (if (literal-path-p path)
      (return-from fit-path path))
  (if (not-split path)
      (split-path path))
  (if (not (path-closed path))
      ;;
      ;; for an open path
      ;;
      (if (> (length (path-x path)) 2)
	  (multiple-value-bind (n p d)
	      (calculate-fit-open path)
	    (loop 
		for i from 0 to (- n 1)
		collect 
		  (list (aref p 0 i)
			(+ (aref p 0 i)(aref d 0 i))
			(- (aref p 0 (+ i 1))(aref d 0 (+ i 1)))
			(aref p 0 (+ i 1))) into xc
		collect 
		  (list (aref p 1 i)
			(+ (aref p 1 i)(aref d 1 i))
			(- (aref p 1 (+ i 1))(aref d 1 (+ i 1)))
			(aref p 1 (+ i 1))) into yc
		finally (setf (path-xc path) xc
			      (path-yc path) yc
			      ;; as fit changed unrender path
			      (path-xr path) nil
			      (path-yr path) nil
			      (path-tr path) nil)))
	(let* ((x1 (first (path-x path)))
	       (x2 (second (path-x path)))
	       (y1 (first (path-y path)))
	       (y2 (second (path-y path))))
	  (warn "[fit-path:open-path] only two points in path")
	  (setf (path-xc path)(list (list x1 x1 x2 x2))
		(path-yc path)(list (list y1 y1 y2 y2))
		;; as fit changed unrender path
		(path-xr path) nil
		(path-yr path) nil
		(path-tr path) nil)))
    ;;
    ;; for a closed path
    ;;
    (if (> (length (path-x path)) 4)
	(multiple-value-bind (n p d)
	    (calculate-fit-closed path)
	  (loop 
	      for i from 0 to (- n 1)
	      collect 
		(list (aref p 0 i)
		      (+ (aref p 0 i)(aref d 0 i))
		      (- (aref p 0 (+ i 1))(aref d 0 (+ i 1)))
		      (aref p 0 (+ i 1))) into xc
	      collect 
		(list (aref p 1 i)
		      (+ (aref p 1 i)(aref d 1 i))
		      (- (aref p 1 (+ i 1))(aref d 1 (+ i 1)))
		      (aref p 1 (+ i 1))) into yc
	      finally 
		(setf 
		  (path-xc path)
		  (nconc xc (list (list (aref p 0 n)
					(+ (aref p 0 n)(aref d 0 n))
					(- (aref p 0 0)(aref d 0 0))
					(aref p 0 0))))
		  (path-yc path)
		  (nconc yc (list (list (aref p 1 n)
					(+ (aref p 1 n)(aref d 1 n))
					(- (aref p 1 0)(aref d 1 0))
					(aref p 1 0))))
		  ;; as fit changed unrender path
		  (path-xr path) nil
		  (path-yr path) nil
		  (path-tr path) nil)))
      (loop 
	  for x1 in (path-x path)
          for x2 in (cdr (path-x path))
          for y1 in (path-y path)
          for y2 in (cdr (path-y path))
	  collect (list x1 x1 x2 x2) into xc
	  collect (list y1 y1 y2 y2) into yc
	  finally (warn "[fit-path:closed-path] not enough points to do bezier fit (~s points)"
			(length (path-x path)))
	  finally (setf (path-xc path) xc
			(path-yc path) yc
			;; as fit changed unrender path
			(path-xr path) nil
			(path-yr path) nil
			(path-tr path) nil))))
  path)

;;; Transform a Bezier control point fit to a linear segment approximation

(defun render-path (path)
  (if (literal-path-p path)
      (return-from render-path path))
  (if (not-fitted path)
      (fit-path path))
  (labels (;;
	   ;; Evaluate a point at parameter u in bezier segment
	   ;;
	   (bezier-point (u c)
	     (let* ((u1 (- 1 u))
		    (cr (make-array '(2 3) :element-type 'float)))
	       (loop for j from 0 to 2 do
		 (setf (aref cr 0 j)(+ (* u1 (aref c 0 j))(* u (aref c 0 (+ j 1))))
		       (aref cr 1 j)(+ (* u1 (aref c 1 j))(* u (aref c 1 (+ j 1))))))
	       (loop for i from 1 downto 0 do
		 (loop for j from 0 to i do
		   (setf (aref cr 0 j)(+ (* u1 (aref cr 0 j))(* u (aref cr 0 (+ j 1))))
			 (aref cr 1 j)(+ (* u1 (aref cr 1 j))(* u (aref cr 1 (+ j 1)))))))
	       (values (aref cr 0 0)
		       (aref cr 1 0))))
	   ;;
	   ;; Create a linear segment rendering of the bezier segment
	   ;;
	   (berny (xl yl xh yh ul u uh c err)
	     (multiple-value-bind (x y)
		 (bezier-point u c)
	       (multiple-value-bind (xn yn)
		   (nearest-point xl yl xh yh x y)
		 (if (> (distance (- xn x)(- yn y)) err)
		     (multiple-value-bind (xi yi)
			 (berny xl yl x y ul (/ (+ ul u) 2) u c err)
		       (multiple-value-bind (xj yj)
			   (berny x y xh yh u (/ (+ u uh) 2) uh c err)
			 (values (nconc xi (list x) xj)
				 (nconc yi (list y) yj))))
		   (values nil nil)))))
	   ;;
	   ;; Calculate the total distance travelled in a segment
	   ;;
	   (segment-distance (xs ys)
	     (loop 
		 for xsi in xs
		 for ysi in ys
		 for xsf in (cdr xs)
		 for ysf in (cdr ys)
		 sum (distance (- xsf xsi)(- ysf ysi)))))
    ;;
    ;; Render each segment of the fitted path
    ;;
    (if (not (first (path-time path)))
	(setf (first (path-time path)) 0
	      (car (last (path-time path))) 1))
    (loop 
	;; accumulate rendered path
	with xr
	with yr
	with tr
	;; initial and final coords
	with xi 
	with yi 
	with ti 
	with xf 
	with yf 
	with tf
	;; iterate over the control points
	for x-bz in (path-xc path)
	for y-bz in (path-yc path)
	;; iterate over the times
	for ti-bz in (path-time path)
	for tf-bz in (cdr (path-time path))
	;; initial and final values
	for xi-bz = (first x-bz)
	for xf-bz = (car (last x-bz))
	for yi-bz = (first y-bz)
	for yf-bz = (car (last y-bz))
	do
	  ;; Calculate lists of x and y components and add to current segment
	  ;;
	  (if (and (null xr)(null yr))
	      (setf ti ti-bz
		    xi xi-bz
		    yi yi-bz))
	  (setf tf tf-bz
		xf xf-bz
		yf yf-bz)
	  ;; approximate the bezier curve with linear segments
	  (multiple-value-bind (xs ys)
	      (berny xi-bz yi-bz xf-bz yf-bz 0.0 0.5 1.0 
		     (make-array '(2 4) :initial-contents (list x-bz y-bz))
		     (path-error path))
	    (setf xr `(,.xr ,xi-bz ,.xs)
		  yr `(,.yr ,yi-bz ,.ys)
		  ;; accumulate intermediate unknown times as nils
		  tr `(,.tr ,ti-bz ,.(loop repeat (length xs) collect nil))))
	finally 
	  ;; add the last point
	  (setf (path-xr path) `(,.xr ,xf-bz)
		(path-yr path) `(,.yr ,yf-bz)
		(path-tr path) `(,.tr ,tf-bz)))
    ;;
    ;; Calculate average velocity for each timed segment
    ;;
    (let* ((avg-v '())
	   (init-v '())
	   (seg-dist '())
	   vi)
      (loop
	  ;; initialize with start of path
	  with xseg = (list (car (path-xr path)))
	  with yseg = (list (car (path-yr path)))
	  with tseg = (list (car (path-tr path)))
	  ;; iterate over rest of coordinates in the path
	  for x in (cdr (path-xr path))
	  for y in (cdr (path-yr path))
	  for time in (cdr (path-tr path))
	  do
	    ;; add the point to the segment
	    (setf xseg `(,.xseg ,x)
		  yseg `(,.yseg ,y)
		  tseg `(,.tseg ,time))
	    ;; process the segment when it gets to an end time
	    (when time
	      (let* ((dist (segment-distance xseg yseg)))
		;; acummulate the average velocity
		(setf avg-v `(,.avg-v ,(/ dist (- (car (last tseg))(car tseg))))
		      seg-dist `(,.seg-dist ,dist))
		;; start the next timed segment
		(setf xseg (list x)
		      yseg (list y)
		      tseg (list time)))))
      ;;
      ;; Determine starting velocities for each segment
      ;;
      (labels (;;
	       ;; Calculate starting velocities for a fixed end chain
	       ;;
	       (fixed-end (vs vf)
		 (loop 
		     with vi
		     for chain on vs
		     collect
		       (loop 
			   for link in chain
			   for mult = 1 then (* mult -1)
			   sum (* 2 link mult) into start
			   finally 
			     (setf vi `(,.vi ,(- start (* vf mult)))))
		     finally
		       (return vi)))
	       ;;
	       ;; Calculate starting velocities for a fixed start chain
	       ;;
	       (fixed-start (vs vi)
		 (let* ((vs (reverse vs)))
		   `(,vi ,.(nbutlast (reverse (fixed-end vs vi))))))
	       ;;
	       ;; Calculate starting velocities for a fixed end and start chanin
	       ;;
	       (fixed-start-end (vs vi vf)
		 (if (> (length vs) 1)
		     ;; create a split segment at the beginning of the chain
		     (let* ((vseg (fixed-end (rest vs) vf)))
		       (list (list vi nil (first vseg)) vseg))
		   ;; one segment means it is just a split segment
		   (list (list vi nil vf))))
	       ;;
	       ;; Final velocity of a chain
	       ;;
	       (final-velocity (vs vi)
		 (let* ((vs (reverse vs)))
		   (loop 
		       for each in vs
		       for mult = 1 then (* mult -1)
		       sum (* 2 each mult) into final
		       finally 
			 (return (- final (* mult vi)))))))
	(loop 
	    with state
	    with starts 
	    with pending
	    with prev
	    with next
	    with this
	    with final
	    with stay
	    with split
	    for rest on avg-v
	    do
	      ;; initialize velocities
	      (setf this (first rest)
		    next (second rest))
	      ;; decide if this segment is constant velocity or has to be split
	      ;; these criteria are very simplistic...
	      (if (not prev)
		  (if (not next)
		      ;; one and only one segment
		      (setf stay t
			    split nil)
		    ;; first segment: cannot be split
		    (setf stay (< this next)
			  split nil))
		;; subsequent segments
		(if next
		    (setf stay (and (> prev this)(> next this))
			  split (if (and (< prev this)(> this next))
				    (list prev this next)
				  nil))
		  ;; last segment: cannot be split
		  (setf stay (> prev this)
			split nil)))
	      ;; process a segment
	      (if (not prev)
		  ;; first segment
		  (if stay
		      ;;
		      ;; if it is a contant velocity segment then define the starting
		      ;; velocity for it and start a fixed start segment.
		      ;; 
		      (setf init-v `(,.init-v ,(list this))
			    state :fixed-start
			    starts this
			    final this
			    pending nil)
		    ;;
		    ;; otherwise nothing is defined and a fixed end chain is started
		    ;;
		    (setf state :fixed-end
			  starts this
			  final nil
			  pending (list this)))
		;; subsequent segments
		(progn
		  (if (eq state :fixed-end)
		      (if stay
			  ;;
			  ;; fixed end chain ends in a constant segment, define velocity
			  ;; for chain ending as the constant segment average velocity.
			  ;; Start a new fixed start chain. 
			  ;;
			  (progn
			    (if pending
				(setf init-v `(,.init-v ,.(fixed-end pending this))))
			    (setf init-v `(,.init-v ,(list this))
				  state :fixed-start
				  starts this
				  final this
				  pending nil))
			(if split
			    ;;
			    ;; fixed end chain ends in a split segment. I don't this can
			    ;; actually happen. We define the starting velocity of the split
			    ;; to be its average velocity and the ending to be the average
			    ;; velocity of the next segment (FIXME: what if no next segment?)
			    ;;
			    (progn
			      (if pending
				  (progn
				    (format t "# fixed-end [split] ~s = ~s~%" pending (fixed-end pending this))
				    (setf init-v `(,.init-v ,.(fixed-end pending this)))))
			      (format t "# split vi=~s vf=~s~%" starts next)
			      (setf init-v `(,.init-v ,(list starts this next)))
			      (setf state :fixed-start
				    starts next
				    final next
				    pending nil))
			  ;;
			  ;; else accumulate segment into fixed end chain
			  ;;
			  (progn
			    (setf pending `(,.pending ,this)))))
		    (if (eq state :fixed-start)
			(if stay
			    ;;
			    ;; fixed start chain ends in a constant segment. That transforms
			    ;; the chain into a fixed end and start, introducing a split in
			    ;; the first segment of the chain to accomodate the fixed bounds. 
			    ;;
			    (progn
			      (if pending
				  (setf init-v `(,.init-v ,.(fixed-start-end pending starts this))))
			      (setf init-v `(,.init-v ,(list this))
				    state :fixed-start
				    starts this
				    final this
				    pending nil))
			  (if split
			      ;;
			      ;; fixed start chain ends in a split segment. The start of the split
			      ;; is defined by the ending velocity of the fixed start chain. If
			      ;; no pending segments then the last set "final" velocity defines
			      ;; the starting velocity for the split. Ending velocity of the split
			      ;; is the average velocity of the next segment, quite arbitrary. 
			      ;;
			      (progn
				(if pending
				    (setf init-v `(,.init-v ,.(fixed-start pending starts))
					  final (final-velocity pending starts)))
				(if (not final) 
				    (error "No pending final velocity."))
				(setf init-v `(,.init-v ,(list final this next))
				      state :fixed-start
				      starts next
				      final next
				      pending nil))
			    ;;
			    ;; else acumulate segment into fixed start chain
			    ;;
			    (progn
			      (setf pending `(,.pending ,this)))))))))
	      ;; remember velocity for next segment
	      (setf prev this)
	    finally
	      (if pending
		  ;;
		  ;; process pending segments when loop terminates. We assume we are in a
		  ;; fixed start chain, can't see how we could be in a fixed end chain. 
		  ;;
		  (if (eq state :fixed-start)
		      (setf init-v `(,.init-v ,.(fixed-start pending starts)))
		    (error "Pending segments [~s] and not a fixed start velocity chain~%" vs)))))
      (flet (;;
	     ;; Calculate times within a segment
	     ;;
	     (calc-times (xs ys ti tf vi vf)
	       (loop 
		   with a = (/ (- (segment-distance xs ys)
				  (* vi (- tf ti)))
			       (* (- tf ti)(- tf ti)))
		   with time = 0
		   with dist = 0
		   with times = (list ti)
		   for xi in xs
		   for yi in ys
		   for xf in (cdr xs)
		   for yf in (cdr ys)
		   do 
		     (incf dist (distance (- xf xi)(- yf yi)))
		     (setf time (+ ti (if (< (abs a) 0.001)
					  (/ dist vi)
					(/ (- (sqrt (+ (* vi vi)
						       (* 4 a dist)))
					      vi)
					   (* 2 a)))))
		     (setf times `(,.times ,time))		 
		   finally 
		     (return times))))
	;;
	;; Calculate all the intermediate times assuming constant acceleration
	;;
	(loop 
	    with split
	    ;; these values are iterated over once per time-limited segment
	    with dist = seg-dist
	    with vm-rest = avg-v
	    with vi-rest = init-v
	    with vi
	    with vf
	    ;; calculated path times
	    with times
	    ;; set up first point
	    with xseg = (list (car (path-xr path)))
	    with yseg = (list (car (path-yr path)))
	    with tseg = (list (car (path-tr path)))
	    ;; iterate over rest of coords
	    for x in (cdr (path-xr path))
	    for y in (cdr (path-yr path))
	    for time in (cdr (path-tr path))
	    do
	      ;; add the point to the segment
	      (setf xseg `(,.xseg ,x)
		    yseg `(,.yseg ,y)
		    tseg `(,.tseg ,time))
	      ;; process the segment when it gets to an end time
	      (when time
		;; determine type of segment and initial and final velocities
		(if (listp (first vi-rest))
		    (if (= (length (first vi-rest)) 1)
			;; constant velocity segment
			(setf split nil
			      vi (first (first vi-rest))
			      vf (first (first vi-rest)))
		      ;; split segment
		      (setf split t
			    vi (first (first vi-rest))
			    vf (third (first vi-rest))))
		  ;; normal segment
		  (setf split nil
			vi (first vi-rest)
			vf (if (second vi-rest)
			       (second vi-rest)
			     (- (* 2 (first vm-rest))
				(first vi-rest)))))
		(if split
		    ;;
		    ;; Calculate times for a split segment
		    ;;
		    (loop
			with xi = (first xseg)
			with xf = (first (last xseg))
			with yi = (first yseg)
			with yf = (first (last yseg))
			with ti = (first tseg)
			with tf = (first (last tseg))
			with ds = (first dist)
			;; first part of split segment
			with xsi
			with ysi
			with tsi
			with dsi
			;; second part of split segment
			with xsf
			with ysf
			with tsf
			with dsf
			;; iterate over the coordinates
			for x in xseg
			for y in yseg
			for time in tseg
			for xe in (cdr xseg)
			for ye in (cdr yseg)
			;; incremental distance travelled in segment
			sum (distance (- xe x)(- ye y)) into sofar
			do
			  (if (< sofar (/ ds 2))
			      ;; collect first part
			      (setf xsi `(,.xsi ,x)
				    ysi `(,.ysi ,y)
				    tsi `(,.tsi ,time))
			    ;; collect second part
			    (setf xsf `(,.xsf ,x)
				  ysf `(,.ysf ,y)
				  tsf `(,.tsf ,time)))
			finally
			  (let* ((tx)
				 (vx))
			    (setf xsi `(,.xsi ,(first xsf))
				  ysi `(,.ysi ,(first ysf))
				  dsi (segment-distance xsi ysi)
				  xsf `(,.xsf ,xf)
				  ysf `(,.ysf ,yf)
				  tsf `(,.tsf ,tf)
				  dsf (segment-distance xsf ysf)
				  tx (/ (- (+ tf ti)
					   (/ (sqrt (+ (* 4 (- (* (+ dsf dsi)(+ dsf dsi))
							       (* (- dsf dsi)
								  (- (* vf (- tf ti))
								     (* vi (- tf ti))))))
						       (* (- tf ti)(- tf ti)
							  (- vf vi)(- vf vi))))
					      (- vi vf))
					   (/ dsf (- vf vi) 0.5)
					   (/ dsi (- vf vi) 0.5)) 2)
				  vx (+ vi (* 2 (/ (- dsi (* vi (- tx ti)))
						   (- tx ti))))
				  tsi `(,.tsi ,tx)
				  (first tsf) tx)
			    (setf times `(,.times ,.(cdr (calc-times xsi ysi ti tx vi vx))))
			    (setf times `(,.times ,.(cdr (calc-times xsf ysf tx tf vx vf))))))
		  ;;
		  ;; Calculate times for a normal segment
		  ;;
		  (setf times `(,.times 
				,.(cdr (calc-times xseg yseg
						   (first tseg)(first (last tseg))
						   vi vf)))))
		;;
		;; Set up values for next segment
		;;
		(setf dist (cdr dist)
		      vm-rest (cdr vm-rest)
		      vi-rest (cdr vi-rest)
		      ;; init coordinate lists with first point
		      xseg (list x)
		      yseg (list y)
		      tseg (list time)))
	    finally 
	      ;;
	      ;; Set times of path to rendered times
	      ;;
	      (setf (path-tr path) `(,(first (path-time path)) ,.times)))))))


#|

(0)  x(t) = xi + vi t + a t^2

At the end of the segment:
(1)  xf = xi + vi (tf - ti) + a (tf - ti)^2

From there we get the acceleration to use (assumed constant):
(2)  a = ((xf - xi) - vi (tf - ti)) / (tf - ti)^2

Velocity as a function of time (first equation's derivative)
(3)  v(t) = vi + 2 a t

The ending velocity will be:
(4)  vf = vi + 2 a (tf - ti)

We need to find time as a function of displacement,
the general solution to (0) is:

(5)  t = ti + (- vi +- sqrt(vi^2 - 4 a (xi - x))) / 2 a
(5a) t = ti + (- vi +- sqrt(vi^2 + 4 a (x - xi))) / 2 a

If a=0 (constant velocity segment)

(6)  xf = xi + vi (tf - ti)
(7)  (tf - ti) = (xf -xi) / vi

|#

      #|
      Now comes the fun part. Segments that have average velocity
      lower than both adjoining segments are set to be constant
      velocity segments. If we were to allow acceleration in those
      segments then to mantain the average velocity we could end up
      with negative velocities which is not good. If a segment is not
      constant velocity then we have to accelerate or deccelerate to
      the starting velocity of the next segment. If a segment is
      between two lower average velocity segments then we have to
      split the segment in two, accelerate in the first part and
      decelerate during the second.

      Using simple algebra we arrive at the following expression for
      calculating the initial velocity of each segment:

      vi = 2*vm1 - 2*vm2 + 2*vm3 ... +/- vmn
      
      where 
        vi = initial velocity for the segment
	vm1 = average velocity for current segment
	vin = initial velocity for last considered segment
	
      evaluation of the series stops when we reach the end of the last
      segment or when we reach a segment that is constant velocity or
      has to be split.

      when deciding whether a segment is constant velocity we take
      into account its relative velocity with respect to the two
      adjoining segments:

      first segment:
      -   c   h
      -   .   l
      entermediate segments:
      lc  s   lc
      l   .   h
      lc  .   h
      h   .   l
      hc  .   l       if "l" becomes constant then we're in trouble. 
      h   c   h
      last segment:
      h   c   -  
      l   .   -

      where 
        c=constant velocity
        l=lower velocity
        h=higher velocity
        s=split in two
      
      |#
#|

Test cases:
One segment:
(setf xxx (make-path :path '((10 10 0)(0 6)(-6 0)(-8 -4)(5 -10)(7 0)(3 7 1)) :error 0.003))

Two timed segments:
-- constant segment + fixed start
(setf xxx (make-path :path '((10 10 0)(0 6)(-6 0)(-8 -4 0.5)(5 -10)(7 0)(3 7 1)) :error 0.003))
-- fixed end + constant segment
(setf xxx (make-path :path '((10 10 0)(0 6)(-6 0)(-8 -4 0.3)(5 -10)(7 0)(3 7 1)) :error 0.003))

Several timed segments:
-- fixed-end:const:fixed-start:fixed-start:fixed-start
(setf xxx (make-path :path '((10 10 0)(0 6)(-6 0 0.2)(-8 -4 0.5)(5 -10 0.8)(7 0 0.94)(3 7 1)) :error 0.003))
-- fixed-end:fixed-end:fixed-end:const:fixed-start
(setf xxx (make-path :path '((10 10 0)(0 6)(-6 0 0.3)(-8 -4 0.4)(5 -10 0.7)(7 0 0.95)(3 7 1)) :error 0.003))



;;; (setf xxx (make-path :path '((10 10 0)(0 6)(-6 0 0.6)(-8 -4)(5 -10)(7 0)(3 7 1))))

;;; (setf xxx (make-path :path '((10 10 0)(0 3 0.5)(-10 10 1))))

;;; (setf xxx (make-path :path '((10 10 0)(0 6)(-6 0 0.5)(-8 -4)(5 -10)(7 0)(3 7 1))))

;;; (setf xxx (make-path :path '((10 10 0)(0 6)(-6 0 0.3)(-8 -4 0.4)(5 -10 0.6)(7 0)(3 7 1)) :error 0.02))

;;; (setf xxx (make-path :path '((10 10 0)(0 6)(-6 0 0.35)(-8 -4 0.5)(5 -10 0.76)(7 0 0.8)(3 7 1)) :error 0.01))

;;; (setf xxx (make-path :path '((10 10 0)(0 6)(-6 0 0.35)(-8 -4 0.5)(5 -10 0.76)(7 0 0.88)(3 7 1)) :error 0.01))

|#

#|

No split (0...1)
  (setf xxx (make-path :path '((10 10 0)(0 6)(-6 0)(-8 -4)(5 -10)(7 0)(3 7 1)) :error 0.005))

|#

;;; Transform a rendered path using scaling, translation and rotation 

(defun transform-path (path &key
			    scaling
			    translation
			    rotation
			    rotation-center)
  (if (not-rendered path)
      (render-path path))
  (if (or scaling translation rotation)
      (let* ((scaling (if scaling scaling '(1.0 1.0)))
	     (translation (if translation translation '(0.0 0.0)))
	     (rotation (if rotation (* 2 pi (/ rotation 360)) 0.0))
	     (rotation-center (if rotation-center rotation-center '(0.0 0.0)))
	     (xcoords (get-x-coords path))
	     (ycoords (get-y-coords path))
	     (tcoords (get-t-coords path)))
	(loop for x in xcoords
	    for y in ycoords
	    for xc = (- x (first rotation-center))
	    for yc = (- y (second rotation-center))
	    for vector = (cis (+ (atan xc yc) rotation))
	    for distance = (distance xc yc)
	    collect (coerce (+ (first rotation-center)
			       (first translation)
			       (* distance (first scaling)(imagpart vector))) 
			    'short-float) into xtransf
	    collect (coerce (+ (second rotation-center)
			       (second translation)
			       (* distance (second scaling)(realpart vector))) 
			    'short-float) into ytransf
	    finally (setf (path-tt path) tcoords
			  (path-xt path) xtransf
			  (path-yt path) ytransf)))
    ;; if there's no transformation just copy the rendered path
    (setf (path-tt path)(copy-list (path-tr path))
	  (path-xt path)(copy-list (path-xr path))
	  (path-yt path)(copy-list (path-yr path))))
  path)

;;; Scale a path

(defun scale-path (path scaling)
  (transform-path path :scaling scaling))

;;; Translate a path

(defun translate-path (path translation)
  (transform-path path :translation translation))

;;; Rotate a path

(defun rotate-path (path rotation &optional rotation-center)
  (transform-path path :rotation rotation :rotation-center rotation-center))

;;; Mirror a path around an axis

(defun mirror-path (path &key (axis 'y)
			      (around 0))
  (if (not-transformed path)
      (transform-path path))
  (if (equal axis 'y)
      (setf (path-xt path)
	(loop for x in (path-xt path) collect (- around x)))
    (setf (path-yt path)
      (loop for y in (path-yt path) collect (- around y))))
  path)

;;; Change the times of the rendered envelope so that the velocity is constant

(defun constant-velocity (path)
  (if (not (path-xr path))
      (render-path path))
  (reset-transformation path)
  (let* ((xcoords (get-x-coords path))
	 (ycoords (get-y-coords path))
	 (tcoords (get-t-coords path))
	 (total-distance (loop 
			     for x1 in xcoords
			     for y1 in ycoords
			     for x2 in (cdr xcoords)
			     for y2 in (cdr ycoords)
			     sum (distance (- x2 x1)(- y2 y1))))
	 (start-time (first tcoords))
	 (end-time (first (last tcoords)))
	 (total-time (- end-time start-time))
	 (velocity (/ total-distance total-time)))
    (loop 
	for xp in xcoords
	for yp in ycoords
	for time in (cdr tcoords)
	for x in (cdr xcoords)
	for y in (cdr ycoords)
	sum (distance (- x xp)(- y yp)) into dist
	collect (/ dist velocity) into now
	finally (setf (path-tr path)(nconc (list start-time) now)
		      (path-xt path)(copy-list (path-xr path))
		      (path-yt path)(copy-list (path-yr path)))))
  path)

;;; Return a list with the xy coordinates of the trajectory

(defun trajectory (path)
  (if (not (path-xr path))
      (render-path path))
  (let* ((xcoords (get-x-coords path))
	 (ycoords (get-y-coords path)))
    (loop 
	for x in xcoords
        for y in ycoords
	collect x collect y)))

;;; Return a list with the xyt coordinates of the trajectory

(defun 3d-trajectory (path)
  (if (not (path-xr path))
      (render-path path))
  (let* ((xcoords (get-x-coords path))
	 (ycoords (get-y-coords path))
	 (zcoords (get-t-coords path)))
    (loop 
	for x in xcoords
        for y in ycoords
	for z in zcoords
	collect x collect y collect z)))

;;; Return a list with the velocity values for a path

(defun velocity (path)
  (if (not (path-xr path))
      (render-path path))
  (let* ((xcoords (get-x-coords path))
	 (ycoords (get-y-coords path))
	 (tcoords (get-t-coords path)))
    (loop 
	for timei in tcoords
	for timej in (cdr tcoords)
	for xi in xcoords
	for yi in ycoords
	for xj in (cdr xcoords)
	for yj in (cdr ycoords)
	for v = (/ (distance (- xj xi)(- yj yi))
		   (- timej timei))
	collect timei
	collect v)))

;;; Return a list with the doppler shift of a path

(defun doppler (path)
  (if (not (path-xr path))
      (render-path path))
  (let* ((xcoords (get-x-coords path))
	 (ycoords (get-y-coords path))
	 (tcoords (get-t-coords path)))
    (loop 
	for timei in tcoords
	for timej in (cdr tcoords)
	for xi in xcoords
	for yi in ycoords
	for xj in (cdr xcoords)
	for yj in (cdr ycoords)
	collect (/ (+ timej timei) 2)
	collect (- (/ (- (distance xj yj)(distance xi yi))
		      (- timej timei))))))

;;; Return a list with the acceleration of a path

(defun acceleration (path)
  (let* ((v (velocity path)))
    (loop 
	for timei in v by #'cddr
	for vi in (cdr v) by #'cddr
	for timej in (cddr v) by #'cddr
	for vj in (cdddr v) by #'cddr
	for am = (/ (- vj vi)
		    (- timej timei))
	collect timei
	collect am
	collect timej
	collect am)))

;;; Return a list with time as a function of distance

(defun time-dist (path)
  (if (not (path-xr path))
      (render-path path))
  (let* ((xcoords (get-x-coords path))
	 (ycoords (get-y-coords path))
	 (tcoords (get-t-coords path)))
    (loop 
	for timei in tcoords
	for timej in (cdr tcoords)
	for xi in xcoords
	for yi in ycoords
	for xj in (cdr xcoords)
	for yj in (cdr ycoords)
	sum (distance (- xj xi)(- yj yi)) into dist
	collect dist 
	collect timej)))

;;; Return a list with distance as a function of time

(defun dist-time (path)
  (if (not (path-xr path))
      (render-path path))
  (let* ((xcoords (get-x-coords path))
	 (ycoords (get-y-coords path))
	 (tcoords (get-t-coords path)))
    (loop 
	for timei in tcoords
	for timej in (cdr tcoords)
	for xi in xcoords
	for yi in ycoords
	for xj in (cdr xcoords)
	for yj in (cdr ycoords)
	sum (distance (- xj xi)(- yj yi)) into dist
	collect timej 
	collect dist)))

;;; Path plotting functions

(defvar plot-stream nil)
(defvar plot-pid nil)

;;; Open a connection to a gnuplot process

(defun open-plot ()
  (multiple-value-bind (input error pid)
      (excl:run-shell-command "/usr/bin/gnuplot" 
			      :wait nil
			      :input :stream
			      :error-output :output)
    (setf plot-stream input
	  plot-pid pid)))

;;; Close and terminate gnuplot

(defun close-plot ()
  (when plot-stream
    (format plot-stream "quit~%")
    (close plot-stream)
    (setf plot-stream nil)
    (multiple-value-bind (status pid)
	(sys:os-wait))))

;;; Plotter routines using the gnuplot process

;;; Set autoscale for selected axes

(defun plot-set-autoscale ()
  (if (not plot-stream)
      (open-plot))
  (format plot-stream "set autoscale~%")
  (finish-output plot-stream))

;;; Set x range

(defun plot-set-x-range (range)
  (if (not plot-stream)
      (open-plot))
  (when range
    (format plot-stream "set xrange [~f:~f]~%"
	    (first range)(second range))
    (finish-output plot-stream)))

;;; Set y range

(defun plot-set-y-range (range)
  (if (not plot-stream)
      (open-plot))
  (when range
    (format plot-stream "set yrange [~f:~f]~%"
	    (first range)(second range))
    (finish-output plot-stream)))

;;; Set z range

(defun plot-set-z-range (range)
  (if (not plot-stream)
      (open-plot))
  (when range
    (format plot-stream "set zrange [~f:~f]~%"
	    (first range)(second range))
    (finish-output plot-stream)))

;;; Set grid

(defun plot-set-grid ()
  (if (not plot-stream)
      (open-plot))
  (format plot-stream "set grid~%")
  (finish-output plot-stream))

;;; Set surface

(defun plot-set-surface ()
  (if (not plot-stream)
      (open-plot))
  (format plot-stream "set surface~%")
  (finish-output plot-stream))

;;; Set parametric mode

(defun plot-set-parametric ()
  (if (not plot-stream)
      (open-plot))
  (format plot-stream "set parametric~%")
  (finish-output plot-stream))

;;; Set ticslevel

(defun plot-set-ticslevel (&optional (level 0))
  (if (not plot-stream)
      (open-plot))
  (format plot-stream "set ticslevel ~s~%" level)
  (finish-output plot-stream))

;;; Plot a supplied curve

(defun plot-curve (curve &key
			 (style "linespoints"))
  (if (not plot-stream)
      (open-plot))
  (let* ((outfile "/zap/plot.dat"))
    (if (probe-file outfile)
	(delete-file outfile))
    (with-open-file (outstream outfile 
		     :direction :output 
		     :if-exists :overwrite
		     :if-does-not-exist :create)
      (loop 
	  for x in curve by #'cddr
	  for y in (cdr curve) by #'cddr
	  do (format outstream "~f ~f~%" x y)))
    (plot-set-grid)
    (format plot-stream "plot ~s" outfile)
    (if style
	(format plot-stream " with ~a~%" style)
      (format plot-stream "~%"))
    (finish-output plot-stream)))

;;; Plot a supplied curve

(defun 3d-plot-curve (3d-curve &key
			       (style "linespoints")
			       (zstyle "impulses")
			       (xrot)
			       (zrot)
			       (scale)
			       (zscale))
  (if (not plot-stream)
      (open-plot))
  (let* ((outfile "/zap/plot.dat")
	 (zmin)
	 (zmax))
    (if (probe-file outfile)
	(delete-file outfile))
    (with-open-file (outstream outfile 
		     :direction :output 
		     :if-exists :overwrite
		     :if-does-not-exist :create)
      (loop 
	  for x in 3d-curve by #'cdddr
	  for y in (cdr 3d-curve) by #'cdddr
	  for z in (cddr 3d-curve) by #'cdddr
	  do (setf zmin (if zmin (if (< z zmin) z) z)
		   zmax (if zmax (if (> z zmax) z) z))
	     (format outstream "~f ~f ~f~%" x y z)))
    (plot-set-grid)
    (plot-set-surface)
    (plot-set-parametric)
    (plot-set-ticslevel 0)
    (if (or xrot zrot scale zscale)
	(format plot-stream "set view ~a,~a,~a,~a~%"
		(if xrot xrot "")
		(if zrot zrot "")
		(if scale scale "")
		(if zscale zscale "")))
    (format plot-stream "splot ~s" outfile)
    (if style
	(format plot-stream " with ~a 1" style))
    (if zstyle
	(format plot-stream ", ~s with ~a 1" outfile zstyle))
    (format plot-stream "~%")
    (finish-output plot-stream)))

;;; Plot the path of a path structure

(defun plot (path)
  (if (not (path-xr path))
      (render-path path))
  (plot-set-autoscale)
  (plot-curve (trajectory path)))

;;; Plot the space-time path of a path structure

(defun 3d-plot (path &key
		     (xrot)
		     (zrot)
		     (scale)
		     (zscale))
  (if (not (path-xr path))
      (render-path path))
  (plot-set-autoscale)
  (3d-plot-curve (3d-trajectory path) 
		 :xrot xrot
		 :zrot zrot
		 :scale scale
		 :zscale zscale))

;;; Plot the velocity of a path

(defun plot-velocity (path)
  (let* ((velocity (velocity path)))
    (plot-set-autoscale)
    
    (plot-set-x-range '(0 1))
    
    (plot-curve velocity :style "steps")))

;;; Plot the doppler shift of a path

(defun plot-doppler (path)
  (if (not (path-xr path))
      (render-path path))
  (plot-set-autoscale)
  (plot-curve (doppler path) :style "steps"))

;;; Plot the acceleration of a path

(defun plot-acceleration (path)
  (if (not (path-xr path))
      (render-path path))
  (plot-set-autoscale)
  (plot-curve (acceleration path) :style "steps"))

;;; Plot time as a function of distance

(defun plot-time (path)
  (if (not (path-xr path))
      (render-path path))
  (plot-set-autoscale)
  (plot-curve (time-dist path) :style "steps"))

;;; Plot distance as a function of time

(defun plot-distance (path)
  (if (not (path-xr path))
      (render-path path))
  (plot-set-autoscale)
  (plot-curve (dist-time path) :style "steps"))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Quad version of NREV (the most popular Samson box reverb)

(defun prime (val)
  (or (= val 2)
      (and (oddp val)
	   (do ((i 3 (+ i 2))
		(lim (sqrt val)))
	       ((or (= 0 (mod val i)) (> i lim))
		(> i lim))))))
       
(defvar dlocnrev-reverb-factor 1.09)
(defvar dlocnrev-lp-coeff 0.7)
(defvar dlocnrev-lp-out-coeff 0.85)
(defvar dlocnrev-output-scale 1.0)
(defvar dlocnrev-volume 1.0)

(definstrument dlocnrev
    (start-time duration 
		&key 
		(reverb-factor dlocnrev-reverb-factor)
		(lp-coeff dlocnrev-lp-coeff)
		(lp-out-coeff dlocnrev-lp-out-coeff)
		(output-scale dlocnrev-output-scale)
		(volume dlocnrev-volume))

  ;; reverb-factor controls the length of the decay -- it should not exceed (/ 1.0 .823)
  ;; lp-coeff controls the strength of the low pass filter inserted in the feedback loop
  ;; output-scale can be used to boost the reverb output

  (let* ((srscale (/ sampling-rate 25641))
	 (val 0)
	 (dly-len (make-array 15 
			      :element-type 'fixnum 
			      :initial-contents '(1433 1601 1867 2053 2251 2399 347 113 37 59 53 43 37 29 19))))
    (loop 
	for i below 15 do
	  (setf val (floor (* srscale (aref dly-len i))))
	  (if (= 0 (mod val 2)) (incf val))
	  (loop while (not (prime val)) do (incf val 2))
	  (setf (aref dly-len i) val))
    (let* ((comb1 (make-comb (* .822 reverb-factor)(aref dly-len 0)))
	   (comb2 (make-comb (* .802 reverb-factor)(aref dly-len 1)))
	   (comb3 (make-comb (* .773 reverb-factor)(aref dly-len 2)))
	   (comb4 (make-comb (* .753 reverb-factor)(aref dly-len 3)))
	   (comb5 (make-comb (* .753 reverb-factor)(aref dly-len 4)))
	   (comb6 (make-comb (* .733 reverb-factor)(aref dly-len 5)))
	   (low (make-one-pole lp-coeff (- lp-coeff 1.0)))
	   (low-a (make-one-pole lp-out-coeff (- lp-coeff 1.0)))
	   (low-b (make-one-pole lp-out-coeff (- lp-coeff 1.0)))
	   (low-c (make-one-pole lp-out-coeff (- lp-coeff 1.0)))
	   (low-d (make-one-pole lp-out-coeff (- lp-coeff 1.0)))
	   (allpass1 (make-all-pass -0.700 0.700 (aref dly-len 6)))
	   (allpass2 (make-all-pass -0.700 0.700 (aref dly-len 7)))
	   (allpass3 (make-all-pass -0.700 0.700 (aref dly-len 8)))
	   (allpass4 (make-all-pass -0.700 0.700 (aref dly-len 10))) ; 9 for stereo 
	   (allpass5 (make-all-pass -0.700 0.700 (aref dly-len 11)))
	   (allpass6 (make-all-pass -0.700 0.700 (aref dly-len 12)))
	   (allpass7 (make-all-pass -0.700 0.700 (aref dly-len 13)))
	   (allpass8 (make-all-pass -0.700 0.700 (aref dly-len 14)))
	   (amp-env (make-env :envelope '(0 1 1 1)
			      :start-time start-time
			      :duration duration))
	   (sample-a 0.0)
	   (sample-b 0.0)
	   (sample-c 0.0)
	   (sample-d 0.0)
	   (rev 0.0)
	   (outrev 0.0)
	   (beg (floor (* start-time sampling-rate)))
	   (end (+ beg (floor (* duration sampling-rate)))))
      (run
       (loop for i from beg to end do
	     (setf rev (* volume (env amp-env)(revin i)))
	     (setf outrev
	       (all-pass allpass4
			 (one-pole low
				   (all-pass allpass3
					     (all-pass allpass2
						       (all-pass allpass1
								 (+ (comb comb1 rev)
								    (comb comb2 rev)
								    (comb comb3 rev)
								    (comb comb4 rev)
								    (comb comb5 rev)
								    (comb comb6 rev))))))))
	     (setf sample-a (* output-scale (one-pole low-a (all-pass allpass5 outrev)))
		   sample-b (* output-scale (one-pole low-b (all-pass allpass6 outrev)))
		   sample-c (* output-scale (one-pole low-c (all-pass allpass7 outrev)))
		   sample-d (* output-scale (one-pole low-d (all-pass allpass8 outrev))))
	     ;; send reverb to output streams, fold back channels if stereo
	     (outa i (if (stereo) 
			 (/ (+ sample-a sample-d) 2)
		       sample-a))
	     (if (or (stereo)(quad))
		 (outb i (if (stereo)
			     (/ (+ sample-b sample-c) 2)
			   sample-b)))
	     (if (quad)
		 (progn (outc i sample-c)
			(outd i sample-d)))))
      (end-run))))