[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:
   * Millero, F. J. and Poisson, A. 1981 International one-atmosphere equation of state of seawater. Deep-Sea Research 28A, 625-629.
   * http://cran.r-project.org/src/contrib/Descriptions/seacarb.html

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:
 (Formula) 97 % CharlesVol info
 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! 
----
[[[Category Physics]|[Category Science]]]