Commit 043a16fd authored by Liam M. Healy's avatar Liam M. Healy
Browse files

Physical quantities with grids

Start implementation of physical quantities with grid magnitudes.
Compiles and loads cleanly but mostly non-functional or only
accidentally functional.
parent 832bb03f
;; Antik system definition
;; Liam Healy 2010-12-24 09:43:28EST antik.asd
;; Time-stamp: <2011-01-01 12:49:52EST antik.asd>
;; Time-stamp: <2011-01-02 18:52:11EST antik.asd>
;; Copyright 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -96,7 +96,8 @@
(:file "unit-definitions" :depends-on ("units"))
(:file "sysunit-definitions"
:depends-on ("units" "unit-definitions"))
(:file "physical-quantities" :depends-on ("units"))
(:file "class")
(:file "physical-quantities" :depends-on ("units" "class"))
(:file "undimension" :depends-on ("physical-quantities"))
(:file "arithmetic" :depends-on ("units" "undimension"))
(:file "iso8601")
......
;; Rotations in 2D and 3D
;; Liam Healy Mon Sep 4 2000 - 17:05
;; Time-stamp: <2010-12-26 23:02:36EST rotation.lisp>
;; Time-stamp: <2011-01-01 12:05:05EST rotation.lisp>
;; Copyright 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -82,7 +82,6 @@
(let ((ang (if passive angle (- angle))))
(+ (grid:identity-matrix 3 (cos ang))
(* (- 1.0 (cos ang))
;; There's got to be a rectangle way of doing this.
;; Dot product.
(grid:matrix-from-columns
(* rotation-axis (grid:gref rotation-axis 0))
......
;; Methods for grid functions
;; Liam Healy 2009-12-21 11:19:00EST methods.lisp
;; Time-stamp: <2010-11-22 10:41:05EST methods.lisp>
;; Time-stamp: <2011-01-02 14:15:18EST methods.lisp>
;;
;; Copyright 2009, 2010 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -443,3 +443,11 @@ error while parsing arguments to DEFINE-COMPILER-MACRO (SETF GREF*):
(cons 'progn
(loop for (thing value) on args by #'cddr
collect (my-setf-pair thing value))))
;;;;****************************************************************************
;;;; Simple math functions
;;;;****************************************************************************
(defmethod antik:zerop ((grid foreign-array))
(loop for ind from 0 below (total-size grid)
always (antik:zerop (gref* grid ind))))
;; Functions on grids
;; Liam Healy 2010-06-19 22:25:17EDT functions.lisp
;; Time-stamp: <2010-12-26 17:15:19EST functions.lisp>
;; Time-stamp: <2011-01-02 14:13:28EST functions.lisp>
;; Copyright 2010 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -130,6 +130,14 @@
"Make the object(s) required for storage of the data,
if it consists of CL arrays and lists."))
;;;;****************************************************************************
;;;; Simple math functions
;;;;****************************************************************************
(defmethod antik:zerop ((grid array))
(loop for ind from 0 below (total-size grid)
always (antik:zerop (gref* grid ind))))
;;;;****************************************************************************
;;;; Element reference
;;;;****************************************************************************
......
;; Generic functions for generalized numbers
;; Liam Healy Tue Feb 9 1999 - 21:38
;; Time-stamp: <2010-12-25 12:33:35EST generic.lisp>
;; Time-stamp: <2011-01-02 18:49:57EST generic.lisp>
;; Copyright 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -382,15 +382,15 @@
(every #'zerop x))
(:method (a) (cl:zerop (numbcomp a))))
(defgeneric floor (a b)
(:documentation "Greatest multiple of b less than a.")
(:method ((a number) (b number))
(cl:floor a b)))
(defgeneric floor (number &optional divisor)
(:documentation "Greatest multiple of divisor less than number.")
(:method ((number number) &optional (divisor 1))
(cl:floor number divisor)))
(defgeneric round (a b)
(:documentation "Nearest multiple of b to a.")
(:method ((a number) (b number))
(cl:round a b)))
(defgeneric round (number &optional divisor)
(:documentation "Nearest multiple of divisor to number.")
(:method ((number number) &optional (divisor 1))
(cl:round number divisor)))
(defgeneric signum (a)
(:documentation "Signum(a).")
......
;; Floating point numbers.
;; Liam Healy Wed May 8 2002 - 10:47
;; Time-stamp: <2010-12-24 16:51:31EST float.lisp>
;; $Id$
;; Time-stamp: <2011-01-02 18:31:44EST float.lisp>
;; Copyright 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -107,18 +106,18 @@
'(#\ )
(format nil
(if print-sign "~v,v@g" "~v,vg")
(+ significant-figures
(if (or print-sign (minusp float)) 7 6))
(+ significant-figures
;; compensate for misdesign of ~g and subtract
;; an extra digit of precision in exponential form,
;; so the same number of significant digits are printed in fixed
;; and exponential forms.
(if (or (zerop float)
(<= 0
(1+ (floor (log (abs float) 10)))
significant-figures))
0 -1))
(cl:+ significant-figures
(if (or print-sign (cl:minusp float)) 7 6))
(cl:+ significant-figures
;; compensate for misdesign of ~g and subtract
;; an extra digit of precision in exponential form,
;; so the same number of significant digits are printed in fixed
;; and exponential forms.
(if (or (zerop float)
(cl:<= 0
(1+ (cl:floor (cl:log (cl:abs float) 10)))
significant-figures))
0 -1))
float))))
;;; (loop for exp from -8 upto 8 do (print (format-float-switch (* pi (expt 10 exp)) 6)))
......
;; Format numerical output
;; Liam Healy Wed Dec 11 2002 - 16:37
;; Time-stamp: <2010-12-24 17:15:38EST format-output.lisp>
;; Time-stamp: <2011-01-02 18:20:33EST format-output.lisp>
;; Copyright 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -70,7 +70,7 @@
`(if ,stream (princ ,object ,stream) ,object))
(defgeneric nf (object &optional stream)
(:documentation "Format output for numerica objects.
(:documentation "Format output for numerical objects.
If stream is nil, return the formatted output as a string.")
(:method ((object t) &optional (stream *standard-output*))
;; by default, just print to stream
......
;; Arithmetic with physical quantities.
;; Liam Healy Fri Apr 22 2005 - 09:58
;; Time-stamp: <2010-12-25 16:08:36EST arithmetic.lisp>
;; Time-stamp: <2011-01-02 18:55:35EST arithmetic.lisp>
;; Copyright 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -202,7 +202,7 @@
(undimensioned-should-not-be-here)
(antik:plusp (pqval x)))
(defmethod floor ((x physical-quantity) y)
(defmethod floor ((x physical-quantity) &optional (y 1))
(undimensioned-should-not-be-here)
(antik:floor (pqval x) (pqval y)))
......@@ -217,15 +217,3 @@
(equal-dimension x y)
(-nz (pqval x) (pqval y) relative-difference))
;;; Special handling of arrays with pqs
(defun *pq
(matrix vector &optional return-units)
"Multiple the matrix by the vector, taking care to remove the units
from the matrix and vector and add units on the result vector,
by default the same units the original vector had."
(undimensioned-should-not-be-here)
(multiple-value-bind (value units allsame)
(pqval vector)
(if allsame
(make-pq (antik:* (pqval matrix) value) (or return-units units))
(make-pq-seq (antik:* (pqval matrix) value) (or return-units units)))))
;; Physical quantity with a magnitude that is a grid.
;; Liam Healy 2011-01-01 12:05:37EST grid.lisp
;; Time-stamp: <2011-01-02 19:08:19EST class.lisp>
;; Copyright 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program. If not, see <http://www.gnu.org/licenses/>.
;;; Physical quantities with a magnitude that is a grid (non-scalar).
(in-package :antik)
(export '(make-pq *zero-is-dimensionless*))
;;;;****************************************************************************
;;;; The class and functions to make instances
;;;;****************************************************************************
(defclass physical-quantity ()
((magnitude :initarg :magnitude :reader pq-magnitude)
(dimension :initarg :dimension :reader pq-dimension :type list)
(scalar-units :initarg :scalar-units :reader scalar-units))
(:documentation "A quantity with a dimension and possibly a unit, e.g. length."))
(defparameter *zero-is-dimensionless* t
"Numbers with zero magnitude are made dimensionless.")
;;; (make-pq 33.4 "meters/second")
;;; (make-pq 88 '(/ feet minute))
;;; (make-pq 5 'speed)
;;; #_5.0_KM/S
;;; Map into structures?
(defun make-pq-OBSOLETE
(magnitude unit-expression &optional (sysunits *system-of-units*) seq)
"Make the physical quantity. "
(if (and seq (typep magnitude 'sequence))
(make-pq-seq magnitude unit-expression sysunits)
(make-pq-int magnitude unit-expression sysunits)))
(defun make-pq-zero (unit-expression &optional (sysunits *system-of-units*))
(let ((*zero-is-dimensionless* nil)) (make-pq 0 unit-expression sysunits)))
(defgeneric make-pq
(magnitude unit-expression &optional sysunits scalar-units)
(:documentation "Make a physical quantity. Allow for the possiblity
that magnitude is already a pq; in this case, check that the
physical dimension is correct, and pass it on. In this case, this
function acts as a way of assigning default units to a number. If
sysunits is specified, unitless dimensions (e.g., 'length) are
interpeted as units in that system. If scalar-units is T and
magnitude is not a scalar, the unit-expression applies to each
element of the magnitude. Otherwise it must be a grid of the shape
as the magnitude."))
(defun make-pq-dimel (mag dimel &optional (scalar-units t))
"Make a physical quantity directly from dimel, unless dimensionless, in which
case the plain number is returned."
(if (or (equal dimel (dimension 'dimensionless))
(and *zero-is-dimensionless* (zerop mag)))
mag
(nth-value
0
(pqwu
(make-instance
'physical-quantity
:magnitude mag :dimension dimel :scalar-units scalar-units)))))
(defun make-pq-ue (magnitude unit-expression sysunits scalar-units)
(if (not unit-expression)
magnitude
(multiple-value-bind (el convfact)
(parse-units unit-expression sysunits)
(make-pq-dimel
(if convfact
(antik:* convfact magnitude)
(error
"~a is not a unit of measure but a physical dimension ~
and cannot be interpreted as a unit of measure ~
without a system of units."
unit-expression))
el
scalar-units))))
(defmethod make-pq
((magnitude number) unit-expression
&optional (sysunits *system-of-units*) (scalar-units t))
(declare (ignore scalar-units)) ; they're irrelevant for scalars
(when (complexp magnitude)
(cerror "Accept" "Magnitude ~a is complex." magnitude))
(make-pq-ue magnitude unit-expression sysunits t))
(defmethod make-pq
((pq physical-quantity) unit-expression
&optional (sysunits *system-of-units*) (scalar-units t))
(declare (ignore sysunits scalar-units))
(check-dimension pq unit-expression)
pq)
(defmethod make-pq
(magnitude unit-expression
&optional (sysunits *system-of-units*) (scalar-units t))
;; This method is for making a pq from a grid; can't be captured
;; otherwise in a single method because grids don't have a single
;; class.
(if (grid:gridp magnitude)
(if scalar-units
(make-pq-ue magnitude unit-expression sysunits t)
(make-pq-ue
magnitude
(make-array (grid:dimensions magnitude)
:initial-contents
unit-expression)
sysunits nil))
(error "Don't know how to make physical quantity.")))
(defmethod grid::make-grid-data :around
(type dimensions rest-spec
&key (initial-element nil initial-element-p)
(initial-contents nil initial-contents-p))
;; Check for pqs in the :initial-element and/or :initial-contents.
(cond (initial-contents-p
(multiple-value-bind (mags units allsame)
(pqval initial-contents)
(make-pq-ue
(call-next-method type dimensions rest-spec :initial-contents mags)
units
nil
allsame)))
(initial-element-p
(multiple-value-bind (mag units)
(pqval initial-element)
(make-pq-ue
(call-next-method type dimensions rest-spec :initial-element mag)
units
nil
t)))
(t (call-next-method))))
(defun replace-pq
(object new-units
&optional allsame
(new-system *system-of-units*) (old-system *system-of-units*))
"Replace the units on the object, using the numerical value given
in the current system of units."
(let ((*zero-is-dimensionless* nil))
(make-pq (pqval object old-system) new-units new-system (not allsame))))
;;;;****************************************************************************
;;;; Grid methods
;;;;****************************************************************************
(defmethod grid:gref ((object physical-quantity) &rest indices)
(make-pq (apply 'grid:gref (pq-magnitude object) indices)
(if (scalar-units object)
(pq-dimension object)
(apply 'grid:gref (pq-dimension object) indices))))
;;;;****************************************************************************
;;;; Bare sequences
;;;;****************************************************************************
(defun make-pq-seq
(magnitude-seq unit-expression-seq &optional (sysunits *system-of-units*))
"Make a sequence of PQs when unit-expression-seq is a sequence
of units that corresponds to the sequence of magnitudes."
(map
(if (listp magnitude-seq) 'list 'vector)
(lambda (m u) (make-pq m u sysunits))
magnitude-seq unit-expression-seq))
(defmethod make-pq-int
((magnitude sequence) unit-expression &optional (sysunits *system-of-units*))
(map (let ((type (type-of magnitude)))
(if (listp type)
(if (eql (first type) 'simple-array)
'simple-vector
(first type))
type))
(lambda (m) (make-pq m unit-expression sysunits)) magnitude))
(defmethod make-pq-int
((magnitude array) unit-expression
&optional (sysunits *system-of-units*))
(case (array-rank magnitude)
(1 (call-next-method))
(2
(let ((arr (make-array (array-dimensions magnitude))))
(loop for i from 0 below (array-dimension magnitude 0)
do
(loop for j from 0 below (array-dimension magnitude 1)
do
(setf (aref arr i j)
(make-pq (aref magnitude i j) unit-expression sysunits))))
arr))
(t (error "Can't make-pq arrays of rank higher than 2 yet."))))
#|
;;; (undimensioned-should-not-be-here) not defined yet
;;; Special handling of arrays with pqs
(defun *pq
(matrix vector &optional return-units)
"Multiple the matrix by the vector, taking care to remove the units
from the matrix and vector and add units on the result vector,
by default the same units the original vector had."
(undimensioned-should-not-be-here)
(multiple-value-bind (value units allsame)
(pqval vector)
(if allsame
(make-pq (antik:* (pqval matrix) value) (or return-units units))
(make-pq-seq (antik:* (pqval matrix) value) (or return-units units)))))
|#
;; Datime formats, compatibility with non-iso8601.
;; Liam Healy Fri Oct 12 2001 - 16:27
;; Time-stamp: <2010-12-25 22:15:18EST formats.lisp>
;; Time-stamp: <2011-01-01 17:14:16EST formats.lisp>
;; Copyright 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -81,7 +81,7 @@
"Find the timepoint from the Julian date.
Assumes a fixed number of seconds per day (86400)."
(multiple-value-bind (jdint jdfrac)
(floor jd)
(cl:floor jd)
(datime+
(make-timepoint
:dtspec (dtspec-from-julian-day-number jdint modified)
......
;; Objects that represent physical measurements.
;; Liam Healy Wed Mar 6 2002 - 09:04
;; Time-stamp: <2010-12-25 16:55:54EST physical-quantities.lisp>
;; Time-stamp: <2011-01-02 18:44:59EST physical-quantities.lisp>
;; Copyright 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -20,123 +20,7 @@
(in-package :antik)
(export '(make-pq pqval check-dimension
new-units with-pq define-physical-constant))
;;;;****************************************************************************
;;;; The object and making it
;;;;****************************************************************************
(defclass physical-quantity ()
((magnitude :initarg :magnitude :reader pq-magnitude :type number)
(dimension :initarg :dimension :reader pq-dimension :type list))
(:documentation "A quantity with a dimension and possibly a unit, e.g. length."))
(defparameter *zero-is-dimensionless* t
"Numbers with zero magnitude are made dimensionless.")
(defun make-pq-dimel (mag dimel)
"Make a physical quantity directly from dimel, unless dimensionless, in which
case the plain number is returned."
(if (or (equal dimel (dimension 'dimensionless))
(when *zero-is-dimensionless* (zerop mag)))
mag
(nth-value
0
(pqwu
(make-instance 'physical-quantity
:magnitude mag
:dimension dimel)))))
;;; (make-pq 33.4 "meters/second")
;;; (make-pq 88 '(/ feet minute))
;;; (make-pq 5 'speed)
;;; #_5.0_KM/S
;;; Use rectangles for arrays?
;;; Map into structures?
(defun make-pq
(magnitude unit-expression &optional (sysunits *system-of-units*) seq)
"Make the physical quantity. Allow for the possiblity that magnitude
is already a pq; in this case, check that the physical dimension
is correct, and pass it on. In this case, this function acts as a way
of assigning default units to a number. If sysunits is specified,
use units in that system instead of what are given in unit-expression."
(when (complexp magnitude)
(cerror "Accept" "Magnitude ~a is complex." magnitude))
(if (and seq (typep magnitude 'sequence))
(make-pq-seq magnitude unit-expression sysunits)
(make-pq-int magnitude unit-expression sysunits)))
(defun make-pq-zero (unit-expression &optional (sysunits *system-of-units*))
(let ((*zero-is-dimensionless* nil)) (make-pq 0 unit-expression sysunits)))
(defgeneric make-pq-int
(magnitude unit-expression &optional sysunits)
(:documentation "Make a physical quantity by class of magnitude."))
(defmethod make-pq-int
((magnitude number) unit-expression &optional (sysunits *system-of-units*))
(multiple-value-bind (el convfact)
(parse-units unit-expression sysunits)
(make-pq-dimel
(if convfact
(antik:* convfact magnitude)
(error
"~a is not a unit of measure but a physical dimension ~
and cannot be interpreted as a unit of measure ~
without a system of units."
unit-expression))
el)))
(defmethod make-pq-int
((pq physical-quantity) unit-expression &optional (sysunits *system-of-units*))
(check-dimension pq unit-expression)
pq)
(defmethod make-pq-int
((magnitude sequence) unit-expression &optional (sysunits *system-of-units*))
(map (let ((type (type-of magnitude)))
(if (listp type)
(if (eql (first type) 'simple-array)
'simple-vector
(first type))
type))
(lambda (m) (make-pq m unit-expression sysunits)) magnitude))
(defmethod make-pq-int
((magnitude array) unit-expression
&optional (sysunits *system-of-units*))
(case (array-rank magnitude)
;; These are hacks and should use dynamic rectangles
(1 (call-next-method))
(2 ; hack! only two dims
(let ((arr (make-array (array-dimensions magnitude))))
(loop for i from 0 below (array-dimension magnitude 0)
do
(loop for j from 0 below (array-dimension magnitude 1)
do
(setf (aref arr i j)
(make-pq (aref magnitude i j) unit-expression sysunits))))
arr))
(t (error "Can't make-pq arrays of rank higher than 2 yet."))))
(defun make-pq-seq
(magnitude-seq unit-expression-seq &optional (sysunits *system-of-units*))
"Make a sequence of PQs when unit-expression-seq is a sequence
of units that corresponds to the sequence of magnitudes."
(map
(if (listp magnitude-seq) 'list 'vector)
(lambda (m u) (make-pq m u sysunits))
magnitude-seq unit-expression-seq))
(defun replace-pq
(object new-units
&optional allsame
(new-system *system-of-units*) (old-system *system-of-units*))
"Replace the units on the object, using the numerical value given
in the current system of units."
(let ((*zero-is-dimensionless* nil))
(make-pq (pqval object old-system) new-units new-system (not allsame))))
(export '(pqval check-dimension new-units with-pq define-physical-constant))
;;;;****************************************************************************
;;;; Parsing unit names
......@@ -251,8 +135,11 @@
(nf val stream) ; format the number only
(format-units val units stream)))))
;;; eliminate?
(defmethod nf :around
((object sequence) &optional (stream *standard-output*))
(call-next-method)
#+(or)
(multiple-value-bind (values units all-same)
(pqval object)
(if (and units all-same)
......@@ -269,7 +156,7 @@
(:documentation
"Get the numerical value and the units of the physical quantity.
Returns the magnitude, the units, and whether units were
all the same for a sequence or array.")
all the same for a sequence or grid.")
(:method (pq &rest sysunits-additional-units)
(declare (ignore sysunits-additional-units))
(values pq nil nil)))
......
;; Do calculations without physical dimension.
;; Liam Healy Wed Jun 22 2005 - 22:35
;; Time-stamp: <2010-12-25 16:02:12EST undimension.lisp>
;; $Id$
;; Time-stamp: <2011-01-02 11:27:36EST undimension.lisp>
;; Copyright 2011 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -78,10 +77,11 @@
(*remove-units-sysunits* nil))
(make-pq magnitude unit-expression sysunits T)))
;;; Generalize to multiple value returns?
(defmacro dimensionless-call (function answer-dimension &rest arguments)
"Call a function that is incapable of handling physical quantities.
This macro will consistently remove the units from all arguments,
then the answer with the specified units."
then create the answer with the specified units."
`(without-units ()