share/physics/dimension.mac

   1 /*
   2 Copyright (C) 2002, 2008 Barton Willis
   3 Brief Description: Maxima code for dimensional analysis.
   4 I release this file under the terms of the GNU General Public License, version 2
   5 */
   6
   7 put('dimension,2,'version)$
   8
   9 /*
  10 If e is a list, return true if either e is empty or if e[i] = e[j]
  11 for all 1 <= i,j <= length(e). If e isn't a list, return false.
  12 */
  13
  14 all_equalp (e) := listp(e) and cardinality(setify(e))< 2;
  15
  16 /*
  17 fundamental_dimensions is a Maxima list that is used by the functions
  18 dimensionless, natural_unit, and dimension_as_list; the function
  19 dimension doesn't use this list.  A user may insert or remove elements
  20 from the list; for example, to use quantities that involve
  21 temperature, you could append temp onto fundamental_dimensions.
  22
  23 You'll get into trouble if you call any one of the functions
  24 dimensionless, natural_unit, or  dimension_as_list with expressions
  25 that have dimensions not in the list fundamental_dimensions.
  26 Be careful.
  27
  28 The default value of fundamental_dimensions is [mass, length, time].
  29 You can change this list to anything you like.
  30 */
  31
  32 assume('mass > 0);
  33 assume('length > 0);
  34 assume('time > 0);
  35 fundamental_dimensions : ['mass, 'length, 'time];
  36
  37 /* Notes
  38
  39 0.  When e is a list, map dimension over the list.
  40
  41 1.  Use logcontract so that log(a) - log(b) is okay when a and b have the
  42     same dimensions.
  43
  44 2.5  dimension(a[x]) = dimension(a) regardless of what x is.
  45
  46 3.  In a sum, check that all terms have the same dimensions; if not,
  47     signal the error "The expression is dimensionally inconsistent."
  48
  49 4.  abs, sqrt, max, and min  are the only prefix functions that may
  50     take an argument that isn't dimensionless.
  51
  52 5.  Check that all function arguments are dimensionless; if not,
  53     signal the error "The expression is dimensionally inconsistent."
  54
  55 Return the dimensions of the expression e; if e is dimensionally
  56 inconsistent, signal an error "The expression is dimensionally inconsistent."
  57
  58 The infix binary operators +, *, /, ^, ^^, ., = and ~ are supported.
  59 (Maxima's vector package defines ~ to be the cross product
  60 operator.)  For the = operator, the expression is dimensionally
  61 inconsistent if the dimension of the right and left sides differ.
  62
  63 It supports the prefix unary negation and the abs function. The nary
  64 min and max functions are supported; for the min and max functions,
  65 all arguments must have the same dimension, otherwise, the
  66 expression is dimensionally inconsistent.
  67 */
  68
  69
  70 dimension (e) := block([inflag : true, e_op, dim, var, i, n],
  71
  72   if listp(e) then (                     /* 0 */
  73      map ('dimension, e))
  74    else (
  75      e : logcontract (e),                 /* 1 */
  76      if not mapatom(e) then e_op : op(e),
  77      if constantp (e) then (
  78        1)
  79      else if symbolp (e) then (
  80        if member(e, fundamental_dimensions) then (
  81          e)
  82        else (
  83          dim : get (e, 'dimension),
  84          if dim = false then (
  85            funmake ('dimension, [e]))
  86          else (
  87            dim)))
  88      else if subvarp (e) then (          /* 2.5 */
  89        dimension(inpart(e,0)))
  90      else if e_op = "*" or e_op = "." or e_op = "~" then (
  91        xreduce("*", map ('dimension, args (e))))
  92      else if e_op = "/" then (
  93        dimension(num(e)) / dimension(denom(e)))
  94      /* apply (e_op, map ('dimension, args (e)))) */
  95      else if (e_op = "^" or e_op = "^^") and dimension (inpart (e,2)) = 1 then (
  96        dimension(inpart(e,1)) ^ inpart(e,2))
  97      else if e_op = "-" then (
  98        dimension(inpart(e,1)))
  99      else if e_op = "+" or e_op = 'max or e_op =  'min then ( /*3*/
 100        dim : map ('dimension, args (e)),
 101        if all_equalp (dim) then (
 102          first (dim))
 103        else (
 104          error ("The expression is dimensionally inconsistent.")))
 105      else if e_op = 'matrix then (
 106        dim : map ('dimension, args (e)),
 107        dim : apply ('append, dim),
 108        if all_equalp (dim) then (
 109          first(dim))
 110        else (
 111          error ("The expression is dimensionally inconsistent.")))
 112      else if (e_op = "="  or e_op = "#" or e_op = "<" or e_op = ">" or e_op = "<=" or e_op = ">=" ) then (
 113        if inpart(e,1) = 0 then (
 114          dimension (inpart(e,2)))
 115        else if inpart (e,2) = 0 then (
 116          dimension (inpart(e,1)))
 117        else (
 118          dim : dimension(inpart(e,1)),
 119          if dim = dimension(inpart(e,2)) then (
 120            dim )
 121          else (
 122            error ("The expression is dimensionally inconsistent."))))
 123      else if e_op = 'abs then (     /*  4 */
 124        dimension (inpart (e, 1)))
 125      else if e_op = 'sqrt then (
 126        dimension (inpart(e,1)) ^ (1/2))
 127      else if e_op = nounify('diff) then (
 128        var : args(e),
 129        n : length(var) - 1,
 130        dim : dimension (inpart(var,1)),
 131        for i : 2 thru n step 2 do (
 132          dim : dim / dimension (inpart(var,i)) ^ inpart(var,i+1)),
 133        dim)
 134      else if e_op = nounify('integrate) then (
 135        var : args(e),
 136        dim : dimension(inpart(var,1)) * dimension(inpart(var,2)),
 137        if length (e) = 4 and not all_equalp([dimension(inpart(var,2)),
 138          dimension(inpart(var,3)),dimension(inpart(var,4))]) then (
 139          error ("The expression is dimensionally inconsistent."))
 140        else (
 141          dim))
 142      else if e_op = nounify('sum) or e_op = nounify('product) then (
 143        error("The function dimension doesn't handle sums or products."))
 144      else if e_op = "DIV" or  e_op = "div" or e_op = "GRAD" or e_op = "grad"  or
 145      e_op = "CURL" or e_op = "curl" then (
 146        if listp(inpart(e,1)) then (
 147          dim : map ('dimension, inpart(e,1)),
 148          if all_equalp(dim) then (
 149            dim : first(dim))
 150          else (
 151            error ("The expression is dimensionally inconsistent.")))
 152        else (
 153          dim : dimension (inpart(e,1))),
 154        dim / 'length)
 155      else if e_op = "laplacian" then (
 156        dimension(inpart(e,1) / 'length^2))
 157      else (                                /*  5 */
 158        dim : map ('dimension, args(e)),
 159        if all_equalp (cons (1,dim)) then (
 160          dim : dimension(inpart(e,0)),
 161          if not atom(dim) and ?equal(inpart(dim,0),'dimension) then (
 162            1)
 163          else (
 164            dim ))
 165        else (
 166          error ("The expression is dimensionally inconsistent.")))));
 167
 168 /*
 169 For the expression e, return a list of the exponents of the
 170 fundamental_dimensions. Thus if  fundamental_dimensions =
 171 [mass, length, time] and c is a velocity, dimension_as_list(c)
 172 returns [0,1,-1].
 173 */
 174
 175 dimension_as_list(e) := (
 176   if listp(e) then map('dimension_as_list, e)
 177   else (
 178     e : dimension(e),
 179     if polynomialp(e, fundamental_dimensions, lambda([x], numberp(x)), lambda([x], integerp(x))) then
 180     makelist(hipow(e, dk), dk, fundamental_dimensions)
 181     else error("Unable to determine the dimension of", e)));
 182
 183 /*
 184 Return a basis for the dimensionless quantities that can be formed from a product of powers of elements
 185 of the list e. The basis excludes constants.
 186 */
 187
 188 dimensionless(e) := (
 189   require_list(e,"first", "dimensionless"),
 190   args(map(lambda([s], xreduce("*", map("^", e, first(transpose(s))))),
 191       nullspace(transpose(funmake('matrix, map('dimension_as_list, e)))))));
 192
 193 natural_unit(dim,e) := block([vars, s, linsolve_params : true,
 194   back_subst : true, globalsolve : false],
 195   require_list(e,"second", "natural_unit"),
 196
 197  dim : dimension_as_list(dim),
 198  s : map ('dimension_as_list, e),
 199  vars : makelist(?gensym(),k,1,length(e)),
 200  s : linsolve(s . vars - dim, vars),
 201  s : xreduce("*", map("^", e, map('rhs, s))),
 202  s : subst(map("=", %rnum_list, makelist(0,i,1,length(%rnum_list))), s),
 203  map(lambda([x], s * x), cons(1, dimensionless(e))));
 204
 205
 206
 207