; File: ucErrorProp.rkt
; Date: 9/20/2010
; Author: Christine Alvarado

;; This program implements a unit calculator, with uncertainties in the 
;; measurements. It does simple arithmetic with the units as well as 
;; error propagation.

;; A new datastructure: the UncertainQuantity.  It is a list of three elements 
;; The first element in the list is a list containing the quantity and the error.
;; The second element is a list of the units in the numerator
;; The third element is a list of the units in the denominator
(define (make-QL quant num den)
  (list (list quant 0) num den))

(define (make-UQL quant err num den)
  (list (list quant err) num den))

;; get-quant
;; input: a QL or UQL
;; output: the quantity part of the quantity list
;; Can handle either the old QL representation or the new UQL one
(define get-quant 
  (lambda (x) 
    (if (list? (first x))
        (first (first x))
        (first x))))

;; get-error
;; input: a QL or UQL
;; output: the error part of the quantity list (0 for a QL)
;; Can handle either the old QL representation or the new UQL one
(define get-error 
  (lambda (x) 
    (if (list? (first x))
        (second (first x))
        0)))

; Get the numerator and denominator (respectively) of the 
; input QL or UQL
(define get-num second)
(define get-den third)



;; function: cancel
;; Input: two lists, L and M
;; Output: A list containing the units in M 
;;     that are not in L (i.e., cancels L out of M)
(define (cancel L M)
  (cond ((null? L) M)
        ((null? M) M)
        ((string-ci<? (symbol->string (first L)) 
                      (symbol->string (first M))) (cancel (rest L) M))
        ((string-ci<? (symbol->string (first M))
                      (symbol->string (first L))) (cons (first M) (cancel L (rest M))))
        (else (cancel (rest L) (rest M)))))

;; function: simplify
;; input: A quantity list
;; output: A simplified quantity list where 
;;     units shared between the numerator and denominator
;;     have been eliminated.
(define (simplify QL)
  (make-UQL (get-quant QL) (get-error QL)
        (cancel (get-den QL) (get-num QL))
        (cancel (get-num QL) (get-den QL))))

;; sum-rule
;; input:
;;   dx, dy: two error terms
;; output: the result of applying the error propagation formula
;; for sums to dx and dy
(define (sum-rule dx dy)
  (sqrt (+ (* dx dx) (* dy dy))))

;; prod-rule
;; input:
;;   z: the output of the product
;;   x, y: the inputs
;;   dx, dy: two error terms of the input
;; output: the result of applying the error propagation formula
;;   of the product.
(define (prod-rule z x y dx dy)
  (* z (sqrt (+ (* (/ dx x) (/ dx x)) (* (/ dy y) (/ dy y))))))

;; power-rule
;; input:
;;   z: the output of the product
;;   a: an input value
;;   c: a non-negative integer power
;;   da: the input error
;; output: the result of applying the error propagation formula
;;   for taking the power of an input.
(define (power-rule z a c da)
  (* z (* c (/ da a))))


;; function: multiply
;; input: two uncertain quantity lists.  
;;     If the inputs are normalized 
;;     then the output will be as well
;; output: the UQL QL1*QL2, including error propagation
(define (multiply QL1 QL2)
  (let* (
         (q1 (get-quant QL1))
         (q2 (get-quant QL2))
         (d1 (get-error QL1))
         (d2 (get-error QL2))
         (new-quant (* q1 q2))
         (new-error (prod-rule new-quant q1 q2 d1 d2))
         )
  (simplify (make-UQL new-quant new-error 
                      (sort (append (get-num QL1) (get-num QL2)) 
                            (lambda (s1 s2) (string-ci<? 
                                             (symbol->string s1)
                                             (symbol->string s2))))
                             (sort (append (get-den QL1) (get-den QL2)) 
                                   (lambda (s1 s2) 
                                     (string-ci<? (symbol->string s1)
                                                  (symbol->string s2))))))))

;; function: power
;; input: A uncertain quantity list QL1 and an integer p.
;; output: the UQL QL1^p (with error propagation).  If QL1 is normalized 
;;     then the output will be as well
(define (power QL1 p)
  (let* (
         (q1 (get-quant QL1))
         (d1 (get-error QL1))
         (exact-QL (exact-power (make-QL q1 (get-num QL1) (get-den QL1)) p))
         (new-error (power-rule (get-quant exact-QL) q1 p d1))
         )
    (make-UQL (get-quant exact-QL) new-error (get-num exact-QL) (get-den exact-QL))))
   
;; function: exact-power
;; input: A uncertain quantity list QL1 with 0 error and an integer p.
;; output: the QL QL^p (with error propagation).  
;;  Returns a QL raised to a power, but returns an error if the uncertainty is not 0
(define (exact-power QL p)
  (if (not (= 0 (get-error QL)))
      (list 'error "exact-power QL must have 0 uncertainty" QL)
      (cond
        ((< p 0) (exact-power (divide (make-UQL 1.0 0.0 '() '()) QL) (* -1 p)))
        ((= p 0) (make-UQL 1.0 0.0 '() '()))
        (else (multiply QL (exact-power QL (- p 1)))))))
      
;; function: divide
;; input: two uncertain quantity lists.  If the inputs are normalized 
;;     then the output will be as well
;; output: the QL QL1/QL2
(define (divide QL1 QL2)
  (let* (
         (q1 (get-quant QL1))
         (q2 (get-quant QL2))
         (d1 (get-error QL1))
         (d2 (get-error QL2))
         (new-quant (/ q1 q2))
         (new-error (prod-rule new-quant q1 q2 d1 d2))
         )
  (simplify (make-UQL new-quant new-error 
                      (sort (append (get-num QL1) (get-den QL2)) 
                            (lambda (s1 s2) (string-ci<? (symbol->string s1)
                                                         (symbol->string s2))))
                      (sort (append (get-den QL1) (get-num QL2)) 
                            (lambda (s1 s2) 
                              (string-ci<? (symbol->string s1)
                                           (symbol->string s2))))))))

;; function: conv-unit
;; input: a unit, s
;; output: looks up a single unit in the database
;; and returns a UQL for that unit if it is basic,
;; or the conversion entry in the database if it is not.
(define (conv-unit s)
  (cond ((assoc s unicalc-db) (second (assoc s unicalc-db)))
        (else (make-UQL 1.0 0 (list s) ()))))

;; function: normalize-unit
;; input: a unit, s
;; output: A normalized quantity list equivalent to the 
;;   input unit.
(define (normalize-unit s)
  (cond ((not (assoc s unicalc-db)) (make-UQL 1.0 0 (list s) ()))
        (else (normalize (conv-unit s)))))

;; function: normalize
;; input: an uncertain Quantity list
;; output: A normalized version of the input QL
(define (normalize QL)
  (cond ((null? (append (get-num QL) (get-den QL))) QL)
        ((null? (get-num QL)) (divide (make-UQL (get-quant QL) 
                                                (get-error QL) '() '())
                                      (normalize (make-QL 1.0 (get-den QL) '()))))
        (else (multiply (normalize-unit (first (get-num QL)))
                        (normalize (make-UQL (get-quant QL) (get-error QL)
                                    (rest (get-num QL))
                                    (get-den QL)))))))

;; function: add
;; input: Two uncertain quantity lists.  This version does not 
;;   assume the inputs are normalized.
;; output: The UQL resulting from adding the two input QLs (with error prop),
;;   or an error if the QLs are not in terms of the same units.
;;   The output is always normalized.
(define (add Q1 Q2)
  (let ((norm-Q1 (normalize Q1))
        (norm-Q2 (normalize Q2)))
    (cond 
      ((units-match norm-Q1 norm-Q2) 
       (let* ((q1 (get-quant norm-Q1))
              (q2 (get-quant norm-Q2))
              (d1 (get-error norm-Q1))
              (d2 (get-error norm-Q2))
              (new-quant (+ q1 q2))
              (new-error (sum-rule d1 d2)))
         (make-UQL new-quant new-error (get-num norm-Q1) (get-den norm-Q1))))
      (else (list 'error "illegal add" Q1 Q2))))) 


;; function: subtract
;; input: Two uncertain quantity lists, QL1 and QL2  
;;   This version does not 
;;   assume the inputs are normalized.
;; output: The UQL resulting from QL1-QL2, with error prop,
;;   or an error if the QLs are not in terms of the same units.
;;   The output is always normalized.
(define (subtract QL1 QL2)
  (let* ((negQL2 (make-UQL (* -1 (get-quant QL2)) (get-error QL2) 
                                 (get-num QL2) (get-den QL2)))
         (sumQL (add QL1 negQL2)))
    (if (equal? (first sumQL) 'error)
        (list 'error "illegal subtract" QL1 QL2)
        sumQL)))

;; function: units-match
;; input: two UQLs, QL1 and QL2
;; output: #t if the units in QL1 are identical to the units
;;   in QL2, #f otherwise.
(define (units-match QL1 QL2)
  (let (
        (sortedQ1num (sort (get-num QL1)
                           (lambda (s1 s2) (string-ci<? (symbol->string s1)
                                                        (symbol->string s2)))))
        (sortedQ2num (sort (get-num QL2)
                           (lambda (s1 s2) (string-ci<? (symbol->string s1)
                                                        (symbol->string s2)))))
        (sortedQ1den (sort (get-den QL1)
                           (lambda (s1 s2) (string-ci<? (symbol->string s1)
                                                        (symbol->string s2)))))
        (sortedQ2den (sort (get-den QL2)
                           (lambda (s1 s2) (string-ci<? (symbol->string s1)
                                                        (symbol->string s2)))))
        )
    (and (equal? sortedQ1num sortedQ2num) (equal? sortedQ1den sortedQ2den))
    ))
                     
         
(load "unicalc-db.rkt")
(load "tester.rkt") 
(load "unicalc-tests.rkt")
(load "unicalc-with-error-tests.rkt")