Updated 2013-01-18 04:29:05 by pooryorick

Arjen Markus (4 may 2006) It occurred to me that we do not yet have a decent module in Tcllib yet that collects the various physical empirical formulae. So, right now, if I want to compute the water content of air, I have to hunt the Internet for the right formula. Similarly, the density of water as a function of temperature and salinity.

I am probably not the only one with that little problem, so I propose to collect such formulae here (with a reference) so that later we can simply wrap them up in a package/module for Tcllib.

Note:

As kbk pointed out, some or all of these formulae will have to be rearranged so that they can be evaluated in a fast and numerically stable way. For now, we are just gathering them in one place.
```# Auxiliary procedures:
#
# Check the range of a parameter
#
proc check {descr var range} {
foreach {min max} \$range break
if { \$var < \$min || \$var > \$max } {
return -code error "Range error: \$descr should be between \$min and \$max"
}
}

# Efficiently evaluate a polynomial -
# coefs is a list with the highest power first
proc poleval {coefs x} {
set y 0
foreach c \$coefs {
set y [expr {\$y * \$x + \$c}]
}
}```

A formula for the humidity content of air:
```check "dew-point temperature" \$td {-20 100}
set vapordens [expr {5.018+0.32321*\$td+8.1847e-3*\$td*\$td+3.1243e-4*\$td*\$td*\$td}]```

Symbols:

• td the dewpoint temperature in degrees C - between -20 and +100 degrees
• vapordens the water vapour density in g water/m3

Reference: http://hyperphysics.phy-astr.gsu.edu/hbase/kinetic/relhum.html#c3

TR - Here comes a formula for the density of seawater as a function of salinity, temperature, and pressure:
```proc rho {Salinity Temperature Pressure} {
#
# compute the density of seawater in kg/m3 as a function of:
#
# Salinity    -> salinity (in psu, practical salinity units)
# Temperature -> temperature (in ℃, degrees Celsius)
# Pressure    -> pressure (in bar)
#
check Temperature \$Temperature {-2 40}
check Salinity \$Salinity {0 40}
check Pressure \$Pressure {0 12000}
set rhow [expr {
999.842594
+ 0.06793952 * \$Temperature
- 0.00909529 * pow(\$Temperature,2)
+ 0.0001001685 * pow(\$Temperature,3)
- 1.120083e-06 * pow(\$Temperature,4)
+ 6.536332e-09 * pow(\$Temperature,5)
}]
set A [expr {
0.824493
- 0.0040899 * \$Temperature
+ 7.6438e-05 * pow(\$Temperature,2)
- 8.2467e-07 * pow(\$Temperature,3)
+ 5.3875e-09 * pow(\$Temperature,4)
}]
set B [expr {
-0.00572466
+ 0.00010227 * \$Temperature
- 1.6546e-06 * pow(\$Temperature,2)
}]
set C 0.00048314
set rho0 [expr {
\$rhow
+ \$A * \$Salinity
+ \$B * pow(\$Salinity, 3.0/2)
+ \$C * pow(\$Salinity, 2)
}]
set Ksbmw [expr {
19652.21
+ 148.4206 * \$Temperature
- 2.327105 * pow(\$Temperature,2)
+ 0.01360477 * pow(\$Temperature,3)
- 5.155288e-05 * pow(\$Temperature,4)
}]
set Ksbm0 [expr {
\$Ksbmw
+ \$Salinity * (
54.6746
- 0.603459 * \$Temperature
+ 0.0109987 * pow(\$Temperature,2)
- 6.167e-05 * pow(\$Temperature,3))
+ pow(\$Salinity,3.0/2) * (
0.07944
+ 0.016483 * \$Temperature
- 0.00053009 * pow(\$Temperature,2))
}]
set Ksbm [expr {
\$Ksbm0
+ \$Pressure * (
3.239908
+ 0.00143713 * \$Temperature
+ 0.000116092 * pow(\$Temperature,2)
- 5.77905e-07 * pow(\$Temperature,3))
+ \$Pressure * \$Salinity * (
0.0022838
- 1.0981e-05 * \$Temperature
- 1.6078e-06 * pow(\$Temperature,2))
+ \$Pressure * pow(\$Salinity,3.0/2) * 0.000191075
+ \$Pressure * \$Pressure * (
8.50935e-05
- 6.12293e-06 * \$Temperature
+ 5.2787e-08 * pow(\$Temperature,2))
+ pow(\$Pressure,2) * \$Salinity * (
-9.9348e-07
+ 2.0816e-08 * \$Temperature
+ 9.1697e-10 * pow(\$Temperature,2))
}]
return [expr {\$rho0/double(1 - double(\$Pressure)/\$Ksbm)}]
}```

or using poleval
```proc rho {Salinity Temperature Pressure} {
set rhow [poleval {+6.536332e-09 -1.120083e-06 +0.0001001685 -0.00909529 +0.06793952 999.842594} \$Temperature]
set A [poleval {+5.3875e-09 -8.2467e-07 +7.6438e-05 -0.0040899 0.824493} \$Temperature]
set B [poleval {-1.6546e-06 +0.00010227 -0.00572466} \$Temperature]
set C 0.00048314
set rho0 [expr {\$rhow + \$A * \$Salinity + \$B * pow(\$Salinity,3.0/2) + \$C * pow(\$Salinity,2)}]
set Ksbmw [poleval {-5.155288e-05 +0.01360477 -2.327105 +148.4206 +19652.21} \$Temperature]
set Ksbm0 [expr {
\$Ksbmw
+ \$Salinity * (
54.6746
- 0.603459 * \$Temperature
+ 0.0109987 * pow(\$Temperature,2)
- 6.167e-05 * pow(\$Temperature,3))
+ pow(\$Salinity,3.0/2) * (
0.07944
+ 0.016483 * \$Temperature
- 0.00053009 * pow(\$Temperature,2))
}]
set Ksbm [expr {
\$Ksbm0
+ \$Pressure * (
3.239908
+ 0.00143713 * \$Temperature
+ 0.000116092 * pow(\$Temperature,2)
- 5.77905e-07 * pow(\$Temperature,3))
+ \$Pressure * \$Salinity * (
0.0022838
- 1.0981e-05 * \$Temperature
- 1.6078e-06 * pow(\$Temperature,2))
+ \$Pressure * pow(\$Salinity,3.0/2) * 0.000191075
+ \$Pressure * \$Pressure * (
8.50935e-05
- 6.12293e-06 * \$Temperature
+ 5.2787e-08 * pow(\$Temperature,2))
+ pow(\$Pressure,2) * \$Salinity * (
-9.9348e-07
+ 2.0816e-08 * \$Temperature
+ 9.1697e-10 * pow(\$Temperature,2))
}]
return [expr {\$rho0/double(1 - double(\$Pressure)/\$Ksbm)}]
}```

References:

Note: I am not totally sure about the ranges of the three parameters, so I added some reasonable defaults

Example: [rho 35 10 0] = 1026.95 (quite standard seawater with a salinity of 35 psu, temperature of 10 ℃ and pressure of 0 bar (surface water))

IDG May 04 2006: All empirical formulae need to have their range of validity stated. Note what happens to the above as td -> infinity :-)

AM Ah! Good thinking! I have added the valid range.

TR - Perhaps we should also add something like accuracy. It is often not good to give a result like 3.24345643 when the formula only has significance for 3.24. And: we should stick to SI units, so one result can be fed into the next formula without conversion.

Sarnold - Add something like accuracy? Take a look at math::bigfloat if you need to manage accuracy across computations. I must admit it may obfuscate the source code for such a package. (When I designed it, I remembered my physicist background...)

As a small contribution to this project, I was thinking that it could be useful to have a "formula" proc that created a proc that applied the formula, converting units if necessary. Here's my whack at a "formula" proc:
```package require units

# Author: Eric Kemp-Benedict, 2007
# Use it as you like!
# Usage:
#   formula nameOfFormula "Description" formula units V1 U1 V2 U2 ...
proc formula {name args} {
set fargs \$args
set procbody1 "
set name \$name
set descr {[lindex \$fargs 0]}
set formula {[lindex \$fargs 1]}
set funits {[lindex \$fargs 2]}
foreach {v u} [list [lrange \$fargs 3 end]] {
set units(\\$v) \\$u
}"
set procbody2 {
if {\$args eq "info"} {
set infostr "\$descr\n\nVariables:\n"
foreach v [lsort [array names units]] {
set infostr "\$infostr   \$v (\$units(\$v))\n"
}
return \$infostr
}
foreach {v = val} \$args {
if {[lsearch [array names units] \$v] == -1} {
error "\$name: Variable \$v is not valid; Use \"\$name info\" for more information"
return 0.0
}
set \$v [::units::convert \$val \$units(\$v)]
}
return "[expr \$formula] \$funits"
}

set procbody "\$procbody1\n\$procbody2"

proc \$name args \$procbody
}

##
## Example
##
formula CharlesVol "Apply Charles' Law to calculate volume" \
{1.0 * \$V1 * \$T2/\$T1} L \
V1 L T1 K T2 K```

Once the example is created, you can do:
```Apply Charles' Law to calculate volume

Variables:
T1 (K)
T2 (K)
V1 (L)```

to get a summary of the formula, and something like
```(Formula) 98 % CharlesVol V1 = 3.2m^3 T1 = 300K T2 = 298K
3178.66666667 L```

to apply it.

AM (5 july 2007) Now, that looks very nice, especially the "info" subcommand! EKB Thanks!