share/ezunits/rtestezunits.mac

   1 /* ezunits: yet another units package for Maxima
   2  * This program copyright 2008 by Robert Dodier.
   3  * I release this program under the terms of the
   4  * GNU General Public License.
   5  */
   6
   7 (kill (all), reset (), load (ezunits), 0);
   8 0;
   9
  10 a`m^3;
  11 a`(m^3);
  12
  13 (a`m)^3;
  14 (a^3)`(m^3);
  15
  16 %pi/a`m;
  17 (%pi/a)`m;
  18
  19 %pi/(a`m);
  20 (%pi/a)`(1/m);
  21
  22 unit_kg: 1`kg;
  23 1`kg;
  24
  25 50*unit_kg;
  26 50`kg;
  27
  28 %pi*(a`m);
  29 %pi*a`m;
  30
  31 a*b`m*kg;
  32 (a*b)`(m*kg);
  33
  34 a/b`m/kg;
  35 (a/b)`(m/kg);
  36
  37 a`m + b`m;
  38 (a + b)`m;
  39
  40 a`m - b`m;
  41 (a - b)`m;
  42
  43 a`m - a`m;
  44 0`m;
  45
  46 a`m + b`kg;
  47 (a`m) + (b`kg);
  48
  49 (a`m)*(b`m);
  50 (a*b)`(m^2);
  51
  52 (a`m)/(b`kg);
  53 (a/b)`(m/kg);
  54
  55 (a`m)*(b`m)*(c`m);
  56 (a*b*c)`(m^3);
  57
  58 (a`m)*exp(sin(z))*(b`m);
  59 (exp(sin(z))*a*b)`(m^2);
  60
  61 (a`m)/(b`s)/(c`kg);
  62 (a/(b*c))`(m/(s*kg));
  63
  64 (a`m)/((b`s)*(c`kg));
  65 (a/(b*c))`(m/(s*kg));
  66
  67 a`m + b`m + c`m;
  68 (a + b + c)`m;
  69
  70 (a`m)*(b`m)*(c`m)*(d`m);
  71 a*b*c*d`m^4;
  72
  73 (a`kg/m^2)/(b`kg/m^2);
  74 a/b;
  75
  76 (a`kg*m/s)/(b`kg*m/s);
  77 a/b;
  78
  79 1729 ` 1;
  80 1729;
  81
  82 0 ` kg/m^3;
  83 0 ` kg/m^3;
  84
  85 /* verify different cases are handled in function units */
  86
  87 (kill (x, y, z, a),
  88  declare_constvalue (x, 1000 ` W/m^2),
  89  declare (y, dimensional),
  90  declare (z, dimensional),
  91  declare_units (z, MJ/hour),
  92  0);
  93 0;
  94
  95 units(x);
  96 W/m^2;
  97
  98 units(y);
  99 'units(y);
 100
 101 units(z);
 102 MJ/hour;
 103
 104 units(123);
 105 1;
 106
 107 units(a);
 108 1;
 109
 110 [units(x[1]), units(x[2]), units(x[k])];
 111 [W/m^2, W/m^2, W/m^2];
 112
 113 units(a ` mile/minute^2);
 114 mile/minute^2;
 115
 116 units(y(1, 2, 3));
 117 'units(y);
 118
 119 units(z(a));
 120 MJ/hour;
 121
 122 units(x*y*z);
 123 MJ*W*'units(y)/(hour*m^2);
 124
 125 units(z^a);
 126 (MJ/hour)^a;
 127
 128 units(z + sin(a)`MJ/hour + (10 ` MJ)/(24 ` hour));
 129 MJ/hour;
 130
 131 units(z - 1000`MJ/hour);
 132 MJ/hour;
 133
 134 units(y = (a ` foot)/(60 ` minute));
 135 'units(y) = foot/minute;
 136
 137 units(equal(y, (a ` foot)/(60 ` minute)));
 138 equal ('units(y), foot/minute);
 139
 140 units(y . x);
 141 W*'units(y)/m^2;
 142
 143 (declare (nounify (a), nondimensional),
 144  units(a(x, y, z)));
 145 1;
 146
 147 (remove (nounify (a), nondimensional),
 148  units(a(x, y, z)));
 149 'units(a(x, y, z));
 150
 151 /* diff, integrate, and at are handled below */
 152
 153 /* TODO: verify different cases are handled in function qty */
 154
 155 /* next one has result = unknown;
 156  * substitute narrower test on "=" instead of equal
 157  * pending simplification rule for equal applied to unitful expressions
 158 is (equal (0 ` m, 0 ` kg));
 159 false;
 160  */
 161
 162 is (0 ` m = 0 ` kg);
 163 false;
 164
 165 /* begin addt'l stuff adapted from Cliff Yapp's email
 166  * to the Maxima mailing list 2005/05/21.
 167  */
 168 a`s + c;
 169 (a`s) + c;
 170
 171 a + c`s;
 172 a + (c`s);
 173
 174 a`kg/s + c;
 175 (a`(kg/s)) + c;
 176
 177 a`kg/s + c`kg/s;
 178 (a + c)`(kg/s);
 179
 180 a`kg/s - c`kg/s + d`kg/s;
 181 (a - c + d)`kg/s;
 182
 183 a`kg*m/s^2 + c;
 184 (a`(kg*m/s^2)) + c;
 185
 186 a`kg*m/s^2 + c`kg*m/s^2;
 187 (a + c)`(kg*m/s^2);
 188
 189 a`kg*m/s^2 + c`kg*m/s^2 - d`kg*m/s^2;
 190 (a + c - d)`(kg*m/s^2);
 191
 192 a`kg/s + c`kg/s + d`kg*m/s^2;
 193 (a + c)`(kg/s) + d`(kg*m/s^2);
 194
 195 a`kg/s + c`kg*m/s^2 + d`kg*m/s^2;
 196 a`kg/s + (c + d)`kg*m/s^2;
 197
 198 (N: kg*m/s^2, 0);
 199 0;
 200
 201 (J: N*m, 0);
 202 0;
 203
 204 a`kg/s + b`N + c`kg/s + d`N + e`kg/s;
 205 (a + c + e)`kg/s + (b + d)`kg*m/s^2;
 206
 207 1/60`kg/s + a`kg/s + b`kg/s + c`kg/s^2 + d`kg/s^2;
 208 (1/60 + a + b)`kg/s + (c + d)`kg/s^2;
 209
 210 a/60`kg/s + a`kg/s + b`kg/s + c`kg/s^2 + d`kg/s^2;
 211 (a/60 + a + b)`kg/s + (c + d)`kg/s^2;
 212
 213 (a/60)`kg/s + a*b`kg/s;
 214 (a/60 + a*b)`kg/s;
 215
 216 a`kg*m/s^2 + b`N*J + c`kg/s + d`N*J - e`kg*m/s^2;
 217 (a - e)`kg*m/s^2 + c`kg/s + (b + d)`kg^2*m^3/s^4;
 218
 219 sin(1`s);
 220 sin(1`s);
 221
 222 -1`1/s;
 223 -1`1/s;
 224
 225 sin(1`kg/s);
 226 sin(1`kg/s);
 227 /* end stuff adapted from Cliff Yapp's email
 228  */
 229
 230 (remvalue (J, N), 0);
 231 0;
 232
 233 /* foo(a`s) fails nondimensional_not1 test, due to presence of a`s within.
 234  * => No rule triggered here.
 235  */
 236 (a`b)*(c`d)*foo(e`f);
 237 ''(block ([simp: false], (a`b)*(c`d)*foo(e`f)));
 238
 239 (F: C*(5/9) + 32, 0);
 240 0;
 241
 242 212`F;
 243 100`C;
 244
 245 32`F;
 246 0`C;
 247
 248 -40`F;
 249 -40`C;
 250
 251 /* deg and rad are units ad hoc; neither one is defined */
 252
 253 deg: rad*%pi/180;
 254 rad*%pi/180;
 255
 256 90`deg;
 257 %pi/2`rad;
 258
 259 a`%pi;
 260 a*%pi;
 261
 262 sin(135`deg);
 263 sin(3*%pi/4`rad);
 264
 265 sin(135`deg), rad=1;
 266 sqrt(2)/2;
 267
 268 (declare ([C1, C2], constant), 0);
 269 0;
 270
 271 /* Constants for Planck's law */
 272 [declare_constvalue (%C2, 14390.0), declare_constvalue (%C1, 3.742E+8)];
 273 [14390.0, 3.742E+8];
 274
 275 (C1: %C1`W*micron^4/m^2, C2: %C2`micron*K, l: %l`micron, T: %T`K, 0);
 276 0;
 277
 278 /* "infeval=true" is a work-around for simplification strangeness:
 279  * after load("units.mac")$ kill(all)$ load("units.mac")$
 280  * the simplifier doesn't try so hard to use tellsimpafter rules,
 281  * but it can be coaxed by re-evaluation via ''expr or infeval=true.
 282  */
 283 expr: C1/(l^5*(exp(C2/(l*T)) - 1)), infeval=true;
 284 %C1/(%l^5*(%e^(%C2/(%l*%T))-1)) ` W/(m^2*micron);
 285
 286 ev (expr, constvalue);
 287 3.742E+8/(%l^5*(%e^(14390.0/(%l*%T))-1)) ` W/(m^2*micron);
 288
 289 qty(%), %T=6000;
 290 3.742E+8 / (%l^5 * (%e^(2.398333333333333/%l) - 1));
 291
 292 /* Now try this -- should yield a skewed bump.
 293  * plot2d (%, [%l, 0.05, 3]);
 294  */
 295
 296 declare ([a1, b1, c1], dimensional);
 297 done;
 298
 299 a1 + b1 + c1;
 300 a1 + b1 + c1;
 301
 302 expand_dimensional (a1 + b1 + c1);
 303 'qty(a1)`'units(a1) + 'qty(b1)`'units(b1) + 'qty(c1)`'units(c1);
 304
 305 declare_units ('[a1, b1, c1], mph);
 306 mph;
 307
 308 a1: 50`mph;
 309 50`mph;
 310
 311 catch (b1: 50`gpf);
 312 oops ('b1, 50`gpf);
 313
 314 b1: c1;
 315 c1;
 316
 317 a1 + b1 + c1;
 318 2*c1 + 50`mph;
 319
 320 (declare_qty ('a1, 100), declare_qty ('b1, 200), declare_qty ('c1, 300));
 321 300;
 322
 323 a1 + b1 + c1;
 324 2*c1 + 50 ` mph;
 325
 326 expand_dimensional (a1 + b1 + c1);
 327 650 ` mph;
 328
 329 c1 + 1000`mph;
 330 c1 + 1000`mph;
 331
 332 expand_dimensional (c1 + 1000`mph);
 333 1300`mph;
 334
 335 expand_dimensional (c1 + 1000`kg);
 336 300`mph + 1000`kg;
 337
 338 expand_dimensional ((a1 + 500`mph)/(b1 - 100`mph));
 339 11/4;
 340
 341 map (declare_units, [aa, bb, cc, dd, ee, ff, gg], [m, kg, s, K, feet, in, acre]);
 342 [m, kg, s, K, feet, in, acre];
 343
 344 catch (aa: bb);
 345 oops ('aa, bb);
 346
 347 catch (bb: aa);
 348 oops ('bb, aa);
 349
 350 catch (cc: a1);
 351 ''(oops ('cc, a1));
 352
 353 catch (apply (lambda ([ff], cc^2), [gg]));
 354 oops ('ff, gg);
 355
 356 units (aa*bb*cc*dd*ee*ff*gg);
 357 m*kg*s*K*feet*in*acre;
 358
 359 [aa: %aa`m, bb: %bb`kg, cc: %cc`s, dd: %dd`K, ee: %ee`feet, ff: %ff`in, gg: %gg`acre];
 360 [%aa`m, %bb`kg, %cc`s, %dd`K, %ee`feet, %ff`in, %gg`acre];
 361
 362 ff: ee*(12`in)/(1`feet);
 363 12*%ee`in;
 364
 365 gg: ff^2*((1`yard)/(36`in))^2*(1`acre)/(55*88`yard^2);
 366 %ee^2/43560`acre;
 367
 368 declare_units (foo, m/s^2);
 369 m/s^2;
 370
 371 units (foo (p, q));
 372 m/s^2;
 373
 374 qty (foo (p, q));
 375 'qty (foo (p, q));
 376
 377 aa: expand_dimensional (foo (p, q) * cc * (100`s));
 378 100 * qty (foo (p, q)) * %cc ` m;
 379
 380 declare_units (bar, kg/acre);
 381 kg/acre;
 382
 383 bb: bar[n] * gg;
 384 (%ee^2/43560 ` acre) * bar[n];
 385
 386 expand_dimensional (bb);
 387 %ee^2 * 'qty (bar [n]) / 43560 ` kg;
 388
 389 catch (cc: bar[k]);
 390 oops ('cc, bar[k]);
 391
 392 cc: bar[11] * (57 ` s/kg) * (10 ` acre);
 393 bar[11] * (570 ` acre*s/kg);
 394
 395 expand_dimensional (cc);
 396 570 * qty(bar[11]) ` s;
 397
 398 kill (aa, bb, cc, dd, ee, ff, gg);
 399 done;
 400
 401 /* Unit conversions */
 402
 403 0 ` kg/m^2 `` lbm/acre;
 404 0 ` lbm/acre;
 405
 406 0 ` degC `` degF;
 407 32 ` degF;
 408
 409 100 ` degC `` degF;
 410 212 ` degF;
 411
 412 -40 ` degF `` degC;
 413 -40 ` degC;
 414
 415 1 ` acre*ft `` m^3;
 416 60228605349/48828125 ` m^3;
 417
 418 declare_unit_conversion (h = hour, d = day);
 419 done;
 420
 421 /* Stefan-Boltzmann constant */
 422 56697/10000/10^8 ` W/m^2/K^4 `` Btu/h/ft^2/R^4;
 423 304821971/178031250000000000 ` Btu/h/ft^2/R^4;
 424
 425 declare_unit_conversion (MMBtu = 10^6*Btu, kW = 1000*W, kWh = kW*h, MWh = 1000*kWh);
 426 done;
 427
 428 1 ` MMBtu `` MWh;
 429 211/720 ` MWh;
 430
 431 1 ` MMBtu/d `` kW;
 432 5275/432 ` kW;
 433
 434 (load (physical_constants), 0);
 435 0;
 436
 437 map (lambda ([x], featurep (x, physical_constant)), [%c, %mu_0, %e_0, %Z_0, %G, %h, %h_bar]);
 438 [true, true, true, true, true, true, true];
 439
 440 map (constantp, [%c, %mu_0, %e_0, %Z_0, %G, %h, %h_bar]);
 441 [true, true, true, true, true, true, true];
 442
 443 map (units, [%c, %mu_0, %e_0, %Z_0, %G, %h, %h_bar]);
 444 ''([m/s, N/A^2, s^2*A^2/(m^2*N), m*N/(s*A^2), m^3/(kg*s^2), J*s, J*s]);
 445
 446 [dimensions (1), dimensions (1729), dimensions (%e)];
 447 [1, 1, 1];
 448
 449 [dimensions (0 ` m/s), dimensions (0 ` kg/m^3), dimensions (0 ` W/kg)];
 450 [length/time, mass/length^3, length^2/time^3];
 451
 452 dimensions (foobar);
 453 'dimensions (foobar);
 454
 455 dimensions_as_list (foobar);
 456 'dimensions_as_list ('dimensions (foobar));
 457
 458 fundamental_units (foobar);
 459 'fundamental_units (foobar);
 460
 461 map (dimensions, [%c, %mu_0, %e_0, %Z_0, %G, %h, %h_bar]);
 462 [length/time,
 463  length*mass/(current^2*time^2),
 464  current^2*time^4/(length^3*mass),
 465  length^2*mass/(current^2*time^3),
 466  length^3/(mass*time^2),
 467  length^2*mass/time,
 468  length^2*mass/time];
 469
 470 fundamental_dimensions;
 471 [length, mass, time, current, temperature, quantity, luminous_intensity];
 472
 473 map (dimensions_as_list, [%c, %mu_0, %e_0, %Z_0, %G, %h, %h_bar]);
 474 [[1, 0, - 1, 0, 0, 0, 0],
 475  [1, 1, - 2, - 2, 0, 0, 0],
 476  [- 3, - 1, 4, 2, 0, 0, 0],
 477  [2, 1, - 3, - 2, 0, 0, 0],
 478  [3, - 1, - 2, 0, 0, 0, 0],
 479  [2, 1, - 1, 0, 0, 0, 0],
 480  [2, 1, - 1, 0, 0, 0, 0]];
 481
 482 map (fundamental_units, [%c, %mu_0, %e_0, %Z_0, %G, %h, %h_bar]);
 483 [m/s,
 484  kg*m/(s^2*A^2),
 485  s^4*A^2/(kg*m^3),
 486  kg*m^2/(s^3*A^2),
 487  m^3/(kg*s^2),
 488  kg*m^2/s,
 489  kg*m^2/s];
 490
 491 /* examples from message from ed romana to mailing list 2007-12-02, -12-03 */
 492
 493 e : 1/sqrt (%mu_0 * %e_0);
 494 1/sqrt (%mu_0 * %e_0);
 495
 496 [e, %c], constvalue;
 497 [299792458 ` m/s, 299792458 ` m/s];
 498
 499 /* Helium sound velocity at 27 C, He(Cp/Cv)=5/3, 4gr/mol */
 500 e : sqrt (5/3 * %R/(4 ` g/mol) * (27 ` degC `` K));
 501 sqrt(%R*(2001/16 ` mol*K/g));
 502
 503 expand_dimensional (e);
 504 sqrt(2001)*sqrt('qty(%R))/4 ` sqrt(J)/sqrt(g);
 505
 506 dimensions (e);
 507 length / time;
 508
 509 e : constvalue (e);
 510 sqrt(2079657309)/(125*2^(7/2)) ` sqrt(J)/sqrt(g);
 511
 512 dimensions (e);
 513 length / time;
 514
 515 u : fundamental_units (e);
 516 m/s;
 517
 518 e : e `` u;
 519 10^(3/2)*sqrt(2079657309)/(125*2^(7/2))  ` m/s;
 520
 521 e, numer;
 522 1019.719890214955 ` m/s;
 523
 524 /* decomposing units into multiple units of different "size" */
 525
 526 xx : 1000 ` m `` [foot];
 527 [1250000/381 ` foot];
 528
 529 xx : map (lambda ([x], x `` m), xx);
 530 [1000 ` m];
 531
 532 xx : 1000 ` m `` [foot, inch];
 533 [3280 ` foot, 1280/127 ` inch];
 534
 535 xx : map (lambda ([x], x `` m), xx);
 536 [124968/125 ` m, 32/125 ` m];
 537
 538 apply ("+", xx);
 539 1000 ` m;
 540
 541 (declare_unit_conversion (hairs_breadth = 1/100 * inch),
 542  xx : 1000 ` m `` [foot, inch, hairs_breadth]);
 543 [3280 ` foot, 10 ` inch, 1000/127 ` hairs_breadth];
 544
 545 xx : map (lambda ([x], x `` m), xx);
 546 [124968/125 ` m, 127/500 ` m, 1/500 ` m];
 547
 548 apply ("+", xx);
 549 1000 ` m;
 550
 551 xx : 10^9 ` s `` [year, month, day, hour, minute, s];
 552 [31 ` year, 8 ` month, 7 ` day, 19 ` hour, 46 ` minute, 40 ` s];
 553
 554 xx : map (lambda ([x], x `` s), xx);
 555 [978285600 ` s, 21038400 ` s, 604800 ` s, 68400 ` s, 2760 ` s, 40 ` s];
 556
 557 apply ("+", xx);
 558 1000000000 ` s;
 559
 560 xx : 10^7 ` m^2 `` [mile^2, acre, yard^2, feet^2, inch^2];
 561 [3 ` mile^2, 551 ` acre, 260 ` yard^2, 4 ` feet^2, 388096/16129 ` inch^2];
 562
 563 xx : map (lambda ([x], x `` m^2), xx);
 564 [121405692672/15625 ` m^2, 174204522558/78125 ` m^2, 16983837/78125 ` m^2, 145161/390625 ` m^2, 6064/390625 ` m^2];
 565
 566 apply ("+", xx);
 567 10000000 ` m^2;
 568
 569 86400 ` s `` [hour, minute, s];
 570 [24 ` hour, 0 ` minute, 0 ` s];
 571
 572 60 ` s `` [hour, minute, s];
 573 [0 ` hour, 1 ` minute, 0 ` s];
 574
 575 1 ` s `` [hour, minute, s];
 576 [0 ` hour, 0 ` minute, 1 ` s];
 577
 578 86401 ` s `` [day, hour, minute, s];
 579 [1 ` day, 0 ` hour, 0 ` minute, 1 ` s];
 580
 581 xx : n ` m `` [yard, foot, inch];
 582 [floor(1250*n/1143) ` yard,
 583  floor(3*(1250*n/1143-floor(1250*n/1143))) ` foot,
 584  12*(3*(1250*n/1143-floor(1250*n/1143)) -floor(3*(1250*n/1143-floor(1250*n/1143)))) ` inch];
 585
 586 xx : ev (xx, n = 100);
 587 [109 ` yard, 1 ` foot, 128/127 ` inch];
 588
 589 xx : map (lambda ([x], x `` m), xx);
 590 [124587/1250 ` m, 381/1250 ` m, 16/625 ` m];
 591
 592 apply ("+", xx);
 593 100 ` m;
 594
 595 declare_units (xyz, kg);
 596 kg;
 597
 598 xx : xyz `` [metric_ton, kg, mg];
 599 xyz `` [metric_ton, kg, mg];
 600
 601 expand_dimensional (xx);
 602 [floor('qty(xyz)/1000) ` metric_ton,
 603  floor(1000*('qty(xyz)/1000-floor('qty(xyz)/1000))) ` kg,
 604  1000000*(1000*('qty(xyz)/1000-floor('qty(xyz)/1000))
 605   - floor(1000*('qty(xyz)/1000-floor('qty(xyz)/1000)))) ` mg];
 606
 607 ev (xx, xyz = (1729 + 42/99) ` kg, nouns);
 608 [1 ` metric_ton, 729 ` kg, 14000000/33 ` mg];
 609
 610 /* needs more work to get rule for diff to apply to verb ... oh well. */
 611 'diff (x (t) ` inch, t ` year, 1);
 612 'diff (x (t), t, 1) ` inch / year;
 613
 614 'diff (x (t) ` inch, t ` year, 2);
 615 'diff (x (t), t, 2) ` inch / year^2;
 616
 617 'integrate (f (x) ` N, x ` m, 1 ` m, 10 ` m);
 618 ('integrate (f (x), x, 1, 10)) ` N*m;
 619
 620 'integrate (f (x) ` N, x ` m, 1 ` inch, 10 ` inch);
 621 ('integrate (f (x), x, 127/5000, 127/500)) ` N*m;
 622
 623 /* small example from mma mailing list circa 2009-05-29 */
 624
 625 block ([R, w, L, x, y, z], kill (x), R : x ` Ohm, w : z ` Hz, L : y ` H, (R + w*L) `` Ohm);
 626 (x + y*z) ` Ohm;
 627
 628 /* test new definitions of pound_mass, gallon, slug */
 629
 630 ((1 ` pound_force)/(1 ` pound_mass), [dimensions (%%), %% `` cm/s^2]);
 631 [length/time^2, ''(980665/1000) ` cm/s^2];
 632
 633 (kill (xx), xx ` gallon `` inch^3);
 634 231*xx ` inch^3;
 635
 636 xx ` gallon/acre `` feet;
 637 ''(231*xx/(55*88*3^2)/(12^2)/12) ` feet;
 638
 639 [dimensions (slug), dimensions (slug/pound_mass)];
 640 [mass, 1];
 641
 642 1 ` slug/pound_mass `` 1;
 643 196133/6096; /* approx 32.something */
 644
 645 /* some examples suggested by Jeff Hankin; thanks, Jeff. */
 646
 647 dimensions (W / (K * m));
 648 length*mass/(temperature*time^3);
 649
 650 dimensions (lux);
 651 luminous_intensity / length^2;
 652
 653 dimensions (parsec);
 654 length;
 655
 656 bfloat (1 ` pc `` m);
 657 3.085677581491367b16 ` m;
 658
 659 dimensions (kat);
 660 quantity / time;
 661
 662 bfloat (300 ` mile `` AU);
 663 3.227340053309997b-6 ` AU;
 664
 665 10^6 ` pc `` AU;
 666 648000000000/%pi ` AU;
 667
 668 70 ` mile/hour `` parsec / julian_year;
 669 4762373*%pi/467493345937500 ` parsec/julian_year;
 670
 671 /* still more examples */
 672
 673 /* assignment to a unit symbol causes trouble;
 674  * maybe rhs of "`" should never be evaluated.
 675  * For now, just remvalue all unit symbols.
 676  */
 677 (apply (remvalue, known_units ()), 0);
 678 0;
 679
 680 (kill (xx), xx : 1 ` ml*lux*katal/(gray*Bq), dimensions (xx));
 681 luminous_intensity * quantity * time^2 / length;
 682
 683 fundamental_units (xx);
 684 cd * mol * s^2 / m;
 685
 686 (kill (uu), uu ` horsepower*julian_year/tbsp `` GJ/liter);
 687 50133999444908841*uu/31501953125000 ` GJ/liter;
 688
 689 'f(uu) ` MW `` short_ton*(lbf/lbm)*mile/hour `` Btu/hour;
 690 720000000*'f(uu)/211 ` Btu/hour;
 691
 692 1 ` W `` Btu/julian_year `` newton*foot/hour `` lbf*m/s `` W;
 693 1 ` W;
 694
 695 /* units ad hoc; no fundamental units or dimensions, just conversions */
 696
 697 declare_unit_conversion (1`elephant=10`cow, 1`cow=10`goat, 1`goat = 2`dog, 1`dog=3`cat);
 698 done;
 699
 700 1 ` elephant/m `` cat/inch;
 701 381/25 ` cat/inch;
 702
 703 dimensions (1 ` goat/hour);
 704 'dimensions(cat) / time;
 705
 706 /* convert argument of trig function */
 707
 708 sin (45 ` degree);
 709 1/sqrt(2);
 710
 711 cos (180 ` degree);
 712 -1;
 713
 714 tan (135 ` degree);
 715 -1;
 716
 717 sec (-45 ` degree);
 718 sqrt(2);
 719
 720 csc (30 ` degree);
 721 2;
 722
 723 cot (60 ` degree);
 724 1/sqrt(3);
 725
 726 sin (30 ` kg);
 727 sin (30 ` kg);
 728
 729 cos (100 ` m);
 730 cos (100 ` m);
 731
 732 is (abs (ev (sin (40 ` degree), numer) - 0.642787609686539) < 1e-12);
 733 true;
 734
 735 is (abs (ev (cot (80 ` degree), numer) - 0.176326980708465) < 1e-12);
 736 true;
 737 /* Hydrostatics examples from:
 738  * http://www.ce.siue.edu/mrossow/Worked_examples_Internet_text-only/Data_files-html/Hydrostatics_list_of_examples.html
 739  * via: http://designerunits.com/examples/sites/ce.siue/hydrostatics
 740  */
 741
 742 (gee :(980 + 2/3) ` cm/s^2, p_AB : (10^3`kg/m^3) * gee * (5`m) `` kPa);
 743 1471/30 ` kPa;
 744
 745 w : p_AB * 4 ` m `` kN/m;
 746 2942/15 ` kN/m;
 747
 748 (kill (F_A), e : (F_A`kN) * (3`m) - w*(3`m)*(3`m)/2, solve (first (e), F_A));
 749 [F_A = 1471/5];
 750
 751 p : ((200 ` N)/(%pi*r^2) + (900 ` kg/m^3) * (200 ` mm) * %pi * r^2 * gee / (%pi * r^2)) `` kPa, r = (25 ` mm)/2;
 752 (1280/%pi + 4413/2500) ` kPa;
 753
 754 ratsimp (W : (p * %pi * r^2) `` kN), r = (300 ` mm)/2;
 755 (39717*%pi + 28800000)/1000000 ` kN;
 756
 757 ratsimp (W * 1/(200 ` N) `` 1);
 758 (39717*%pi + 28800000)/200000;
 759
 760 kill (W);
 761 done;
 762
 763 p_B : (10^3 ` kg/m^3) * gee * (2 ` m + 3 ` m) `` kPa;
 764 1471/30 ` kPa;
 765
 766 w_B : p_B * 4 ` m `` kN/m;
 767 2942/15 ` kN/m;
 768
 769 1/2 * w_B * (2`m + 3`m);
 770 1471/3 ` kN;
 771
 772 1/3*(2 ` m + 3 ` m);
 773 5/3 ` m;
 774
 775
 776 /* Miscellaneous examples from: "Metric Units and Conversion Charts"
 777  * by Theodore Wildi (www.wildi-theo.com) ISBN 0-7803-1050-0
 778  * via: http://designerunits.com/examples/sites/wildi-theo/chartschapter04
 779  */
 780
 781 (400 ` lbm) * (3 ` cm/s^2);
 782 1200 ` lbm*cm/s^2;
 783
 784 (400 ` lbm) * (3 ` cm/s^2) `` N;
 785 136077711/25000000 ` N;
 786
 787 (declare_unit_conversion (ozf = lbf/16), (6 ` ozf)/(3 ` kg) `` m/s^2);
 788 8896443230521/16000000000000 ` m/s^2;
 789
 790 (4 ` mile) * (3 ` feet) * (2 ` mm) * (20 ` kg/foot^3) `` kg;
 791 1056000/127 ` kg;
 792
 793 (150 ` lbm) * (2*%pi/(1`day))^2 * (4000 ` mile) `` N;
 794 63366854089*%pi^2/270000000000 ` N;
 795
 796 (if get ('physical_constants, 'version) = false then load (physical_constants),
 797  declare_unit_conversion (1 ` eV = (1 ` V) * constvalue (%%e) `` J),
 798  vv : solve (Ek = 1/2 * m * v^2, v),
 799  v1 : ev (second (vv), Ek = 1`eV, m = 1`lbm));
 800 v = sqrt(2) ` sqrt(eV)/sqrt(lbm);
 801
 802 map (dimensions, v1);
 803 'dimensions(v) = length/time;
 804
 805 map (lambda ([x], x `` cm/year), v1);
 806 v `` cm/year = 39447*sqrt(1602176487)/(15625*2^(5/2)*sqrt(45359237)) ` cm/year;
 807
 808 (declare_unit_conversion (gauss = tesla/10^4, G = gauss, kG = 1000*G, mph = mile/hour), (10 ` kG) * (6 ` inch) * (60 ` mph) `` V);
 809 1596771/390625 ` V;
 810
 811 (declare_unit_conversion (1 ` cmil = %pi * ((1/1000 ` inch) / 2)^2, microohm = Ohm/10^6),
 812  (172/100 ` microohm*cm) * (8/10 ` mile) / (10400 ` cmil) `` Ohm);
 813 544896/(41275*%pi) ` Ohm;
 814
 815 (declare_unit_conversion (cal_tc = 4184/1000 * J, Btu_tc = 105435026444/10^8 * J, ton = 2000 * lbm),
 816  (3 ` ton) * (1000 ` degF `` R - 70 ` degF `` R) * (9/100 ` cal_tc/g/K),
 817  En : %% `` cal_tc);
 818 12655227123/100 ` cal_tc;
 819
 820 En `` MJ;
 821 6618683785329/12500000000 ` MJ;
 822
 823 En `` Btu_tc;
 824 1203397051878000/2396250601 ` Btu_tc;
 825
 826 ((3/10 ` Btu*inch/hour/foot^2/R) * (150 ` m^2) * ((20 ` degC `` K) - (-15 ` degC `` K)) / (8 ` cm),
 827  loss : %% `` Btu/hour);
 828 4921875/508 ` Btu/hour;
 829
 830 loss `` kW;
 831 184625/65024 ` kW;
 832
 833 (kill (p, V, n, R, T, foo),
 834  solve ([p*V = n*R*T, n = m/M], [p, n]),
 835  ev (first (first (%%)), R=constvalue (%R), m=5`oz, M=40`g/mol, V=678`inch^3, T=72`degF``K),
 836  foo : rhs (%%) `` Pa);
 837 2506412710065636911/3199803664896 ` Pa;
 838
 839 (declare_unit_conversion (atm = 101325/1000*kPa),
 840  foo `` atm);
 841 358058958580805273/46317158049369600 ` atm;
 842
 843 (ev (nu * rho * d / eta, nu = 5 ` mph, rho = 1 ` kg/liter, d = 7/8 ` inch, eta = 8*10^-6 ` lbf*s/foot^2),
 844  Re : %% `` 1);
 845 14984331321600000/115538223773;
 846
 847 620 ` degC `` K;
 848 17863/20 ` K;
 849
 850 23 ` degF `` K;
 851 5363/20 ` K;
 852
 853 23 ` degF `` degC;
 854 -5 ` degC;
 855
 856 (87 ` degF `` R - 60 ` degF `` R) `` K;
 857 15 ` K;
 858
 859 /* from: http://www.kc9aop.net/Doc/link_pages/motors_and_mechanical.htm
 860  * via: http://designerunits.com/examples/sites/kc9aop/motors_and_mechanical
 861  */
 862
 863 (declare_unit_conversion ('Hz[S] = 2*%pi*Hz, rpm = 2*%pi/(60*s)),
 864  omega_mech : ev (2/npoles * 60 ` 'Hz[S] `` rpm, npoles = 2));
 865 3600 ` rpm;
 866
 867 (40 ` hp)/(1725 ` rpm) `` lbf*foot;
 868 8800/(23*%pi) ` foot*lbf;
 869
 870 compute_conversion_factor (hp/rpm, lbf*foot);
 871 16500/%pi;
 872
 873 (declare_unit_conversion (radian=1),
 874  (70 ` lbf)*(3 ` feet)/(1 ` radian) `` lbf*foot);
 875 210 ` lbf*foot;
 876
 877 (30 ` hp)/(1725 ` rpm) `` lbf*foot;
 878 6600/(23*%pi) ` foot*lbf;
 879
 880 (230 ` V)*(4 ` A)*82/100 `` hp;
 881 37720000000000000/37284993579113511 ` hp;
 882
 883 (1725 ` rpm)*(31/10 ` lbf*foot) `` hp;
 884 713*%pi/2200 ` hp;
 885
 886 omega_mech : ev (2/npoles * 50 ` 'Hz[S] `` rpm, npoles = 4);
 887 1500 ` rpm;
 888
 889 /* from: http://www.ibiblio.org/harris/500milemail.html
 890  * "The case of the 500-mile email" -- a pleasant diversion
 891  */
 892
 893 declare_unit_conversion (1 ` lightsecond = constvalue (%c) * (1 ` s), millilightsecond = lightsecond/1000);
 894 done;
 895
 896 3 ` millilightsecond `` mile;
 897 149896229/268224 ` mile;
 898
 899 block ([fpprec : 16], 3b0 ` millilightsecond `` mile);
 900 5.588471911536626b2 ` mile;
 901
 902 /* from the mailing list 2015-03-18: "Some ezunits questions" */
 903
 904 (1000 ` C+66.0*10^-9 ` F*V)*(46*10^3 ` Hz*V)``W;
 905 4.6000000003036E+7 ` W;
 906
 907 /* subscripts applied to unitful expressions */
 908
 909 (kill (x, k), (x ` kg/m^2)[k]);
 910 x[k] ` kg/m^2;
 911
 912 ([1, 2, 3] ` m/s^2)[3];
 913 3 ` m/s^2;
 914
 915 (kill (foo, i), declare_units (foo, MJ/cd), foo[i]);
 916 /* verify that foo[i] not expanded to qty(foo[i]) ` MJ/cd or anything else */
 917 foo[i];
 918
 919 /* "@" applied to unitful expressions */
 920
 921 (kill (foo, a, b, c), defstruct (foo (a, b, c)), 0);
 922 0;
 923
 924 (block ([x : new (foo ('x1, %pi, 17)) ` m/s], [x@a, x@b, x@c], ev (%%)));
 925 [x1 ` m/s, %pi ` m/s, 17 ` m/s];
 926
 927 /* realpart/imagpart applied to unitful expressions
 928  * ev(..., nouns) seems to be necessary to work around noun/verb strangeness ... sigh
 929  */
 930
 931 (kill (x), x : (17 - %e*%i) ` s*A/kg, ['realpart (x), 'imagpart (x)], ev (%%, nouns));
 932 [17 ` s*A/kg, -%e ` s*A/kg];
 933
 934 /* examples derived from dissertation of Andrew John Kennedy */
 935
 936 /* statistics of a list of numbers x:
 937  * mean has same units as x
 938  * variance has units = (units(x))^2
 939  * std deviation has same units as x
 940  * skewness is dimensionless
 941  * covariance has units = units(x) * units(y)
 942  * correlation is dimensionless
 943  */
 944
 945 (kill (x, n), x : [2 ` kg/m^2, 3 ` kg/m^2, 5 ` kg/m^2, 7 ` kg/m^2, 11 ` kg/m^2], n : length (x), 0);
 946 0;
 947
 948 (kill (mean_x), mean_x : apply ("+", x) / n);
 949 28/5 ` kg/m^2;
 950
 951 (kill (mean_x2, var_x), mean_x2 : apply ("+", map (lambda ([x1], x1^2), x)) / n, var_x : mean_x2 - mean_x^2);
 952 256/25 ` kg^2/m^4;
 953
 954 (kill (std_x), std_x : sqrt (var_x));
 955 16/5 ` kg/m^2;
 956
 957 (kill (mean_x3, skewness_x),
 958   mean_x3 : apply ("+", map (lambda ([x1], x1^3), x)) / n,
 959   moment3_x : mean_x3 - 3 * mean_x2 * mean_x + 2 * mean_x^3,
 960   skewness_x : moment3_x / std_x^3);
 961 1197/2048;
 962
 963 (kill (y), y : [13 ` s*A, 17 ` s*A, 19 ` s*A, 23 ` s*A, 29 ` s*A], y_length : length (y), is (y_length = n));
 964 true;
 965
 966 (kill (mean_y), mean_y : apply ("+", y) / y_length);
 967 101/5 ` s*A;
 968
 969 (kill (mean_y2, var_y, std_y), mean_y2 : sum (y[i]^2, i, 1, y_length) / y_length, var_y : mean_y2 - mean_y^2, std_y : sqrt (var_y));
 970 2*sqrt(186)/5 ` s*A;
 971
 972 (kill (mean_xy, covar_xy), mean_xy : sum (x[i] * y[i], i, 1, n) / n, covar_xy : mean_xy - mean_x * mean_y);
 973 432/25 ` kg*s*A/m^2;
 974
 975 (kill (corr_xy), corr_xy : covar_xy / std_x / std_y);
 976 27/(2*sqrt(186));
 977
 978 /* differentiation via secant approximation:
 979  * diff(f(x), x) = (f(x + h) - f(x - h))/(2 h)
 980  *  => diff has units = units(f) / units(h) = units(f) / units(x)
 981  *                                            since units(x) = units(h)
 982  */
 983
 984 (kill (f, x, y, h, e, e1),
 985  declare ([f, x, y, h], dimensional),
 986  declare_units (h, units(x)),
 987  units ('diff (f(x), x)));
 988 units(f)/units(x);
 989
 990 units ((f(x + h) - f(x - h))/(2*h));
 991 units(f)/units(x);
 992
 993 /* additional examples for 'diff(...) */
 994
 995 units ('diff (f(x, y), x, 1, y, 1));
 996 units(f)/(units(x)*units(y));
 997
 998 units ('diff (f(x, h), x, 1, h, 1));
 999 units(f)/units(x)^2;
1000
1001 units ('diff (f(x, y), x, 2, y, 4));
1002 units(f)/(units(x)^2*units(y)^4);
1003
1004 (kill (m, n),
1005  units ('diff (f(x, y), x, m, y, n)));
1006 units(f)/(units(x)^m*units(y)^n);
1007
1008 /* second divided difference */
1009 (declare_units (x1, units(x)),
1010  declare_units (x2, units(x)),
1011  e: (f(x + h) - f(x - h))/(2*h),
1012  e1: subst (x = x1, e),
1013  e2: subst (x = x2, e),
1014  units ((e2 - e1)/(x2 - x1)));
1015 units(f)/units(x)^2;
1016
1017 /* integration via trapezoid rule:
1018  * integral(f(x), x, a, b) = (h/2) (f(a) + 2 f(a + h) + ... + 2 f(b - h) + f(b)), h = (b - a)/n
1019  *  => integral has units = units(x) * units(f)
1020  */
1021
1022 (kill (a, b),
1023  declare ([a, b], dimensional),
1024  units ('integrate (f(x), x, a, b)));
1025 units(f)*units(x);
1026
1027 /* TODO: also express summation as a 'sum(...) expression and look at units */
1028 (e1: (h/2)*(f(a) + 2*f(a + h) + 2*f(a + 2*h) + 2*f(b - 2*h) + 2*f(b - h) + f(b)),
1029  units (e1));
1030 units(f)*units(x);
1031
1032 /* additional examples for 'integrate(...) */
1033
1034 (declare_units (f, mile/hour^2),
1035  declare_units (t, hour),
1036  units ('integrate (f(t), t, %a ` hour, %b ` hour)));
1037 mile/hour;
1038
1039 units ('integrate ('integrate (f(t), t), t));
1040 mile;
1041
1042 /* root finding via Newton-Raphson algorithm:
1043  * x[n + 1] = x[n] - f(x[n])/f'(x[n])
1044  *  => each x[n] has same units as x[1]
1045  */
1046
1047 units (x[n + 1] = x[n] - f(x[n])/'at('diff(f(x), x), x = x[n]));
1048 units(x) = units(x);
1049
1050 /* solve for units of variable in equation
1051  * see: http://groups.google.com/group/sage-support/browse_thread/thread/a60c943c224dfe87#
1052  */
1053
1054 (kill (eqn, H_l, h_c, T_a, T_l),
1055  eqn : H_l = h_c * (T_a - T_l),
1056  declare_units (H_l, calorie/cm^2/minute),
1057  declare_units ([T_a, T_l], K),
1058  declare (h_c, dimensional),
1059  eqn : units (eqn));
1060 calorie/(cm^2*minute) = 'units(h_c)*K;
1061
1062 solve (eqn, units (h_c));
1063 ['units(h_c) = calorie/(cm^2*minute*K)];
1064
1065 /* tickles bug in tellsimpafter: this example triggers an error
1066  * unless built-in simplifications are applied so that MRAT is
1067  * simplified away by the time GREAT is called.
1068  */
1069 (kill (a, b, c, x), block ([x : (rat (1) ` a) * (1 ` b + c)], expand (x)));
1070 1 ` a*b + c ` a;
1071
1072 /* "*" rule now splits arguments into 3 types:
1073  * dimensional, nondimensional, and terms not known to be either one.
1074  * Thanks to Stefan Karrmann for suggesting it.
1075  */
1076 (kill (a, b, c, d, e, f, g, h, foo, bar),
1077   exp (a ` b) * (c ` d) * (e ` f) * g * h);
1078 exp (a ` b) * (c * e * g * h ` d * f);
1079
1080 (c ` d) * (e ` f) * g * h;
1081 c * e * g * h ` d * f;
1082
1083 (c ` d) * (e ` f) * (g ` h);
1084 c * e * g ` d * f * h;
1085
1086 foo (a ` b) * bar (c ` d) * (e ` f) * (g ` h) * 2 * %pi;
1087 foo (a ` b) * bar (c ` d) * (2 * %pi * e * g ` f * h);
1088
1089 sin (a ` b) * bar (c ` d) * (e ` f) * (g ` h) * 2 * %pi;
1090 bar (c ` d) * (2 * %pi * e * g * sin (a ` b) ` f * h);
1091
1092 exp (a ` b) + c ` d + e ` d + g + h;
1093 exp (a ` b) + (c + e) ` d  + g + h;
1094
1095 c ` d + e ` d + g + h;
1096 (c + e) ` d  + g + h;
1097
1098 c ` d + e ` d;
1099 (c + e) ` d;
1100
1101 foo (a ` b) + c ` d + e ` d + 2 + %pi;
1102 (c + e) ` d + foo (a ` b) + 2 + %pi;
1103
1104 /* test rules for relational expressions */
1105
1106 (kill (a, b, x),
1107  block ([a : 17 ` m, b : x ` mm],
1108   [a < b, a <= b, a >= b, a > b,
1109    /* SEE NOTE ABOUT SIMPLIFICATION OF DIMENSIONAL "=" IN EZUNITS_FUNCTIONS !!
1110    a = b, a # b,
1111     */
1112    equal (a, b), notequal (a, b)]));
1113 [17 < x/1000, 17 <= x/1000, 17 >= x/1000, 17 > x/1000,
1114  /*
1115  17 = x/1000, 17 # x/1000,
1116   */
1117  equal (17, x/1000), notequal (17, x/1000)];
1118
1119 /* I'M NOT COMFORTABLE WITH IMPLICIT CONVERSION OF 0 TO DIMENSIONAL !!
1120 block ([a : b ` kg],
1121  [a < 0, 0 < a,
1122   a <= 0, 0 <= a,
1123   a >= 0, 0 >= a,
1124   a > 0, 0 > a,
1125   a = 0, 0 = a,
1126   a # 0, 0 # a,
1127   equal (a, 0), equal (0, a),
1128   notequal (a, 0), notequal (0, a)]);
1129 [b < 0, 0 < b,
1130  b <= 0, 0 <= b,
1131  b >= 0, 0 >= b,
1132  b > 0, 0 > b,
1133  b = 0, 0 = b,
1134  b # 0, 0 # b,
1135  equal (b, 0), equal (0, b),
1136  notequal (b, 0), notequal (0, b)];
1137  */
1138
1139 /* sin, cos, etc of dimensionless quantities */
1140
1141 (foo1 : [sin, cos, tan, csc, sec, cot, exp, log],
1142  kill (x, y),
1143  z : (x ` m)/(y ` foot) `` 1);
1144 (1250*x)/(381*y);
1145
1146 map (lambda ([f1], f1((x ` m)/(y ` foot))), foo1);
1147 ''(makelist (f1(z), f1, foo1));
1148
1149 atan2 (x ` m, y ` foot);
1150 atan2 (x, 381*y/1250);
1151
1152 exp (- (t ` hour)/(T ` s));
1153 exp (- 3600*t/T);
1154
1155 /* in this next example, %%k*(x ` K)/(y ` J) is dimensionless,
1156  * but it is a little too complex to be detected by the rule for
1157  * dimensionless argument of exp.
1158  * We'll help out a little.
1159  */
1160 foo: exp (-%%k*(x ` K)/(y ` J));
1161 exp (%%k*((-x/y) ` K/J));
1162
1163 map (lambda ([e], if dimensions(e) = 1 then qty(e) else e), foo);
1164 exp (-qty(%%k)*x/y);
1165
1166 /* fractions of a farad */
1167
1168 1`pF*1`V*1`Hz `` A;
1169 1/1000000000000 ` A$
1170
1171 1 ` nF `` C/V;
1172 1/1000000000 ` C/V;
1173
1174 1 ` microF/nF `` 1;
1175 1000;
1176
1177 /* Reported on mailing list by Antonio Lapira 2015-04-04 - FIXED */
1178 3`ha``m^2;
1179 30000`m^2;
1180
1181 /* ensure nondimensional/(sum of dimensional) `` whatever = nondimensional/((sum of dimensional) `` 1/whatever)
1182  * inspired by mailing list 2018-12-12: "ezunits: strange multiplication behaviour?"
1183  */
1184
1185 1/(sqrt(1`henry*1`farad)) `` 1/s;
1186 1 ` 1/s;
1187
1188 1/(1`s+sqrt(1`henry*1`farad)) `` 1/s;
1189 1/2 ` 1/s;
1190
1191 /* other examples */
1192
1193 (%e + %pi) / (100 ` hour - 100 ` s + x ` ms) `` Hz;
1194 (%e + %pi) / (x/1000 + 359900) ` Hz;
1195
1196 /* SF bug #3789: "package ezunits: ev(dimensions(u), nouns) stack overflow" */
1197
1198 (foo: 81*sin(acos((77-x)/81)-atan(7/20)), 0);
1199 0;
1200
1201 /* primary test here is just the absence of stack overflow error;
1202  * numerical value not really crucial
1203  */
1204 is (abs (find_root (foo - 50.58515344111066, x, 5, 55, abserr = 1e-6) - 34.0) <= 1e-6);
1205 true;
1206
1207 (kill(u), ev (dimensions(u), nouns));
1208 'dimensions(u);
1209
1210 /* now try this:
1211  * plot2d (foo, [x, 5, 55]);
1212  */