# File SIunits12.maple for SI units in physics, phys. chemistry, fluid mech.
# J F Harper, School of Mathematics, Statistics and Operations Research,
# Victoria University, Wellington, New Zealand, 31 Jan 2012 
# References:
# NIST10 http://physics.nist.gov/cuu/Constants/index.html as at 31/1/2010 
# HCP10  Handbook of Chemistry and Physics, on-line  , 2010
# JFH10  Harper, J.Colloid Interface Sci. 350, 361-367, 2010
# Ke96b  Kelsall et al., J.Chem.Soc.Faraday Trans. 92, 3887-3893, 1996
# RS70   Robinson & Stokes, Electrolyte Solutions, Butterworth, 1970
# Known problems:
#  If something has numerical value 0 its dimensions can't be found. 
# Recent changes:
# 31/01/12 G,Planck included, NIST updated to NIST10 from NIST07
# 05/05/11 HCP06 data updated to HCP10
# 23/03/11 showvolts, showohms, usevolts, useohms to use volt or ohm
#          instead of A as the electrical unit (show uses A as before)
# 24/11/10 replacing cat by nprintf (frac. powers like 0.333333 not 1/3)
#          Most HCP06 references replaced by HCP10
# 18/11/10 format like m^-1 s instead of s/m; order kg m s A K mol cal
# 08/11/10 Electrothermocapillary calculations
# 15/09/09 qe changed to qE for consistency with JFH10
# 10/08/09 RT, cal (calorie), Rcal, RcalT, Qstar, rhostar, FG added.
num := x->subs(kg=1,m=1,s=1,A=1,volt=1,ohm=1,K=1,mol=1,cal=1,x):
pwr := proc(x,unitx) # returns unitx^its power in x unless it's 0 or 1
# e.g. pwr(m/s,s) is s^-1, pwr(m/s,m) is m, pwr(m/s,kg) is the empty string
  local pwrx:
  if(is(num(x)='0') or is(num(x)='0.'))
  then nprintf("")
  else pwrx := simplify(diff(x,unitx)*unitx/x) assuming  
                   kg>0,m>0,s>0,K>0,mol>0,cal>0,ohm>0,volt>0,A>0;
#       printf("%s%a"," Debug pwrx = """,pwrx,"""");
#       lprint():
       if(is(pwrx='0') or is(pwrx='0.'))
       then nprintf("")	   
       elif(is(pwrx='1') or is(pwrx='1.')) 
       then nprintf(" %s",unitx)
       else nprintf(" %s^%g",unitx,pwrx)
       end if
  end if	
end proc:	   

usevolts := x->subs(A=kg*m^2*s^(-3)/volt,x):
useohms  := x->subs(A=sqrt(kg*m^2*s^(-3)/ohm),x)
            assuming kg>0,m>0,s>0,K>0,mol>0,ohm>0:
# Note: x is, by default, in terms of A    and kg,m,s,K,mol
# but usevolts(x) gives x in terms of volt and kg,m,s,K,mol
# and useohms(x)  gives x in terms of ohm  and kg,m,s,K,mol

units := proc(x,elec) # returns the units of x, not the numerical value
  local localx,ux:    # elec must be v or o for volts or ohms resp.
  if(is(elec='v'))
  then localx := usevolts(x) 
  elif(is(elec='o'))
  then localx := useohms(x)
  else localx := x
  end if:
  if(is(num(localx)='0') or is(num(localx)='0.'))
  then nprintf("units unknown")
  else	   	   
      ux := simplify(localx/num(localx)):
      if (ux=ux^2)
      then nprintf(" dimensionless")
      else 
           if(is(elec='v'))
           then nprintf("%s%s%s%s%s%s%s",pwr(ux,kg),pwr(ux,m),\
             pwr(ux,s),pwr(ux,volt),pwr(ux,K),pwr(ux,mol),pwr(ux,cal))
           elif(is(elec='o'))
	   then nprintf("%s%s%s%s%s%s%s",pwr(ux,kg),pwr(ux,m),\
             pwr(ux,s),pwr(ux,ohm),pwr(ux,K),pwr(ux,mol),pwr(ux,cal))
           else nprintf("%s%s%s%s%s%s%s",pwr(ux,kg),pwr(ux,m),\
             pwr(ux,s),pwr(ux,A),pwr(ux,K),pwr(ux,mol),pwr(ux,cal))
           end if
      end if
  end if
end proc:

show      := x->printf("%12.4e%s \n",num(x),units(x,'A')):
showvolts := x->printf("%12.4e%s \n",num(x),units(x,'v')):
showohms  := x->printf("%12.4e%s \n",num(x),units(x,'o')):
# show gives number in Fortran95 ES12.4 format, units on same line,
# using A for the electrical one. Showvolts, showohms use volt, ohm. 

pi       := evalf(Pi)        :show(pi)      ;# 3.14159...
N        := kg*m/s^2         :show(N)       ;# newton
Pa       := N/m^2            :show(Pa)      ;# pascal
J        := N*m              :show(J)       ;# joule
W        := J/s              :show(W)       ;# watt
C        := A*s              :show(C)       ;# coulomb
V        := W/A              :show(V)       ;# volt    using kg,m,s,A   
Omega    := V/A              :show(Omega)   ;# ohm     using kg,m,s,A   
S        := 1/Omega          :show(S)       ;# siemens using kg,m,s,A
F        := C/V              :show(F)       ;# farad   using kg,m,s,A
H        := Omega*s          :show(H)       ;# henry   using kg,m,s,A
M        := 1e3*mol/m^3      :show(M)       ;# molarity
c        := 299792458*m/s    :show(c)       ;# speed of light   NIST10
faraday  := 96485.3365*C/mol :show(faraday) ;# 2 uncertain figs NIST10
e        := 1.602176565e-19*C:show(e)       ;# 2 uncertain figs NIST10
kBoltz   := 1.3806488e-23*J/K:show(kBoltz)  ;# 2 uncertain figs NIST10
h        := 6.62606957e-34*J*s:show(h)      ;# Planck const NIST10
hbar     := h/(2*pi)         :show(hbar)    ;# reduced Planck const
G        := 6.67384e-11*N*m^2/kg^2:show(G)  ;# Gravitational const NIST10
lP       :=sqrt(hbar*G/c^3)  :show(lP)      ;# Planck length 
mP       :=sqrt(hbar*c/G)    :show(mP)      ;# Planck mass   
TP       :=sqrt(hbar*c^5/G)/kBoltz:show(TP) ;# Planck temperature
tP       :=sqrt(hbar*G/c^5)  :show(tP)      ;# Planck time
avogadro := faraday/e        :show(avogadro);# Avogadro's constant 
gasconst := kBoltz*avogadro  :show(gasconst);# usually called R
mu0      := 4e-7*pi*N/A^2    :show(mu0)     ;# vacuum permeability
epsilon0 := 1.0/(c^2*mu0)    :show(epsilon0);# vacuum permittivity
mobunits := m^2/(e*faraday*Omega*mol):show(mobunits);#mobility RS70 p43
protect('kg','m','s','A','K','mol','volt','ohm','pi','N','Pa','J','W','C',
'V','Omega','S','F','H','M','c','faraday','e','kBoltz','h','hbar','G','lP',
'mP','TP','tP','avogadro','gasconst','mu0','epsilon0','mobunits','cal'):
# All those are basic units or physical constants.
T        := 298.15*K         :show(T)       ;# 25 C as absolute temp.
kT       := kBoltz*T         :show(kT)      ;
RT       := gasconst*T       :show(RT)      ;
Rcal     := gasconst*cal/4.184/J:show(Rcal) ;# Gas constant using calories
RcalT    := Rcal*T           :show(RcalT)   ;# Gas constant*temp. (25 C)
epsH2O   := 78.408*epsilon0  :show(epsH2O);  #H2O permittivity  25 C HCP10   
etaH2O   := 0.89002e-3*Pa*s  :show(etaH2O);  #H2O dynamic visc. 25 C HCP10
rho      := 0.9970470e3*kg/m^3:show(rho);    #H2O density       25 C HCP10
g        := 9.80665*m/s^2    :show(g) ;      #std.grav.accn.    Wikipedia
nu       := etaH2O/rho       :show(nu) ;     #H2O kinemat.visc. 25 C
gamma0   := 71.97e-3*N/m     :show(gamma0);  #H2O sfc tension   25 C HCP10
etaAir   := 20.6e-6*Pa*s     :show(etaAir);  # air dyn.visc. 25 C
mobNa    := 50.08e-4*mobunits:show(mobNa);   # absolute mobility 25C HCP10
mobClO4  := 67.3e-4 *mobunits:show(mobClO4); # absolute mobility 25C HCP10  
DiffNa     := mobNa*kT       :show(DiffNa);  # ion diffusivity
DiffClO4   := mobClO4*kT     :show(DiffClO4);# ion diffusivity
DiffNaClO4 := 2*DiffNa*DiffClO4/(DiffNa+DiffClO4): 
              show(DiffNaClO4);              # NaClO4 diffusivity