;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; bidirection.lisp ; bidirection calculation using macro ; by T. Shido ; on January 04, 2005 ; modified on January 11, 2005 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (eval-when (:compile-toplevel :load-toplevel :execute) (defun single (obj) (and (car obj) (not (cdr obj)))) (defun group (source n) (if (zerop n) (error "zero length")) (labels ((rec (source acc) (let ((rest (nthcdr n source))) (if (consp rest) (rec rest (cons (subseq source 0 n) acc)) (nreverse (cons source acc)))))) (if source (rec source nil) nil))) (defun find-tree (obj tree) (cond ((null tree) nil) ((atom tree) (eql obj tree)) (t (or (find-tree obj (car tree)) (find-tree obj (cdr tree)))))) (defun position-tree (target lst) (do ((i 0 (1+ i)) (ls1 lst (cdr ls1))) ((not ls1)) (if (find-tree target (car ls1)) (return i)))) (defun remove-nth (n lst) (cond ((not lst) nil) ((zerop n) (cdr lst)) (t (cons (car lst) (remove-nth (1- n) (cdr lst)))))) (defun get-par (form) (labels ((rec (x acc) (cond ((not x) acc) ((atom x) (if (or (eq x 'pi) (numberp x) (fboundp x)) acc (cons x acc))) (t (rec (car x) (rec (cdr x) acc)))))) (rec form nil))) (let ((fun-reverse (make-hash-table :test #'eq))) (dolist (fpair '((exp . log) (sin . asin) (cos . acos) (tan . atan) (sinh . asinh) (cosh . acosh) (tanh . atanh))) (setf (gethash (car fpair) fun-reverse) (cdr fpair) (gethash (cdr fpair) fun-reverse) (car fpair))) (defun solve (target formula) (labels ((rec (this that) (if (eq target this) that (let* ((f (car this)) (term0 (second this)) (terms (cdr this)) (terms-rest (cdr terms)) (p (position-tree target terms))) (case f ((* +) (rec (nth p terms) `(,(if (eq f '*) '/ '-) ,(or that (if (eq f '*) 1.0 0.0)) ,@(remove-nth p terms)))) ((/ -) (if (and (not that) (single terms-rest)) (rec (nth p terms) (car (remove-nth p terms))) (if (eql p 0) (rec term0 `(,(if (eq f '/) '* '+) ,@terms-rest ,(if that that (values)))) (rec (nth p terms) `(,f ,@(remove-nth p terms) ,(if that that (values))))))) (expt (if (eql p 0) (rec term0 (if (eql (car terms-rest) 2) (if that `(sqrt ,that) 1.0) (if that `(expt ,that (/ 1.0 ,(car terms-rest))) 1.0))) (rec (car terms-rest) (if that `(log ,that ,term0) 0.0)))) (sqrt (rec term0 (if that `(expt ,that 2) 1.0))) (log (if (single terms) (rec term0 `(exp ,(or that 1))) (if (eql p 0) (rec term0 (if that `(expt ,(car terms-rest) ,that) (car terms-rest))) (rec (car terms-rest) (if that `(expt ,that (/ 1.0 ,term0)) 1.0))))) ((exp sin cos tan asin acos atan sinh cosh tanh asinh acosh atanh) (rec term0 (list (gethash f fun-reverse) (or that 1.0))))))))) (rec formula nil)))) (defun devide-a-n (ls) (labels ((rec (l0 la ln) (if l0 (let ((o (car l0))) (if (member o ln) (rec (cdr l0) la ln) (if (member o la) (rec (cdr l0) (remove o la) (cons o ln)) (rec (cdr l0) (cons o la) ln)))) (list la ln)))) (rec ls nil nil))) (defun bisec (f min max epsilon) (let ((m (* 0.5 (+ min max)))) (if (< (- max min) epsilon) m (let ((fmin (funcall f min)) (fmax (funcall f max)) (fm (funcall f m))) (cond ((< 0 (* fmin fmax)) (error "wrong range")) ((= 0 fm) m) ((< 0 (* fmin fm)) (bisec f m max epsilon)) ((< 0 (* fmax fm)) (bisec f min m epsilon)) (t nil)))))) (defun round-float (fl c) (declare (long-float fl c)) (multiple-value-bind (r d)(round fl (* 1.0L-1 c)) (declare (ignore r) (long-float d)) (- fl d))) (defun nansw (f v) (let ((p (progn (format t "give r-min, r-max, and epsilon for ~A~%" v) (read-from-string (concatenate 'string "(" (read-line) ")"))))) (round-float (bisec f (first p) (second p) (third p)) (third p)))) (defmacro ave (&rest av) `(/ (+ ,@av) ,(length av))) (defmacro defformula (name form) (let* ((p-in-form0 (get-par form)) (a-n-parms (devide-a-n p-in-form0)) (p-in-form (remove-duplicates p-in-form0))) `(defmacro ,name (&rest av) (if (equal av '(-h)) '',form (let* ((parms (group av 2)) (cparms (mapcar #'car parms)) (unknown (set-difference ',p-in-form cparms))) (unless (single unknown) (error "Invalid unknown parameters.")) `(let ((*WARN-ON-FLOATING-POINT-CONTAGION* nil)) (let ,parms (declare (double-float ,@cparms)) ,(case (car unknown) ,@(mapcar #'(lambda (x) `(,x `(values ,',(solve x form) ',',x))) (first a-n-parms)) ,@(mapcar #'(lambda (x) `(,x `(values (nansw #'(lambda (,',x) ,',form) ',',x) ',',x))) (second a-n-parms)))))))))) ) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; for example ;;; gas equation ; equation of state of ideal gases, PV = nRT (defformula igas (- (* p v) (* n 8.31451 d))) ; equation of state of real gases (van der Waals equation) ; (P + an^2/v^2)(V-nb) = nRT (defformula rgas (- (* (+ p (/ (* a (expt n 2)) (expt v 2))) (- v (* n b))) (* n 8.31451 d))) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Constants of van der Waals equation for several gases ; gases air He H2 O2 N2 CO2 CH4 C2H2 C6H6 ; ak10^-8 Pa m^6 mol^-2l 13.5 0.345 2.485 13.5 14.1 36.5 22.9 45.3 182 ; bk10^-6 m^3 mol^-1l 36.6 23.8 26.7 31.9 39.2 42.8 43.0 57.3 115 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (defmacro air (&rest av) `(rgas a 1.35e-7 b 3.66e-5 ,@av)) (defmacro co2 (&rest av) `(rgas a 3.65e-7 b 4.28e-5 ,@av)) ; Fahrenheit Celsius conversion (defformula fc (- (* 5.0 (- f 32.0)) (* 9.0 c)))