share/ezunits/ezunits_functions.mac

   1 infix ("`", 118, 118);
   2 infix ("``", 117, 117);
   3
   4 declare (constvalue, evfun);
   5
   6 constvalue (e) :=
   7   if mapatom(e)
   8     then
   9      (if symbolp(e)
  10         then (get (e, 'constvalue), if %% # false then %% else e)
  11         else e)
  12     else map (constvalue, e);
  13
  14 declare_constvalue (a, x) ::= buildq ([a, x], put (a, x, 'constvalue));
  15
  16 remove_constvalue (a) ::= buildq ([a], rem (a, 'constvalue));
  17
  18 /* Hold off on this display symbol business.
  19  * :lisp (put '$\` '(#\space #\space) 'dissym)
  20  */
  21
  22 declare (dimensional, feature);
  23 declare (nondimensional, feature);
  24
  25 declare_units (a%, u%) :=
  26   block (local (a%, u%),
  27   if listp(a%)
  28     then map (lambda ([a1%], local(a1%), declare1_units (a1%, u%)), a%)
  29     else declare1_units (a%, u%), u%);
  30
  31 declare1_units (a%, u%) :=
  32   block (local (a%, u%),
  33   apply ('define_variable, [a%, a%, any_check]),
  34   put (a%, units_check (a%), value_check),
  35   put (a%, u%, units),
  36   if featurep (a%, 'dimensional) = false
  37     then apply (declare, [a%, 'dimensional]));
  38
  39 declare_qty (a%, q%) :=
  40   block (local (a%, q%),
  41   if listp(a%)
  42     then map (lambda ([a1%], local(a1%), declare1_qty (a1%, q%)), a%)
  43     else declare1_qty (a%, q%), q%);
  44
  45 /* Maybe we want to work with constvalue in declare1_qty
  46  * (instead of the indicator qty).
  47  */
  48 declare1_qty (a%, q%) :=
  49   block (local (a%, q%),
  50   if featurep (a%, 'dimensional) = false
  51     then apply (declare, [a%, 'dimensional]), put (a%, q%, qty));
  52
  53 units_check (zz%) :=
  54   block (local (zz%),
  55   buildq ([zz%: zz%], lambda([yy%],
  56     if units('zz%) # 1 and units(yy%) # 1 and units(yy%) # units('zz%)
  57     then throw (oops ('zz%, yy%)))));
  58
  59 /* PROBABLY QTY SHOULD BE A SIMPLIFYING FUNCTION !! OH WELL. KEEP STUMBLING AHEAD. */
  60
  61 qty(e%) :=
  62   block (local (e%),
  63   if atom(e%) then
  64     if featurep (e%, 'dimensional) then
  65       if get (e%, qty) = false then 'qty(e%) else get (e%, qty)
  66     else e%
  67   else if unitop_p(e%) then first(e%)
  68   else if featurep (op(e%), 'dimensional) then 'qty(e%)
  69   else if op(e%) = 'qty or op(e%) = nounify('qty) then e%
  70   else apply (op(e%), map (qty, args(e%))));
  71
  72 units(e%) :=
  73   block (local (e%),
  74   if atom(e%) then
  75     if symbolp (e%) and get (e%, 'constvalue) # false
  76       then units (get (e%, 'constvalue))
  77     elseif featurep (e%, 'dimensional) then
  78       if get (e%, units) = false then 'units(e%) else get (e%, units)
  79     else 1
  80   else if mapatom (e%) then units (op (e%))  /* mapatom but not atom => subscripted */
  81   else if unitop_p(e%) then second(e%)
  82   else if featurep (op(e%), 'dimensional) then units(op(e%))
  83   else if op(e%) = "*" or op(e%) = "/" then map (units, e%)
  84   else if op(e%) = "^" then units (first (e%)) ^ second (e%)
  85   else if op(e%) = 'sqrt then sqrt (units (e%))
  86   else if op(e%) = "+" or op(e%) = "-" then apply ("+", unique (map (units, args (e%))))
  87   else if member (op(e%), ["=", 'equal, "[", 'matrix, "{"]) then map (units, e%)
  88   else if op(e%) = "." then apply ("*", map (units, args (e%)))
  89   else if op(e%) = nounify(diff)
  90     then units_diff(e%)
  91   else if op(e%) = nounify(integrate)
  92     then units_integrate(e%)
  93   else if op(e%) = nounify(at)
  94     then units_at(e%)
  95   else if nondimensionalp (e%) then 1
  96   else 'units(e%));
  97
  98 units_diff(e%) :=
  99   block ([yy: first(e%), xnpairs: by_twos (rest (args (e%)))],
 100          lproduct (units(xn[1])^xn[2], xn, xnpairs),
 101          units(yy) / %%);
 102
 103 /* by_twos might want to be a built-in function ... */
 104 by_twos (l) :=
 105   block ([n: length(l)],
 106          if mod (n, 2) = 0
 107            then makelist ([l[2*k - 1], l[2*k]], k, 1, n/2)
 108            else error ("by_twos: length of list must be a multiple of 2, not", n));
 109
 110 /* define lproduct as a macro to work around funny evaluation of makelist ...
 111  * lproduct might want to be a built-in function ...
 112  */
 113 lproduct (ee, xx, l) ::=
 114   buildq ([ee, xx, l], apply ("*", makelist (ee, xx, l)));
 115
 116 units_integrate(e%) :=
 117   block ([ii: first(e%), xx: second(e%)],
 118          units(ii) * units(xx));
 119
 120 units_at(e%) := units (first (e%));
 121
 122 d(e%) := dimension (units (e%));
 123
 124 matchdeclare
 125   ([aa%, bb%], true,
 126     cc%, constantp,
 127     dd%, lambda ([e], atom(e) and featurep (e, 'dimensional)),
 128     ccn0%, constantp_not0,
 129     ccn1%, constantp_not1,
 130     /* dl = "dimensionless" */
 131     dl%, lambda ([e], dimensions (e) = 1),
 132     mm%, mult_expr_nontrivialconstfactorsp,
 133     nd%, nondimensionalp,
 134     pp%, lambda ([e], not atom(e) and op(e) = "+"),
 135     qq%, lambda ([e], not atom(e) and (op(e) = "=" or listp(e) or matrixp(e) or setp(e))),
 136     uu%, unitop_p,
 137     xx%, nonconstantp);
 138
 139 fooexpr : nd% / pp%; /* simplify this, to be used later */
 140 simp: false$
 141 /* tellsimpafter rules are applied in order of declaration; tellsimp, the other way around */
 142 tellsimp (nd% `` dl%, apply_units_conversion (nd%, dl%));
 143 tellsimp (uu% `` aa%, apply_units_conversion (uu%, aa%));
 144 tellsimp (aa%*pp% `` bb% , expand (aa%*pp%) `` bb%); /* HOW ABOUT THIS ONE ?? */
 145 tellsimp (pp% `` aa%, apply ("+", map (buildq ([aa%], lambda ([e], e `` aa%)), args (pp%))));
 146 tellsimp (dd% `` aa%, (qty(dd%) ` units(dd%)) `` aa%);
 147 tellsimp (aa%*dd% `` bb%, aa% * (qty(dd%) ` units(dd%)) `` bb%);
 148 tellsimp (qq% `` aa%, map (lambda ([e], e `` aa%), qq%));
 149 tellsimp (''fooexpr `` aa%, nd%/(pp% `` 1/aa%));
 150 tellsimpafter (uu%^ccn1%, (qty(uu%)^ccn1%) ` simplify_units (units(uu%)^ccn1%));
 151 tellsimpafter (aa% * nd% * uu%, aa% * multiply_with_units (nd%, uu%));
 152 tellsimpafter (aa% + uu%, add_with_units (aa%, uu%));
 153
 154 tellsimpafter (aa% ` (xx% + ccn0%), (aa% - ccn0%)`xx%);
 155 tellsimpafter (aa% ` mm%, (constant_factors(mm%) * aa%) ` everything_else(mm%));
 156 tellsimpafter (aa% ` cc%, aa%*cc%);
 157 tellsimpafter (uu% ` nd%, qty (uu%) ` (units (uu%) * nd%));
 158 tellsimpafter (abs (aa% ` bb%), abs (aa%) ` bb%);
 159 tellsimpafter (aa% ` abs (bb%), aa% ` bb%);
 160
 161 /* RULES FOR RELATIONAL EXPRESSIONS HERE !! NOTIFY GUNTER KOENIGSMANN */
 162 matchdeclare (uu%, unitp, vv%, lambda ([e], unitp(e) and dimensions(e) = dimensions(uu%)));
 163 tellsimpafter (uu% < vv%, qty (uu%) < qty (vv% `` units (uu%)));
 164 tellsimpafter (uu% <= vv%, qty (uu%) <= qty (vv% `` units (uu%)));
 165 tellsimpafter (uu% >= vv%, qty (uu%) >= qty (vv% `` units (uu%)));
 166 tellsimpafter (uu% > vv%, qty (uu%) > qty (vv% `` units (uu%)));
 167 /* SIMPLIFICATION OF DIMENSIONAL "=" CAUSES TROUBLE FOR EXISTING FUNCTIONALITY !!
 168  * ALSO IT SEEMS REALLY HEAVYHANDED TO RULE OUT FOO`M = BAR`FOOT !!
 169 tellsimpafter (uu% = vv%, qty (uu%) = qty (vv% `` units (uu%)));
 170 tellsimpafter (uu% # vv%, qty (uu%) # qty (vv% `` units (uu%)));
 171  */
 172 tellsimpafter (equal (uu%, vv%), equal (qty (uu%), qty (vv% `` units (uu%))));
 173 tellsimpafter (notequal (uu%, vv%), notequal (qty (uu%), qty (vv% `` units (uu%))));
 174 /* I'M NOT COMFORTABLE WITH IMPLICIT CONVERSION OF 0 TO DIMENSIONAL !!
 175  * I REMEMBER NOW THAT THE JUSTIFICATION IS THAT 0 IS AN ELEMENT OF EVERY DIMENSIONAL LINE !!
 176  * IS THAT ENOUGH TO GO ON ??
 177 tellsimpafter (uu% < 0, qty (uu%) < 0);
 178 tellsimpafter (0 < uu%, 0 < qty (uu%));
 179 tellsimpafter (uu% <= 0, qty (uu%) <= 0);
 180 tellsimpafter (0 <= uu%, 0 <= qty (uu%));
 181 tellsimpafter (uu% >= 0, qty (uu%) >= 0);
 182 tellsimpafter (0 >= uu%, 0 >= qty (uu%));
 183 tellsimpafter (uu% > 0, qty (uu%) > 0);
 184 tellsimpafter (0 > uu%, 0 > qty (uu%));
 185 tellsimpafter (uu% = 0, qty (uu%) = 0);
 186 tellsimpafter (0 = uu%, 0 = qty (uu%));
 187 tellsimpafter (uu% # 0, qty (uu%) # 0);
 188 tellsimpafter (0 # uu%, 0 # qty (uu%));
 189 tellsimpafter (equal (uu%, 0), equal (qty (uu%), 0));
 190 tellsimpafter (equal (0, uu%), equal (0, qty (uu%)));
 191 tellsimpafter (notequal (uu%, 0), notequal (qty (uu%), 0));
 192 tellsimpafter (notequal (0, uu%), notequal (0, qty (uu%)));
 193  */
 194
 195 simp: true$
 196
 197 matchdeclare (uu%, lambda ([e], unitp (e) and not unitop_p (e)));
 198 defrule (rexpand_dimensional, uu%, qty (uu%) ` units (uu%));
 199 expand_dimensional (e) := applyb1 (e, rexpand_dimensional);
 200
 201 /* expand_dimensional('qty(xyz)) => 'qty('qty(xyz) ` units(xyz)) without this rule ... sigh !!
 202  * does the need for it go away if qty is a simplifying function ??
 203  */
 204 matchdeclare (uu%, unitop_p);
 205 tellsimpafter ('qty (uu%), qty (uu%));
 206
 207 unitop_p (e%) := not atom(e%) and atom (op(e%)) and is (nounify (op(e%)) = nounify ("`"));
 208 unitp (e%) :=
 209   if atom(e%) then featurep (e%, 'dimensional)
 210   else if unitop_p (e%) then true
 211   else featurep (op(e%), 'dimensional);
 212
 213 constantp_not1 (x%) := constantp(x%) and x% # 1;
 214 constantp_not0 (x%) := constantp(x%) and x% # 0;
 215 nonconstantp (x%) := not constantp (x%);
 216 mult_exprp (e%) := not atom(e%) and atom (op(e%)) and (nounify (op(e%)) = nounify ("*") or nounify (op(e%)) = nounify ("/"));
 217 mult_expr_nontrivialconstfactorsp (e%) := mult_exprp (e%) and constant_factors (e%) # 1;
 218 nondimensionalp (e%) :=
 219   if atom (e%) then (not symbolp (e%) or not featurep (e%, 'dimensional))
 220   elseif featurep (nounify (op (e%)), 'nondimensional) then true
 221   else
 222    (get (nounify (op (e%)), 'nondimensionalp),
 223     if %% # false then apply (%%, args (e%)));
 224 nondimensional_not1 (e%) := e% # 1 and nondimensionalp (e%);
 225 nondimensional_not0 (e%) := e% # 0 and nondimensionalp (e%);
 226
 227 (map (nounify, [qty, sin, cos, tan, asin, acos, atan, atan2, log, exp, sinh, cosh, tanh, asinh, acosh, atanh]),
 228   apply (declare, [%%, 'nondimensional]));
 229 /* HMM. MAYBE THESE RULES FOR "^" AND "/" ARE TOO STRICT. LET IT STAND FOR NOW. */
 230 put (nounify ("^"), lambda ([a%, b%], nondimensionalp (a%) and nondimensionalp (b%)), 'nondimensionalp);
 231 put (nounify ("/"), lambda ([a%, b%], nondimensionalp (a%) and nondimensionalp (b%)), 'nondimensionalp);
 232 put (nounify (sqrt), lambda ([a%], nondimensionalp (a%)), 'nondimensionalp);
 233 put (nounify ("+"), lambda ([[a%]], every (nondimensionalp, a%)), 'nondimensionalp);
 234 put (nounify ("*"), lambda ([[a%]], every (nondimensionalp, a%)), 'nondimensionalp);
 235 put (nounify ("-"), lambda ([[a%]], every (nondimensionalp, a%)), 'nondimensionalp);
 236
 237 add_with_units (x%, u%) :=
 238   if unitp (u%) then
 239     x% + u%
 240   else
 241     /* identify terms with like units and add them up. */
 242     block ([r%: [], a%: args(u%)],
 243     while a% # [] do
 244       block ([s%: qty(a%[1]), b%: []],
 245       for i%:2 thru length(a%) do
 246         if units (a%[i%]) = units (a%[1]) then
 247           s%: s% + qty (a%[i%])
 248         else
 249           b%: cons (a%[i%], b%),
 250       r%: cons (s%`units(a%[1]), r%),
 251       a%: b%),
 252     if length (r%) = 1 then
 253       x% + r%[1]
 254     else
 255       cons (x%, apply ("+", r%)));
 256
 257 multiply_with_units(x%, u%) :=
 258   if unitp (u%) then
 259     (x%*qty(u%))`units(u%)
 260   else
 261     /* Work around effect of *AFTERFLAG in tellsimpafter rules:
 262      * if units = 1, return foo instead of foo ` 1.
 263      */
 264     block ([qq, uu],
 265       qq : x%*apply ("*", map (qty, args (u%))),
 266       uu : simplify_units (apply ("*", map (units, args (u%)))),
 267       if uu = 1 then qq else qq ` uu);
 268
 269 matchdeclare (a%, all);
 270 defrule (runits1, 'units (a%), 1);
 271
 272 constant_factors (e%) := block ([L%, eqns%],
 273     e% : apply1 (e%, runits1),
 274     L%: listofvars(e%),
 275     eqns%: map ("=", L%, L%*0+1),
 276     subst (eqns%, e%));
 277
 278 everything_else(e%) := e%/constant_factors(e%);
 279
 280 simplify_units (e%) := block
 281   ([L% : listofvars (e%), mycontext],
 282     L% : makelist (x% > 0, x%, L%),
 283     mycontext : newcontext (),
 284     apply (assume, L%),
 285     e% : expand (e%, 0, 0),
 286     apply (forget, L%),
 287     killcontext (mycontext),
 288     e%);
 289
 290 known_unit_conversions : {};
 291
 292 get_unit_conversion (u) := assoc (u, known_unit_conversions);
 293
 294 declare_unit_conversion ([L]) := block ([S, ctxt],
 295     if not every (equationp, L)
 296         then error ("declare_unit_conversion: all arguments must be equations."),
 297     setify (ratsimp (L)),
 298     S : map (flatten_conversion_equation, %%),
 299     ctxt : context,
 300     context : ezunits_context,
 301     apply (assume, map (lambda ([x], x > 0), maplist (rhs, S))),
 302     known_unit_conversions : union (known_unit_conversions, S),
 303     kill (conversions_to_base_units), /* simple, heavy-handed */
 304     reset_convert_units (), /* simple, heavy-handed */
 305     context : ctxt,
 306     done);
 307
 308 flatten_conversion_equation (e%) := block
 309   ([l% : lhs (e%), r% : rhs (e%)],
 310     if unitop_p (l%) and unitop_p (r%)
 311         then qty (l%) * units (l%) = qty (r%) * units (r%)
 312         elseif not unitop_p (l%) and not unitop_p (r%) then e%
 313         else error ("declare_unit_conversion: both sides must be unit expressions or both not unit expressions; found: ", e%));
 314
 315 equationp (e) := not atom (e) and op (e) = "=";
 316
 317 compute_conversion_factor (u1%, u2%) :=
 318  (if conversions_to_base_units = 'conversions_to_base_units
 319     then compute_conversions_to_base_units (),
 320   ev (simplify_units (u1% / u2%), conversions_to_base_units));
 321
 322 compute_conversions_to_base_units () :=
 323   conversions_to_base_units : listify (compute_conversions_to_base_units0 (known_unit_conversions));
 324
 325 compute_conversions_to_base_units0 (conversion_eqns) :=
 326 block ([derived_units, primitive_units, defined_so_far, undefined_so_far, definitions, definablep, now_definable],
 327
 328   derived_units : map (lhs, conversion_eqns),
 329   primitive_units : setdifference (setify (listofvars (map (rhs, conversion_eqns))), derived_units),
 330   defined_so_far : primitive_units,
 331   undefined_so_far : conversion_eqns,
 332   definitions : {},
 333
 334   definablep : lambda ([e], every (lambda ([x], elementp (x, defined_so_far)), setify (listofvars (rhs (e))))),
 335   now_definable : subset (undefined_so_far, definablep),
 336
 337   while length (now_definable) > 0
 338     do block ([make_definition, new_definitions],
 339       make_definition : lambda ([e], block ([L : listofvars (rhs (e))],
 340                        for x in definitions do if member (lhs (x), L) then e : subst (x, e)), e),
 341
 342       new_definitions : map (make_definition, now_definable),
 343       definitions : union (definitions, new_definitions),
 344       defined_so_far : union (defined_so_far, map (lhs, new_definitions)),
 345       undefined_so_far : setdifference (undefined_so_far, now_definable),
 346       now_definable : subset (undefined_so_far, definablep)),
 347
 348   if undefined_so_far # {} then print ("ezunits: hmm, these are still undefined:", map (lhs, undefined_so_far)),
 349   definitions);
 350
 351 apply_units_conversion (uu%, nd%) :=
 352   if listp (nd%)
 353     /* Handle conversions of the form foo `` [u1, u2, u3].
 354      * May want to ensure dimensions of u1, u2, u3 are all the same.
 355      * Might change list to some other operator, e.g. noncommutative plus
 356      * (which would have to be invented for the purpose).
 357      */
 358     then block ([q%, r%],
 359         q% : convert_units [units (uu%), first (nd%)] (qty (uu%)),
 360         r% : [(if rest (nd%) = [] then q% else floor (q%)) ` first (nd%)],
 361         while rest (nd%) # []
 362             do
 363                (q% : convert_units [first (nd%), second (nd%)] (q% - floor (q%)),
 364                 nd% : rest (nd%),
 365                 r% : cons ((if rest (nd%) = [] then q% else floor (q%)) ` first (nd%), r%)),
 366         reverse (r%))
 367
 368     else convert_units [units(uu%), nd%] (qty(uu%)) ` nd%;
 369
 370 /* kill the convert_units memoizing function (including all memoized values),
 371  * then redefine it, and prep it with conversions for temperature units,
 372  * which can't be computed by a simple multiplication.
 373  * Hmm, I wonder if the temperature conversions should be declared in ezunits_functions.mac.
 374  */
 375
 376 reset_convert_units () :=
 377    (kill (convert_units),
 378
 379     convert_units [u1%, u2%] := block ([a%],
 380         a% : compute_conversion_factor (u1%, u2%),
 381         buildq ([a%], lambda ([x], a% * x))),
 382
 383     convert_units [degF, degC] : lambda ([x], 5/9*(x - 32)),
 384     convert_units [degF, K] : lambda ([x], 27315/100 + 5/9*(x - 32)),
 385     convert_units [degF, R] : lambda ([x], 45967/100 + x),
 386
 387     convert_units [degC, degF] : lambda ([x], 32 + 9/5*x),
 388     convert_units [degC, R] : lambda ([x], 45967/100 + 32 + 9/5*x),
 389     convert_units [degC, K] : lambda ([x], 27315/100 + x),
 390
 391     convert_units [K, degC] : lambda ([x], x - 27315/100),
 392     convert_units [K, degF] : lambda ([x], 32 + 9/5*(x - 27315/100)),
 393
 394     convert_units [R, degC] : lambda ([x], 5/9*(x - 45967/100 - 32)),
 395     convert_units [R, degF] : lambda ([x], x - 45967/100));
 396
 397 /* Dimensional analysis
 398  */
 399
 400 declare_dimensions ([L%]) :=
 401     if length (L%) > 2
 402         then (declare_dimensions1 (L%[1], L%[2]), apply (declare_dimensions, rest (rest (L%))))
 403         elseif length (L%) = 2
 404             then (declare_dimensions1 (L%[1], L%[2]), 'done)
 405             else error ("declare_dimensions: expected an even number of arguments.");
 406
 407 declare_dimensions1 (a%, b%) :=
 408     if listp (a%)
 409         then map (lambda ([a1%], declare_dimensions1 (a1%, b%)), a%)
 410         elseif symbolp (b%)
 411             then put (a%, b%, 'dimension)
 412             else error ("declare_dimensions: second argument must be a symbol.");
 413
 414 remove_dimensions ([L%]) :=
 415    (if not every (symbolp, L)
 416         then error ("remove_dimensions: all arguments must be symbols."),
 417     map (lambda ([a%], put (a%, false, 'dimension)), L%));
 418
 419 declare_fundamental_dimensions ([L]) :=
 420    (if not every (symbolp, L)
 421         then error ("declare_fundamental_dimensions: all arguments must be symbols."),
 422     L : sublist (L, lambda ([a], not member (a, fundamental_dimensions))),
 423     map (lambda ([a], assume (a > 0)), L),
 424     fundamental_dimensions : append (fundamental_dimensions, L),
 425     'done);
 426
 427 remove_fundamental_dimensions ([L]) :=
 428     if L = '[all]
 429         then apply (remove_fundamental_dimensions, fundamental_dimensions)
 430         else
 431            (fundamental_dimensions : sublist (fundamental_dimensions, lambda ([x], not member (x, L))),
 432             map (lambda ([a], forget (a > 0)), L),
 433             'done);
 434
 435 declare_fundamental_units ([L]) :=
 436    (if not every (symbolp, L)
 437         then error ("declare_fundamental_units: all arguments must be symbols."),
 438     map (lambda ([x, y], put (x, y, 'dimension)), odds (L), evens (L)),
 439     map (lambda ([x, y], put (y, x, 'unit)), odds (L), evens (L)));
 440
 441 remove_fundamental_units ([L]) :=
 442    (if not every (symbolp, L)
 443         then error ("remove_fundamental_units: all arguments must be symbols."),
 444     map (lambda ([x], put (get (x, 'dimension), false, 'unit)), L),
 445     map (lambda ([x], put (x, false, 'dimension)), L));
 446
 447 fundamental_dimensions : [];
 448 block ([L],
 449     L : [[m, length], [kg, mass], [s, time], [A, current], [K, temperature], [mol, quantity], [cd, luminous_intensity]],
 450     apply (declare_fundamental_dimensions, map (second, L)),
 451     apply (declare_fundamental_units, flatten (L)));
 452
 453 /* I GUESS RULES FOR DIMENSIONS COULD RECOGNIZE EQUATIONS, DERIVATIVES, AND INTEGRALS
 454  */
 455
 456 matchdeclare
 457    (aa%, lambda ([e%], numberp (e%) or (atom (e%) and not symbolp (e%))),
 458     bb%, lambda ([e%], symbolp (e%) and member (e%, fundamental_dimensions)),
 459     cc%, lambda ([e%], symbolp (e%) and get (e%, 'units) # false),
 460     dd%, lambda ([e%], symbolp (e%) and get (e%, 'constvalue) # false),
 461     ee%, lambda ([e%], symbolp (e%) and get (e%, 'dimension) # false),
 462     ee1%, lambda ([e%], symbolp (e%) and ?mget (e%, 'numer) # false),
 463     ff%, lambda ([e%], symbolp (e%)),
 464     gg%, lambda ([e%], not atom (e%) and op (e%) = "`"),
 465     hh%, lambda ([e%], not atom (e%) and op (e%) = "^"),
 466     ii%, lambda ([e%], not mapatom (e%) and member (op (e%), ["+", "-", "*", "/"])));
 467
 468 simp : false;
 469 defrule (raa, 'dimensions (aa%), 1);
 470 defrule (rbb, 'dimensions (bb%), bb%);
 471 defrule (rcc, 'dimensions (cc%), dimensions (get (cc%, 'units)));
 472 defrule (rdd, 'dimensions (dd%), dimensions (get (dd%, 'constvalue)));
 473 defrule (ree, 'dimensions (ee%), get (ee%, 'dimension));
 474 defrule (ree1, 'dimensions (ee1%), dimensions (?mget (ee1%, 'numer)));
 475 defrule
 476    (rff,
 477     'dimensions (ff%),
 478     block
 479       ([ff%2 : ev (ff%, args (known_unit_conversions), infeval)],
 480         if ff%2 # ff% then dimensions (ff%2) else 'dimensions (ff%)));
 481 defrule (rgg, 'dimensions (gg%), (constvalue (first (gg%)), if unitp (%%) then dimensions (second (%%)) else 1) * dimensions (second (gg%)));
 482 matchdeclare (aa%, all, bb%, constantp_not1);
 483 defrule (rhh, 'dimensions (aa%^bb%), dimensions (aa%) ^ bb%);
 484 defrule (rii, 'dimensions (ii%), map (dimensions, ii%));
 485 /* ??? defrule (rjj, 'dimensions (sqrt (aa%)), sqrt (dimensions (aa%))); */
 486 simp : true;
 487
 488 dimensions (expr) :=
 489     /* Binding *NOUNSFLAG* is pretty terrible,
 490      * but it's the only way to disable ev(..., nouns);
 491      * ev(..., nouns=false) has no effect!
 492      */
 493     block ([?\*nounsflag\*: false],
 494            dimensions1 (expr));
 495
 496 dimensions1 (expr) :=
 497     if not atom (expr) and (setp (expr) or listp (expr) or matrixp (expr) or op (expr) = "=")
 498         then map (dimensions, expr)
 499         else apply (apply1, ['dimensions (expr), raa, rbb, rcc, rdd, ree, ree1, rff, rgg, rhh, rii]);
 500
 501 /* Adapted from dimensions_as_list in share/physics/dimension.mac.
 502  * Thanks to Barton Willis.
 503  */
 504 /* FOLLOWING DEFINITION HANDLES DIMENSION NOUNS INCORRECTLY
 505  * NEED TO REVISIT THIS
 506  */
 507 dimensions_as_list (e) :=
 508     if listp (e)
 509         then map ('dimensions_as_list, e)
 510         else
 511            (e : dimensions (e),
 512             if polynomialp (e, fundamental_dimensions, lambda ([x], numberp (x)), lambda ([x], integerp (x)))
 513                 then makelist (hipow (e, d), d, fundamental_dimensions)
 514                 else 'dimensions_as_list (e));
 515
 516 fundamental_units ([L]) :=
 517     if L = []
 518         then
 519             map
 520                (lambda ([x], (get (x, 'unit), if %% = false then 'unit (x) else %%)),
 521                 fundamental_dimensions)
 522         else
 523             if listp (L [1])
 524                 then map ('fundamental_units, L [1])
 525                 else block ([d],
 526                     d : dimensions_as_list (L [1]),
 527                     if listp (d)
 528                         then apply ("*", map ("^", fundamental_units (), d))
 529                         else 'fundamental_units (L [1]));
 530
 531 /* Adapted from dimensionless in share/physics/dimension.mac.
 532  * Thanks to Barton Willis.
 533  */
 534
 535 dimensionless (e) :=
 536     if not listp (e)
 537         then 'dimensionless (e)
 538         else
 539             args
 540                (map
 541                    (lambda ([s], xreduce ("*", map("^", e, first (transpose (s))))),
 542                     nullspace (transpose (funmake ('matrix, map ('dimensions_as_list, e))))));
 543
 544 /* Adapted from natural_unit in share/physics/dimension.mac.
 545  * Thanks to Barton Willis.
 546  */
 547
 548 natural_unit (d, e) :=
 549 block
 550   ([vars, s],
 551
 552     if not listp (e)
 553         then return ('natural_unit (d, e)),
 554
 555     d : dimensions_as_list (d),
 556     s : map ('dimensions_as_list, e),
 557     vars : makelist (?gensym(), k, 1, length (e)),
 558
 559     block
 560       ([linsolve_params : true,
 561         back_subst : true,
 562         globalsolve : false],
 563         s : linsolve (s . vars - d, vars)),
 564
 565     s : xreduce ("*", map ("^", e, map ('rhs, s))),
 566     s : subst (map ("=", %rnum_list, makelist (0, i, 1, length (%rnum_list))), s),
 567     map (lambda ([x], s * x), cons (1, dimensionless (e))));
 568
 569 /* A couple of convenience functions.
 570  */
 571
 572 derived_units () :=
 573    (if conversions_to_base_units = 'conversions_to_base_units
 574         then compute_conversions_to_base_units (),
 575     listify (map (lhs, conversions_to_base_units)));
 576
 577 known_units () :=
 578     sort (append (fundamental_units (), derived_units ()));
 579
 580 display_known_unit_conversions () :=
 581     for x in known_unit_conversions do display (x);
 582
 583 /* Rules to unit-ify 'diff and 'integrate expressions.
 584  */
 585
 586 simp : false;
 587
 588 matchdeclare ([uu%, vv%, aa%, bb%], unitop_p);
 589 matchdeclare (nn%, integerp);
 590
 591 tellsimp ('diff (uu%, vv%, nn%), 'diff (qty (uu%), qty (vv%), nn%) ` (units (uu%) / (units (vv%))^nn%));
 592
 593 tellsimp ('integrate (uu%, vv%), 'integrate (qty (uu%), qty (vv%)) ` (units (uu%) * units (vv%)));
 594
 595 tellsimp
 596    ('integrate (uu%, vv%, aa%, bb%),
 597     if units (aa%) = units (vv%) and units (bb%) = units (vv%)
 598         then 'integrate (qty (uu%), qty (vv%), qty (aa%), qty (bb%)) ` (units (uu%) * units (vv%))
 599         else 'integrate (qty (uu%), qty (vv%), qty (aa% `` units (vv%)), qty (bb% `` units (vv%))) ` (units (uu%) * units (vv%)));
 600
 601 simp : true;
 602
 603 /* Rules to handle arbitrary Maxima expressions dimensionally
 604  * by replacing expressions a`b with a*gensymized(b).
 605  * When a = 0, replace it with a gensym as well, to prevent b from getting lost.
 606  */
 607
 608 matchdeclare (uu%, unitop_p);
 609 defrule (rununitify, uu%, maybe_subst_qty (uu%) * subst_units (uu%));
 610
 611 maybe_subst_qty (u) :=
 612   if qty (u) = 0 or qty (u) = 0.0 or qty (u) = 0b0
 613     then (if qtys_gensyms [units (u)] = false
 614             then (qtys_gensyms [units (u)] : gensym ("zero"),
 615                   assume (equal (qtys_gensyms [units (u)], 0))),
 616           qtys_gensyms [units (u)])
 617     else qty (u);
 618
 619 subst_units (u) := block ([v : listofvars (units (u)), substs],
 620   for v1 in v do
 621    (if units_gensyms [v1] = false then units_gensyms [v1] : gensym (),
 622     assume (units_gensyms [v1] > 0)),
 623   substs : makelist (v1 = units_gensyms [v1], v1, v),
 624   subst (substs, units (u)));
 625
 626 dimensionally (e%) ::=
 627   buildq ([e%_op : op (e%), e%_args : args (e%), context_name],
 628     block
 629      ([qtys_gensyms, units_gensyms, unitless_args, reunitification_eqs, reunitified, rezeroified],
 630       context_name : supcontext (),
 631       context : context_name,
 632       qtys_gensyms : make_array ('hashed),
 633       units_gensyms : make_array ('hashed),
 634       unitless_args : applyb1 (e%_args, rununitify),
 635       reunitification_eqs : map (lambda ([e], lhs (e) = 1 ` rhs (e)), inverse_gensym_map (units_gensyms)),
 636       reunitified : subst (reunitification_eqs, apply (e%_op, unitless_args)),
 637       rezeroified : subst (makelist (x = 0, x, listarray (qtys_gensyms)), reunitified),
 638       killcontext (context_name),
 639       rezeroified));
 640
 641 gensym_map (h) := map ("=", rest (arrayinfo (h), 2), listarray (h));
 642 inverse_gensym_map (h) := map ("=", listarray (h), rest (arrayinfo (h), 2));
 643
 644 /* Nondimensional to pure number conversion for trig functions, exp, and log.
 645  * Conversion is applied only for explicit "`" expressions;
 646  * previous implementation in terms of dimensions(e) could be time-consuming.
 647  */
 648
 649 simp : false;
 650 declare_unit_conversion (degree = %pi/180);
 651 matchdeclare (uu, lambda ([e], unitop_p(e) and dimensions(units(e)) = 1));
 652 map (lambda ([f], apply (tellsimp, [apply (f, [uu]), apply (f, [uu `` 1])])), '[sin, cos, tan, csc, sec, cot, atan]);
 653
 654 tellsimp (%e^uu, %e^(uu `` 1));
 655 tellsimp (log (uu), log (uu `` 1));
 656
 657 /* Convert atan2 arguments to pure numbers when arguments have same dimensions.
 658  * Determining dimensions are the same is potentially time-consuming.
 659  */
 660
 661 matchdeclare (aa, unitp, bb, lambda ([e], dimensions (e) = dimensions (aa)));
 662 tellsimp (atan2 (aa, bb), atan2 (qty (aa), qty (bb `` units (aa))));
 663
 664 /* Subscript applied to "`" expression => exchange "`" with subscript
 665  * Note that this applies only to literal "`" expressions.
 666  */
 667
 668 matchdeclare ([aa, bb, cc], all);
 669 tellsimpafter ((aa ` bb) [cc], aa [cc] ` bb);
 670
 671 /* "@" applied to "`" expression => exchange "`" with "@"
 672  * Note that this applies only to literal "`" expressions.
 673  */
 674
 675 matchdeclare ([aa, bb, cc], all);
 676 /* substitute into "@" expression because "@" quotes its second argument */
 677 tellsimpafter ((aa ` bb) @ cc, buildq ([aa, bb, cc], (aa @ cc) ` bb));
 678
 679 /* realpart/imagpart applied to "`" expression => exchange "`" with realpart/imagpart
 680  * Note that this applies only to literal "`" expressions.
 681  * Note further that this applies only to realpart/imagpart nouns.
 682  */
 683
 684 matchdeclare ([aa, bb], all);
 685 tellsimpafter ('realpart (aa ` bb), 'realpart (aa) ` bb);
 686 tellsimpafter ('imagpart (aa ` bb), 'imagpart (aa) ` bb);
 687
 688 simp : true;