;;; strad.scm -- Translation CLM -> Snd

;; Bowed string physical model with stiffness.  CLM version adapted
;; from the Matlab and C versions courtesy of JOS and Stefania Serafin
;; from code revised on 7/14/01

;; CLM version by Juan Reyes
;; SND version by Michael Scholz (based on strad.ins)
;; revised by Bill to suit the run macro

(provide 'snd-strad.scm)
(if (not (provided? 'snd-ws.scm)) (load "ws.scm"))
;;; (if (not (provided? 'snd-jcrev.scm)) (load "jcrev.scm"))


(definstrument (bow beg dur frq amplitude
		    (bufsize 2205)
		    (fb 0.2) ;; bow force: between 0.0 and 1.0
		    (vb 0.05)	;; bow velocity: between 0.0 and 0.8
		    (bp 0.08)	;; bow position: 0.0=bridge; 0.5=middle of string; 1.0=Nut 
		    (inharm 0.1) ;; inharmonicity: 0.0 harmonic; 1.0 not harmonico
		    (ampenv '(0 1 15 1 95 1 100 0))
		    (degree 45) (dist 0.0025) (reverb 0))
  (let* ((beg (seconds->samples beg))
	 (len (seconds->samples dur))
	 (end (+ beg len))
	 (ampf (make-env :envelope ampenv :scaler amplitude :duration dur))
	 (loc (make-locsig degree dist reverb))
	 (vinut      (make-vct bufsize)) 
	 (vinbridge  (make-vct bufsize))
	 (vinutt     (make-vct bufsize))
	 (vinbridget (make-vct bufsize))
	 (vib 0.0) (vin 0.0) (vibt 0.0) (vint 0.0)
	 (freq frq)
	 (mus 0.8)
	 (twavespeedfactor 5.2)
	 (posl 0) (posr 0) 
	 (poslt 0) (posrt 0)
	 (indexl 0) (indexr 0)
	 (indexlt 0) (indexrt 0)
	 (indexl_1 0) (indexr_1 0)
	 (indexlt_1 0) (indexrt_1 0)
	 (indexl_2 0) (indexr_2 0)
	 (indexlt_2 0) (indexrt_2 0)
	 (updl 0) (updr 0)
	 (updlt 0) (updrt 0)
	 (b0b 0.859210)
	 (b1b -0.704922)
	 (b2b 0.022502)
	 (a1b -0.943639)
	 (a2b  0.120665)
	 (b0n  7.0580050e-001)
	 (b1n -5.3168461e-001)
	 (b2n  1.4579750e-002)
	 (a1n  -9.9142489e-001)
	 (a2n   1.8012052e-001)
	 (b0bt 9.9157155e-001)
	 (b1bt -8.2342890e-001)
	 (b2bt  8.8441749e-002)
	 (a1bt -8.3628218e-001)
	 (a2bt  9.2866585e-002)
	 (b0nt 4.3721359e-001)
	 (b1nt -2.7034968e-001)
	 (b2nt -5.7147560e-002)
	 (a1nt -1.2158343e+000)
	 (a2nt  3.2555068e-001)
	 (xm1bt 0.0d0) (xm2bt 0.0d0)
	 (xm1nt 0.0d0)(xm2nt 0.0d0)
	 (ym1bt 0.0d0) (ym2bt 0.0d0) 
	 (ym1nt 0.0d0) (ym2nt 0.0d0)
	 (xm1b 0.0d0) (xm2b 0.0d0)(ym1b 0.0d0) (ym2b 0.0d0)
	 (xm1n 0.0d0) (xm2n 0.0d0)(ym1n 0.0d0) (ym2n 0.0d0)
	 (ynb 0.0d0) (ynbt 0.0d0)
	 (ynn 0.0d0) (ynnt 0.0d0)
	 (ya1nb 0.0d0) (ynba1 0.0d0)
	 (y1nb 0.0d0) 
	 (vh 0.0d0)
	 (aa 0.0d0) (bb1 0.0d0) (cc1 0.0d0) (delta1 0.0d0)
	 (bb2 0.0d0) (cc2 0.0d0) (delta2 0.0d0)
	 (v 0.0d0) (v1 0.0d0) (v2 0.0d0)
	 (rhs #f) (lhs #f)
	 (vtemp 0.0d0)
	 (f 0.0d0)
	 (stringImpedance 0.55)
	 (stringImpedancet 1.8)
	 (stick 0)
	 (zslope (/ 1 (+ (/ 1 (* 2 stringImpedance)) (/ 1 (* 2 stringImpedancet)))))
	 (xnn 0.0d0) (xnb 0.0d0)
	 (xnnt 0.0d0) (xnbt 0.0d0)
	 (alphar 0) (alphal 0)
	 (alphart 0) (alphalt 0)
	 (len (- (/ (mus-srate) freq ) 2))
	 (lent (/ (- (/ (mus-srate) freq) 2) twavespeedfactor))
	 (del_right (* len bp))
	 (del_left (* len (- 1 bp)))
	 (del_leftt (* lent (- 1 bp)))
	 (del_rightt (* lent bp))
	 (samp_rperiod (floor del_right))
	 (samp_lperiod  (floor del_left))
	 (samp_lperiodt (floor del_leftt))
	 (samp_rperiodt (floor del_rightt)))

    (define (bowfilt inharmon)
      (set! ynb (- (- (+ (* b0b vib) (* b1b xm1b) (* b2b xm2b)) (* a1b ym1b)) (* a2b ym2b)))
      (set! xm2b  xm1b)
      (set! xm1b  vib)
      (set! ym2b  ym1b)
      (set! ym1b  ynb)
      (set! ynn (- (- (+ (* b0n vin) (* b1n xm1n) (* b2n xm2n)) (* a1n ym1n)) (* a2n ym2n)))
      (set! xm2n  xm1n)
      (set! xm1n  vin)
      (set! ym2n  ym1n)
      (set! ym1n  ynn)
      (set! ynbt (- (- (+ (* b0bt vibt) (* b1bt xm1bt) (* b2bt xm2bt))
		       (* a1bt ym1bt)) (* a2bt ym2bt)))
      (set! xm2bt xm1bt)
      (set! xm1bt vibt)
      (set! ym2bt ym1bt)
      (set! ym1bt ynbt)
      (set! ynnt (- (- (+ (* b0nt vint) (* b1nt xm1nt) (* b2nt xm2nt))
		       (* a1nt ym1nt)) (* a2nt ym2nt)))
      (set! xm2nt  xm1nt)
      (set! xm1nt  vint)
      (set! ym2nt  ym1nt)
      (set! ym1nt  ynnt)
      (if (<= inharmon 0.00001) (set! inharmon 0.00001))
      (if (>= inharmon 0.9999) (set! inharmon 0.9999))
      (set! y1nb (+ (* -1 (* inharmon ynb)) ynba1 (* inharmon ya1nb)))
      (set! ya1nb y1nb)
      (set! ynba1 ynb)
      (set! y1nb (* -1 y1nb))
      (set! ynn (* -1 ynn))
      (set! ynbt (* -1 ynbt)))
    
    (if (< samp_rperiod 0) (set! samp_rperiod 0))
    (if (> samp_rperiod (- bufsize 1)) (set! samp_rperiod (- bufsize 1)))
    (if (< samp_lperiod 0) (set! samp_lperiod 0))
    (if (> samp_lperiod (- bufsize 1)) (set! samp_lperiod (- bufsize 1)))
    (set! alphar (exact->inexact (- del_right samp_rperiod)))
    (set! alphal (exact->inexact (- del_left samp_lperiod)))
    (if (< samp_rperiodt 0)(set! samp_rperiodt 0))
    (if (> samp_rperiodt (- bufsize 1)) (set! samp_rperiodt (- bufsize 1)))
    (if (< samp_lperiodt 0) (set! samp_lperiodt 0))
    (if (> samp_lperiodt (- bufsize 1)) (set! samp_lperiodt (- bufsize 1)))
    (set! alphart (exact->inexact (- del_rightt samp_rperiodt)))
    (set! alphalt (exact->inexact (- del_leftt samp_lperiodt)))
    (set! posr (modulo (inexact->exact (+ end posr)) bufsize))
    (set! posl (modulo (inexact->exact (+ end posl)) bufsize))
    (set! posrt (modulo (inexact->exact (+ end posrt)) bufsize))
    (set! poslt (modulo (inexact->exact (+ end poslt)) bufsize))
    (ws-interrupt?)
    (run
     (do ((i beg (+ i 1)))
	 ((= i end))
       (set! indexl (modulo (inexact->exact (- (+ i posl bufsize) samp_lperiod)) bufsize))
       (set! indexr (modulo (inexact->exact (- (+ i posr bufsize) samp_rperiod)) bufsize))
       (set! indexlt (modulo (inexact->exact (- (+ i poslt bufsize) samp_lperiodt)) bufsize))
       (set! indexrt (modulo (inexact->exact (- (+ i posrt bufsize) samp_rperiodt)) bufsize))
       (set! indexl_1 (modulo
		       (inexact->exact (- (- (+ i posl bufsize) samp_lperiod) 1)) bufsize))
       (set! indexr_1 (modulo
		       (inexact->exact (- (- (+ i posr bufsize) samp_rperiod) 1)) bufsize))
       (set! indexlt_1 (modulo
			(inexact->exact (- (- (+ i poslt bufsize) samp_lperiodt) 1)) bufsize))
       (set! indexrt_1 (modulo
			(inexact->exact (- (- (+ i posrt bufsize) samp_rperiodt) 1)) bufsize))
       (set! indexl_2 (modulo
		       (inexact->exact (- (- (+ i posl bufsize) samp_lperiod) 2)) bufsize))
       (set! indexr_2 (modulo
		       (inexact->exact (- (- (+ i posr bufsize) samp_rperiod) 2)) bufsize))
       (set! indexlt_2 (modulo
			(inexact->exact (- (- (+ i poslt bufsize) samp_lperiodt) 2)) bufsize))
       (set! indexrt_2 (modulo
			(inexact->exact (- (- (+ i posrt bufsize) samp_rperiodt) 2)) bufsize))
       (set! vib (+ (/ (* (vct-ref vinbridge indexl_2) (- alphal 1)(- alphal 2)) 2)
		    (* (vct-ref vinbridge indexl_1) (* alphal -1) (- alphal 2))
		    (/ (* (vct-ref vinbridge indexl) alphal (- alphal 1)) 2)))
       (set! vin (+ (/ (* (vct-ref vinut indexr_2) (- alphar 1)(- alphar 2)) 2)
		    (* (vct-ref vinut indexr_1) (* alphar -1) (- alphar 2))
		    (/ (* (vct-ref vinut indexr) (- alphar 1) alphar) 2)))
       (set! vibt (+ (/ (* (vct-ref vinbridget indexlt_2) (- alphalt 1)(- alphalt 2)) 2)
		     (* (vct-ref vinbridget indexlt_1) (* alphalt -1) (- alphalt 2))
		     (/ (* (vct-ref vinbridget indexlt) alphalt (- alphalt 1)) 2)))
       (set! vint (+ (/ (* (vct-ref vinutt indexrt_2) (- alphart 1)(- alphart 2)) 2)
		     (* (vct-ref vinutt indexrt_1) (* alphart -1) (- alphart 2))
		     (/ (* (vct-ref vinutt indexrt) (- alphart 1) alphart) 2)))
       (bowfilt inharm)
       (set! vh (+ (+ ynn y1nb) (+ ynnt ynbt)))
       
       (set! aa zslope)
       (set! bb1 (- (- (+ (* 0.2 zslope) (* 0.3 fb)) (* zslope vb)) (* zslope vh)))
       (set! cc1 (- (- (+ (* 0.06 fb) (* (* zslope vh) vb)) (* 0.2 zslope vh)) (* 0.3 vb fb)))
       (set! delta1 (- (* bb1 bb1) (* 4 aa cc1)))
       (set! bb2 (- (- (- (* -0.2 zslope) (* 0.3 fb)) (* zslope vb)) (* zslope vh)))
       (set! cc2 (+ (+ (+ (+ (* 0.06 fb) (* zslope vh vb))
			  (* 0.2 zslope vh)) (* 0.3 vb fb)) (* 0.1 fb)))
       (set! delta2 (- (* bb2 bb2) (* 4 aa cc2)))
       (if (or (= vb 0) (= fb 0))
	   (set! v vh)
	   (begin
	     (if (= vh vb)
		 (begin
		   (set! v vb)
		   (set! stick 1))
		 (begin
		   (if (> vh vb)
		       (begin
			 (set! lhs #f)
			 (set! rhs #t))
		       (begin
			 (set! rhs #f)
			 (set! lhs #t)))
		   (if rhs
		       (begin
			 (if (< delta1 0)
			     (begin
			       (set! v vb)
			       (set! stick 1))
			     (begin
			       (if (= stick 1)
				   (begin
				     (set! vtemp vb)
				     (set! f (* (* 2 zslope) (- vtemp vh)))
				     (if (>= f (* -1 (* mus fb)))
					 (set! v vtemp)
					 (begin
					   (set! v1 (/ (+ (* -1 bb1 ) (sqrt delta1)) (* 2 aa)))
					   (set! v2 (/ (- (* -1 bb1) (sqrt delta1)) (* 2 aa)))
					   (set! v (min v1 v2))
					   (set! stick 0))))
				   (begin
				     (set! v1 (/ (+ (* -1 bb1 ) (sqrt delta1)) (* 2 aa)))
				     (set! v2 (/ (- (* -1 bb1) (sqrt delta1)) (* 2 aa)))
				     (set! v (min v1 v2))
				     (set! stick 0))))))
		       (if lhs
			   (begin
			     (if (< delta2 0)
				 (begin
				   (set! v vb)
				   (set! stick 1))
				 (begin
				   (if (= stick 1)
				       (begin
					 (set! vtemp vb)
					 (set! f (* zslope (- vtemp vh)))
					 (if (and (<= f (* mus fb)) (> f 0))
					     (begin
					       (set! v vtemp))
					     (begin
					       (set! v1 (/ (- (* -1 bb2 ) (sqrt delta2)) (* 2 aa)))
					       (set! v2 (/ (+ (* -1 bb2) (sqrt delta2)) (* 2 aa)))
					       (set! vtemp (min v1 v2))
					       (set! stick 0)
					       (if (> vtemp vb)
						   (begin
						     (set! v vb)
						     (set! stick 1))
						   (begin
						     (set! v vtemp)
						     (set! f (* zslope (- v vh) )))))))
				       (begin
					 (set! v1 (/ (- (* -1 bb2 ) (sqrt delta2)) (* 2 aa)))
					 (set! v2 (/ (+ (* -1 bb2) (sqrt delta2)) (* 2 aa)))
					 (set! v (min v1 v2))
					 (set! stick 0)))))
			     (if (> v vb)
				 (begin
				   (set! v vb)
				   (set! stick 1))))))))
	     (set! f (* zslope (- v vh)))
	     (set! xnn (+ y1nb (/ f (* 2 stringImpedance))))
	     (set! xnb (+ ynn (/ f (* 2 stringImpedance))))))
       
       (set! f (* zslope (- v vh)))
       (set! xnnt (+ ynbt (/ f (* 2 stringImpedancet))))
       (set! xnbt (+ ynnt (/ f (* 2 stringImpedancet))))
       (set! updl (modulo (inexact->exact (+ i posl bufsize)) bufsize))
       (set! updr (modulo (inexact->exact (+ i posr bufsize)) bufsize))
       (set! updlt (modulo (inexact->exact (+ i poslt bufsize)) bufsize))
       (set! updrt (modulo (inexact->exact (+ i posrt bufsize)) bufsize))
       (vct-set! vinbridge updl xnb)
       (vct-set! vinut updr xnn)
       (vct-set! vinbridget updlt xnbt)
       (vct-set! vinutt updrt xnnt)
       (locsig loc i (* xnb (env ampf)))
       (set! lhs #f)
       (set! rhs #f)))))

					;(with-sound (:channels 2) (bow 0 3 400 0.5 :vb 0.15 :fb 0.1 :inharm 0.25))
					;(with-sound (:channels 2) (bow 0 2 440 0.5  :fb 0.25))
					;(with-sound (:channels 2) (bow 0 4 600 0.8))
					;(with-sound (:channels 2) (bow 0 6 147 2 :fb 0.035 :vb 0.1))
					;(with-sound (:channels 2) (bow 0 3 1100 0.5 :vb 0.45 :fb 0.9 :inharm 0.3))
					;(with-sound (:channels 2) (bow 0 3 1500 0.5 :vb 0.25 :fb 0.9 :inharm 0.3))
					;(with-sound (:channels 2) (bow 0 3 1525 0.5 :vb 0.25 :fb 0.9 :inharm 0.3))
					;(with-sound (:channels 2 :reverb jc-reverb) (bow 0 1 400 0.5 :reverb 0.0051))
					;
					;(with-sound (:channels 2 :reverb jc-reverb)
					; 	    (bow 0 3 366 0.5  :degree 0)
					; 	    (bow 0 3 422 0.5  :degree 90)
					; 	    (bow 4 6 147 2 :fb 0.035 :vb 0.1 :reverb 0.051))

;; strad.scm ends here
