LAST REVISION: 14 Apr 2012.

By Steven Byrnes, June 2008. This is intended as a package to facilitate the use of units in Mathematica, especially when dimensionful quantities need to be subjected to numerical operations. I'm releasing this notebook to the public domain. I also plan to keep it posted on my personal website, http://sjbyrnes.com/

MOTIVATION / IDEA FOR THIS PACKAGE: First of all, Mathematica has unit-handling capabilities...more or less. With the built-in "Units" package, you can make a quantity equal to "5 Coulomb", and if you multiply it by "5 Volt" and tell Mathematica to express the answer in SI, then you will indeed get "5 Joule", which it can also convert to eV, erg, etc. The problem with this approach, though, is that a quantity equal to "5 Coulomb" is not a numerical quantity (it's equal to five times a symbolic variable called "Coulomb"), so you can't put it inside a "NIntegrate", "FindRoot", or any other numerical operations, unless you first divide out the units to make it unitless, do the operation, and then tack the units back on.

In this package, a different approach is used. Instead of storing "Coulomb" and other units as symbolic variables, they're set to be arbitrary numerical values. Then, all numerical operations can go through -- numerical integration, findroot, and whatever else -- and it is quite easy to both display answers in whatever dimensions you want, and to check that your calculation is free of dimensional errors. (See below.)

TO CREATE THIS PACKAGE: What you're reading right now is a "notebook", not a package. To make it into a package, you go to the "Save As" menu, and tell it to save as a package. In order for Mathematica to be able to find it and load it, you need to save it into a proper directory: Type $Path into a blank Mathematica notebook, and it will tell you which directories will work for this. Or, you can put it in any directory, and then modify $Path to point there, as described on this webpage: http://reference.wolfram.com/mathematica/ref/$Path.html
You can name the package whatever you want, as long as it doesn't happen to conflict with other Mathematica package names.

TO USE THIS PACKAGE IN A PROGRAM: Right at or near the start of your Mathematica program, run the command
  <<XXYYZZ`
where "XXYYZZ.m" is the name that you've given this package. Then, go ahead and do your calculation with units, e.g. writing "5 angstrom^2" for 5 squared angstroms, and so forth. To display a result in a certain unit, you divide the (numerical) answer you get by the (numerical) unit value, e.g. for a mobility you might write μInCmSqPerVoltSecond = μ / (cm^2/(volt second)).

To check that your dimensional analysis is correct, you go back to the start of the program and rerun the notebook, including the <<XXYYZZ` command. The base units will take on different numerical values (you can see below that they are defined to be randomized each time), but the answers will come out exactly the same...if your dimensions are correct. Conversely, if the dimensional analysis is not 100% consistent, it's extremely likely that you'll get a very different final answer when you re-run the notebook.

You can also do straight-up unit conversions using this package. For example, How many torr in a psi? Well, 1 psi = X torr, and that gives X=psi/torr. So to calculate how many torr in a psi, you run the command psi/torr. To put it another way, if you want a quantity (e.g. 1 psi) "expressed in torr", you think, you want this quantity expressed as a multiple of 1 torr. So you replace the phrase "expressed in torr" with "divided by (1 torr)", or "divided by torr". Again, re-running the notebook will not affect these unit-conversion calculations if and only if you're converting between two things that do actually have equivalent dimension.

EXAMPLE 1: If a 1-milligram mass free-falls 1 nm under Earth's gravity, what is its final kinetic energy? Express the answer in eV. 
In[1]:= <<XXYYZZ`
In[2]:= KineticEnergy = (1 mg) * (9.8 meter/sec^2) *(1nm);
In[3]:= KineticEnergy / eV
Out[3]:= 61166.8

EXAMPLE 2: I want to accelerate a proton from rest to 1 m/s in 1 microsecond. How large an electric field do I need? Express the answer in volts/cm.
In[1]:= <<XXYYZZ`
In[2]:= Solve[e * (X volt/cm) * (1 μs) / (1 amu) ==  1 meter/sec, X]
Out[2]:= {{X -> 0.000103643}}

EXAMPLE 3: A car has an engine with maximum power of 100 kW. When the car travels at speed v, it requires a power of K1*v (where K1=100 N) to overcome tire friction, and a power of K2*(v+5m/s)^2 (where K2=100 kg/s) to overcome wind resistance (the "5 m/s" is the wind speed), and a power of 1 kW to run the air conditioner. What is the car's maximum velocity? Express the answer in km per hour.
In[1]:= <<XXYYZZ`
In[2]:= K1 = 100 newton; K2 = 30 kg/sec;
In[3]:= RequiredEnginePower[v_] := K1 * v + K2 * (v+5 meter/sec)^2 + 1 kW;
In[4]:= Solve[RequiredEnginePower[X km/hour] == 100kW]
Out[4]:={{X -> -231.413}, {X -> 183.413}}

(In other words, if you drive "against the wind", you can go 183 km/h, and if you drive "with the wind", you can go 231 km/h.)

EXAMPLE 4: What is 1 psi expressed in torr?
In[1]:= <<XXYYZZ`
In[2]:= 1 psi / torr
Out[2]:= 51.7151

EXAMPLE 5: You measured a voltage as a function of the position of dial: 10 volts when the dial is at 1cm, .... Interpolate from this data to get the expected voltage when the dial is at 41mm, and express the answer in mV.
In[1]:= <<XXYYZZ`
In[2]:= VoltageData = {{1 cm, 10 volt}, {2 cm, 11 volt}, {3 cm, 13 volt}, {4 cm, 16 volt}, {5 cm, 18 volt}};
In[3]:= VoltageFunction = Interpolation[VoltageData];
In[4]:= VoltageFunction[41 mm] / mV
Out[4]:= 16278.
*Another way to type In[2] would be:  In[2]:= {#[[1]] cm, #[[2]] volt}& /@ {{1,10},{2,11},{3,13},{4,16},{5,18}};

TO MODIFY THIS PACKAGE: Simply add to or modify the input "cells" below. Don't add a new cell, though -- there's certain cell formatting that's required for a package to run properly, which I've already set up. One possible modification, for example, would be that if there are units or constants that you plan to use all the time, you may want to shorten their names (e.g. "h" instead of "hPlanck"). Another would be adding more physical constants and units. Another would be that if you want to free up the variable "e" for other uses in your programs, change the variable name here.

NOTE ON ELECTROMAGNETISM UNITS: The distinctions between cgs/Gaussian and SI electromagnetism units are too complicated to capture with simple conversion factors. For  example, the definition of "coulomb" versus "statcoulomb": (1 statcoulomb)^2/(1 meter) has units of energy, whereas (1 Coulomb)^2/(1 meter) does not. Consequently, even though a physical object with a charge of 1 coulomb has a charge of 2997924580 statcoulombs, you can't quite say that 1 coulomb "equals" 2997924580 statcoulombs. It's more like "corresponds to". Any electromagnetism quantity (electric field, etc.) can be defined with underlying cgs/Gaussian conventions, or with underlying SI conventions, but not both simultaneously. In this document, most units and constants are defined with underlying SI conventions; two (statcoulomb, statvolt) are defined with underlying cgs conventions; and four others (the elementary charge, 1 gauss, 1 oersted, Bohr magneton) are defined in both ("e" and "ecgs", "gaussSI" and "gausscgs", "oerstedSI" and "oerstedcgs", "μBohr" and "μBohrcgs").

NOTE ON MOLES: If you're working with moles, you can include "mol" as an extra base unit. In this case, the definition of Avogadro's number below should have a mol^-1 appended onto the end. Also, you may want to re-define "amu" to equal "gram/mol" instead of its current definition. You can add molar, molal, etc. as derived units.

WARNINGS / PITFALLS:
*This package is intended to be instead of, not in addition to, the built-in Mathematica packages "Units" and "PhysicalConstants". They may have conflicting definitions, so I doubt they'll work simultaneously.
*This was created and tested in Mathematica version 6.0. In other versions, I suspect it would work just as well, but no guarantees. If there's a problem, try deleting all the text and comments, so there's just the Mathematica commands and nothing else. Then try again to save as a package file.
*I tried to correctly copy down all the units and constants, but you should check for errors yourself and email me if you find any.
*If your program uses the command SeedRandom[n], you can't check your dimensional analysis as described here -- each time you run the program the base units will take the same random values, leaving no trace of your mistake. Instead, you can for example start your program with the commands
  SeedRandom[];
  <<XXYYZZ`
  SeedRandom[n];
or you could alter the base units by hand, etc.
*You should load the package as described here. Don't use the "Needs" command. Don't load the package anywhere except for once at the start of the notebook.
*Some numerical functions (such as FindRoot) will iterate until EITHER "PrecisionGoal" (which relates to the number of significant figures) OR "AccuracyGoal" (which relates to the absolute error) is satisfied. With the default AccuracyGoal, FindRoot doesn't care whether the answer is 1×(10^(-30)) or 5×(10^(-30)), even though it may make a very big difference to you. (This is especially liable to come up when you're using this package; e.g. "1 eV" may correspond numerically to (10^(-30)).) To fix this error, always put in the option "AccuracyGoal -> Infinity". By the way, other numerical functions, e.g. NIntegrate, have "AccuracyGoal -> Infinity" by default, so no need to worry about it.
*I used some special characters, along with how to input it. For example, to input a μ (Greek letter "mu"), press the escape key, then the letter m, then the escape key again.
*This has nothing to do with this package, but see the following webpage: http://reference.wolfram.com/mathematica/tutorial/UnconstrainedOptimizationSymbolicEvaluation.html
(2nd paragraph below the pink box and below) for a very common bug that arises in Mathematica numerical programming, and how to fix it.

SECTION 1: BASE UNITS

meter=(10^(Random[Real,{-3,3}]));
kg=(10^(Random[Real,{-3,3}]));
sec=(10^(Random[Real,{-3,3}]));
coul=(10^(Random[Real,{-3,3}]));
kelvin=(10^(Random[Real,{-3,3}]));

SECTION 2: UNITS, DERIVED FROM BASE UNITS

cm=(10^(-2)) meter;
mm =(10^(-3))meter;
μm=micron=(10^(-6)) meter; (* character input: esc-m-esc gives μ *)
nm=(10^(-9))meter;
pm=(10^(-12))meter;
fm=(10^(-15))meter;
km=(10^3)meter;
Å=angstrom =(10^(-10))meter; (* character input: esc-Ao-esc gives Å *)
lightyear=9460730472580.8 km;
AstronomicalUnit=149597870691meter;
parsec=(129600/(2π))AstronomicalUnit;
inch=2.54cm;
foot=12inch;
mile=5280foot;
thou=(10^(-3))inch; (* thousandth of an inch; also called mil *)
liter=(10^(-3))(meter^3);
mL=(10^(-3))liter;
μL=(10^(-6))liter;
nL=(10^(-9))liter;
kL=(10^3)liter;
ms=(10^(-3))sec;
μs=(10^(-6))sec;
ns=(10^(-9))sec;
ps=(10^(-12))sec;
fs=(10^(-15))sec;
Hz=(sec^(-1));
mHz=(10^(-3))Hz;
kHz=(10^3)Hz;
MHz=(10^6)Hz;
GHz=(10^9)Hz;
THz=(10^12)Hz;
PHz=(10^15)Hz;
minute=60 sec;
hour=60minute;
SolarDay=24 hour;
SiderealYear=365.2564 SolarDay;
gram=(10^(-3))kg;
amu=1.66053886×(10^(-27))kg;
mg=(10^(-3))gram;
μg=(10^(-6))gram;
tonne=(10^3)kg;
joule=((kg (meter^2))/(sec^2));
mJ=(10^(-3)) joule;
μJ=(10^(-6))joule;
nJ=(10^(-9))joule;
kJ=(10^3)joule;
erg=(10^(-7))joule;
eV=1.60217653×(10^(-19)) joule;
meV=(10^(-3))eV;
keV=(10^3)eV;
MeV=(10^6)eV;
GeV=(10^9)eV;
TeV=(10^12)eV;
btu=1055.056 joule; (*British thermal unit*)
kcal=4184 joule; (*kilocalorie*)
newton=((kg meter)/(sec^2));
dyn=(10^(-5))newton;
pascal=(newton/(meter^2));
kPa=(10^3)pascal;
MPa=(10^6)pascal;
GPa=(10^9)pascal;
bar=(10^5)pascal;
torr=133.322 pascal;
mtorr=(10^(-3))torr;
psi=6894.76 pascal;
atm=101325 pascal;
watt=(joule/sec);
mW=(10^(-3))watt;
μW=(10^(-6))watt;
kW=(10^3)watt;
MW=(10^6)watt;
GW=(10^9)watt;
TW=(10^12)watt;
μC=(10^(-6))coul;
amp=(coul/sec);
mA=(10^(-3))amp;
μA=(10^(-6))amp;
nA=(10^(-9)) amp;
pA=(10^(-12))amp;
fA=(10^(-15))amp;
volt=(joule/coul);
mV=(10^(-3))volt;
μV=(10^(-6))volt;
kV=(10^3)volt;
MV=(10^6)volt;
GV=(10^9)volt;
TV=(10^12)volt;
farad=(coul/volt);
μF=(10^(-6)) farad;
nF=(10^(-9))farad;
pF=(10^(-12))farad;
weber=(joule/amp);
Ω=ohm=(volt/amp); (* character input: esc-O-esc gives Ω (that's the letter O, not the number zero) *)
mΩ=(10^(-3))ohm;
kΩ=(10^3)ohm;
MΩ=(10^6)ohm;
GΩ=(10^9)ohm;
siemens=(ohm^(-1));
mS=(10^(-3))siemens;
μS=(10^(-6))siemens;
nS=(10^(-9))siemens;
tesla=((volt sec)/(meter^2));
mT=(10^(-3))tesla;
μT=(10^(-6))tesla;
nT=(10^(-9))tesla;
gaussSI=(10^(-4))tesla;
henry=(weber/amp);
mH=(10^(-3))henry;
μH=(10^(-6))henry;
nH=(10^(-9))henry;
oerstedSI=(1000/(4π))(amp/meter);
statcoulomb = (((((gram (cm^3))/(sec^2))))^(1/2));
statvolt = (erg/statcoulomb);
gausscgs = (statvolt/cm);
oerstedcgs=gausscgs;

SECTION 3: PHYSICAL CONSTANTS

c0=SpeedOfLight=299792458 (meter/sec); (*speed of light*)
μ0=4π × (10^(-7)) (newton/(amp^2)); (*magnetic constant*)
ϵ0=1/(μ0*(SpeedOfLight^2));(*electric constant*) (* character input: esc-e-esc gives ϵ *)
Z0=μ0*SpeedOfLight; (*impedance of free space*)
hPlanck=6.62606957×(10^(-34))joule sec; (*planck constant*)
ℏ=hbar=(hPlanck/(2π));(*reduced planck constant*) (* character input: esc-hb-esc gives ℏ *)
e=1.60217653×(10^(-19)) coul;(*charge of proton, with underlying SI electromagnetism conventions*)
ecgs = 4.80320427×(10^(-10)) statcoulomb; (*charge of proton, with underlying cgs/Gaussian conventions*)
μBohr=9.27400915×(10^(-24))(joule/tesla);(*Bohr magneton, with underlying SI electromagnetism conventions*)
μBohrcgs=9.27400915×(10^(-21))(erg/oerstedcgs);(*Bohr magneton, with underlying cgs/Gaussian conventions*)
αFS=7.297352568 × (10^(-3)) (*fine-structure constant*) (* character input: esc-a-esc gives α *)
aBohr=0.5291772108×(10^(-10))meter;(*Bohr radius*)
me=9.1093826×(10^(-31))kg;(*electron mass*)
Ry=2.17987197×(10^(-18))joule; (*Rydberg energy, approximately 13.6 eV *)
NA=6.0221415×(10^23);(*Avogadro's number*)
kB=1.3806505×(10^(-23))(joule/kelvin);(*Boltzmann constant*)
GNewton=6.67384×(10^(-11))((meter^3)/(kg (sec^2))); (* "G", the constant of Newtonian gravity*)
σSB=5.670400×(10^(-8))(watt/((meter^2)(kelvin^4)));(*Stefan-Boltzmann constant*) (* character input: esc-s-esc gives σ *)
ARichardson=(4π e me (kB^2))/(hPlanck^3); (*Richardson constant, with underlying SI electromagnetism conventions*)
lPlanck=(Sqrt[ℏ GNewton/(SpeedOfLight^3)]); (*Planck length*)
mPlanck=(Sqrt[ℏ SpeedOfLight/GNewton]); (*Planck mass*)
EPlanck=(Sqrt[ℏ (SpeedOfLight^5)/GNewton]); (*Planck energy*)
tPlanck=(Sqrt[ℏ GNewton/(SpeedOfLight^5)]); (*Planck time*)
TPlanck=(Sqrt[ℏ (SpeedOfLight^5)/GNewton])/kB; (*Planck temperature*)