;;; Constraint networks in ordinary programs

(define (harmonic-oscillator net name)
  (let ((n (make-compound-name name)))
    (let ((alpha (cp:make-connector (n 'alpha) net))
	  (2alpha (cp:make-connector (n '2*alpha) net))
	  (omega (cp:make-connector (n 'omega_0) net))
	  (omega2 (cp:make-connector (n 'omega_0^2) net))
	  (Q (cp:make-connector (n 'Q) net))
	  (delta (cp:make-connector (n 'delta) net))
	  (two (cp:make-connector (n 'two) net))
	  (alpha2 (cp:make-connector (n 'alpha^2) net))
	  (root1 (cp:make-connector (n 'root1) net))
	  (root2 (cp:make-connector (n 'root2) net))
	  (radical (cp:make-connector (n 'radical) net)))
      (cp:assume-value two 2)
      (make-multiplier net (n 'm1) 2alpha two alpha)
      (make-multiplier net (n 'm1) omega 2alpha Q)
      (make-squarer net (n 's1) omega2 omega)
      (make-squarer net (n 's2) alpha2 alpha)

      (make-adder net (n 'a1) alpha2 delta omega2)
      (make-adder net (n 'r1) 2alpha root1 root2)
      (make-multiplier net (n 'r2) omega2 root1 root2)
      (make-squarer net (n 'rad) delta radical)
      (make-adder net (n 'r11) root1 radical alpha)
      (make-adder net (n 'r21) alpha root2 radical))))

(define (series-RLC net name)
  (let ((n (make-compound-name name)))
    (let ((R (cp:make-connector (n 'R) net))
	  (L (cp:make-connector (n 'L) net))
	  (C (cp:make-connector (n 'C) net))
	  (LC (cp:make-connector (n 'LC) net))
	  (2alpha (cp:make-connector (n '2*alpha) net))
	  (omega2 (cp:make-connector (n 'omega_0^2) net))
	  (one (cp:make-connector (n 'one) net)))
      (cp:assume-value one 1)
      (make-multiplier net (n 'm1) LC L C)
      (make-multiplier net (n 'm2) one LC omega2)
      (make-multiplier net (n 'm3) R L 2alpha))))

(define (parallel-RLC net name)
  (let ((n (make-compound-name name)))
    (let ((R (cp:make-connector (n 'R) net))
	  (L (cp:make-connector (n 'L) net))
	  (C (cp:make-connector (n 'C) net))
	  (LC (cp:make-connector (n 'LC) net))
	  (RC (cp:make-connector (n 'RC) net))
	  (2alpha (cp:make-connector (n '2*alpha) net))
	  (omega2 (cp:make-connector (n 'omega_0^2) net))
	  (one (cp:make-connector (n 'one) net)))
      (cp:assume-value one 1)
      (make-multiplier net (n 'm1) LC L C)
      (make-multiplier net (n 'm2) one LC omega2)
      (make-multiplier net (n 'm3) RC R C)
      (make-multiplier net (n 'm4) one RC 2alpha))))

(define (spring-mass-dashpot net name)
  (let ((n (make-compound-name name)))
    (let ((m (cp:make-connector (n 'm) net))
	  (b (cp:make-connector (n 'b) net))
	  (k (cp:make-connector (n 'k) net))
	  (2alpha (cp:make-connector (n '2*alpha) net))
	  (omega2 (cp:make-connector (n 'omega_0^2) net)))
      (make-multiplier net (n 'm1) k m omega2)
      (make-multiplier net (n 'm2) b m 2alpha))))

(define (make-harmonic-oscillator-network name)
  (let ((alg (algebraic-solver #f)))
    (let ((net (cp:make-network name (alg 'value-model))))
      (let ((n (make-compound-name name)))
	(harmonic-oscillator net (n 'harmonic))
	(series-RLC net (n 'series))
	(parallel-RLC net (n 'parallel))
	(spring-mass-dashpot net (n 'smd))

	(make-equal net (n 'e1)
	  (connector-named (n '2*alpha:harmonic) net)
	  (connector-named (n '2*alpha:series) net)
	  (connector-named (n '2*alpha:parallel) net)
	  (connector-named (n '2*alpha:smd) net))
	(make-equal net (n 'e2)
          (connector-named (n 'omega_0^2:harmonic) net)
	  (connector-named (n 'omega_0^2:series) net)
	  (connector-named (n 'omega_0^2:parallel) net)
	  (connector-named (n 'omega_0^2:smd) net)))
      net)))

#|
(define net
  (make-harmonic-oscillator-network 'foo))

(for-each (compose pp-comment car) (make-directory net))
;m:smd:foo
;b:smd:foo
;k:smd:foo
;2*alpha:smd:foo
;omega_0^2:smd:foo
;R:parallel:foo
;L:parallel:foo
;C:parallel:foo
;LC:parallel:foo
;RC:parallel:foo
;2*alpha:parallel:foo
;omega_0^2:parallel:foo
;one:parallel:foo
;R:series:foo
;L:series:foo
;C:series:foo
;LC:series:foo
;2*alpha:series:foo
;omega_0^2:series:foo
;one:series:foo
;alpha:harmonic:foo
;2*alpha:harmonic:foo
;omega_0:harmonic:foo
;omega_0^2:harmonic:foo
;Q:harmonic:foo
;delta:harmonic:foo
;two:harmonic:foo
;alpha^2:harmonic:foo
;root1:harmonic:foo
;root2:harmonic:foo
;radical:harmonic:foo
;Unspecified return value

(cp:assume-value (connector-named 'R:series:foo net)
		 1)

(cp:assume-value (connector-named 'L:series:foo net)
		 1/20)

(cp:assume-value (connector-named 'C:series:foo net)
		 1/500)

(cp:value-of (connector-named 'Q:harmonic:foo net))
;Value: 5


(cp:value-of (connector-named 'root1:harmonic:foo net))
;Value: 10+99.498743710662i

(cp:value-of (connector-named 'root2:harmonic:foo net))
;Value: 10-99.498743710662i
|#

#|
(define net
  (make-harmonic-oscillator-network 'foo))

(cp:assume-value (connector-named 'R:series:foo net)
		 'R)

(cp:assume-value (connector-named 'L:series:foo net)
		 'L)

(cp:assume-value (connector-named 'C:series:foo net)
		 'C)

(pe (cp:value-of (connector-named 'Q:harmonic:foo net)))
;(/ (sqrt L) (* R (sqrt C)))
|#

#|
(define net
  (make-harmonic-oscillator-network 'foo))

(cp:assume-value (connector-named 'L:series:foo net)
		 'L)

(cp:assume-value (connector-named 'R:series:foo net)
		 'R)

(cp:assume-value (connector-named 'Q:harmonic:foo net)
		 'Q)

(pe (cp:value-of (connector-named 'C:series:foo net)))
;(/ L (* (expt Q 2) (expt R 2)))
|#

#|
(define net
  (make-harmonic-oscillator-network 'foo))

(cp:assume-value
 (connector-named 'R:series:foo net)
 (& 'R ohm))

(cp:assume-value
 (connector-named 'L:series:foo net)
 (& 'L henry))

(cp:assume-value
 (connector-named 'C:series:foo net)
 (& 'C farad))

(pe
 (cp:value-of
  (connector-named 'Q:harmonic:foo net)))
;(/ (sqrt L) (* R (sqrt C)))


(pe
 (cp:value-of
  (connector-named 'omega_0^2:harmonic:foo net)))
;(& (/ 1 (* C L)) (/ 1 (expt second 2)))

(pe
 (cp:value-of
  (connector-named 'omega_0:harmonic:foo net)))
;(& (/ 1 (* (sqrt C) (sqrt L))) hertz)
|#

#|
(define ((harmonic-system-derivative 2a w2 drive) state)
  (let ((t (ref state 0))
	(x (ref state 1))
	(x. (ref state 2)))
    (let ((x.. (- (drive t) (+ (* 2a x.) (* w2 x)))))
      (vector 1 x. x..))))

(define (test L-series R-series Q drive)
  (let ((net (make-harmonic-oscillator-network 'foo)))
    (let ((f (net->function net
			    '(L:series:foo
			      R:series:foo
			      Q:harmonic:foo)
			    '(2*alpha:harmonic:foo
			      omega_0^2:harmonic:foo))))
      (let ((vs (f L-series R-series Q)))
	(harmonic-system-derivative (car vs)
				    (cadr vs)
				    drive)))))

(pe ((test 'L 'R 'Q (literal-function 'V_s))
     (vector 't 'v_c 'vdot_c)))
(up
 1
 vdot_c
 (/ (+ (* -1 (expt Q 2) (expt R 2) v_c)
       (* (expt L 2) (V_s t))
       (* -1 L R vdot_c))
    (expt L 2)))
|#

(define (net->function net input-names output-names)
  (lambda input-values
    (for-each (lambda (name value)
		(cp:assume-value
		 (connector-named name net)
		 value))
	      input-names
	      input-values)
    (let ((connectors
	   (connectors-named output-names net)))
      (if (not (for-all? connectors cp:has-value?))
	  (error "Outputs of network not set."))
      (map cp:value-of connectors))))


(define (make-directory net)
  (map (lambda (conn)
	 (list (cp:connector-name conn) conn))
       (cp:network-connectors net)))


(define ((make-compound-name outer) inner)
  (symbol inner ': outer))

(define (connectors-named names net)
  (map (lambda (name) (connector-named name net)) names))

(define (connector-named name net)
  (find-matching-item (cp:network-connectors net)
    (lambda (conn)
      (eq? (cp:connector-name conn) name))))

(define (make-adder net name sum a1 a2)
  (cp:make-constraint cp:adder name net
		      sum a1 a2))

(define (make-multiplier net name product m1 m2)
  (cp:make-constraint cp:multiplier name net
		      product m1 m2))


(define (make-negator net name x y)
  (cp:make-constraint cp:negator name net
		      x y))

(define (make-squarer net name y x)
  (cp:make-constraint cp:squarer name net
		      y x))

(define (make-equal net name x . rest)
  (let lp ((l rest) (i 1))
    (if (pair? l)
	(begin
	  (cp:make-constraint cp:B=A
			      (symbol i '% name) net
			      x (car l))
	  (lp (cdr l) (+ i 1))))))