#lang racket

; 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.

(require "unicalc-db.rkt")

;; Export just these functions to other modules
(provide add subtract multiply divide power normalize)

;; 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))))))))

;;;;;;;;;;;;;;;;;;;
;; Normalization ;;
;;;;;;;;;;;;;;;;;;;

;; 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 UQL)
  (make-UQL (get-quant UQL)
        (get-error UQL)
        (cancel (get-den UQL) (get-num UQL))
        (cancel (get-num UQL) (get-den UQL))))


;; function: conv-unit
;; input: a unit, s
;; output: looks up a single unit in the database
;; and returns a QL 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-QL 1.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-QL 1.0 (list s) '()))
        (else (normalize (conv-unit s)))))

;; function: conversion-factor
;; input: a unit, s
;; output: the dimensional-analysis conversion factor for 
;;   converting this unit into normalized units.
;; i.e., given input 'inch, the output would be the quantity
;;   representing "0.0254 meters/inch."
(define (conversion-factor s)
  (let ((UQL (normalize-unit s)))
    (simplify (make-UQL (get-quant UQL)
                        (get-error UQL)
                        (get-num UQL)
                        (cons s (get-den UQL))))))

; function: normalize
; input: a Quantity list, QL
; output: A normalized version of the input QL
; Multiplies the given QL by the conversion factors for the given
;   numerator units, divides by the conversion factors for the given 
;   denominator units, and then cancels out any normalized units
;   appearing in both the numerator and denominator.
(define (normalize QL)
  (let* ((numerator-conversions (map conversion-factor (get-num QL)))
         (denominator-conversions (map conversion-factor (get-den QL))))
    (simplify 
        ;; Unfortunately, divide cares about the order of its arguments.
        ;; and foldr provides list values as first arguments and 
        ;; "running totals" as the second argument. Since we want
        ;; to divide the running total by each list element in turn,
        ;; this means swapping the order of arguments before calling divide.
        (foldr (lambda (x y) (divide y x))
               (foldr multiply QL numerator-conversions)
               denominator-conversions))))

;; 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 (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)
        (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))
    ))