# physics.tcl --
 #    Provide access to physical constants, relations and
 #    formulas
 #
 #    Simple example:
 #    Conversion of temperature to absolute temperature
 #    In many formulas you need to specify the absolute temperature
 #    rather than the temperature most people are used to (degrees
 #    centigrade or the infamous degrees Fahrenheit).
 #    Mistakes are easily made ... so include the unit in the value.
 #
 #    Note:
 #    In this implementation I have been less than consistent
 #    in the use of unit conversions. I mainly wanted to get
 #    something working ;)
 #

 # physics --
 #    Namespace for the physical computations
 #
 namespace eval ::physics {

    #
    # Add common constants etc.
    #
    variable gas_constant 8.3143 ;# J/mol.K

    namespace export abs_temperature     wien_wavelength \
                     velocity_sound_gas  diffusivePoint1d \
                     superpose
 }

 # Conversions --
 #    Set of conversion procedures
 #
 proc ::physics::Conv_K_K   {t} {return $t}
 proc ::physics::Conv_oC_K  {t} {expr {$t+273.15}}
 proc ::physics::Conv_oF_K  {t} {Conv_oC_K [Conv_oF_oC $t]}
 proc ::physics::Conv_oF_oC {t} {expr {5.0*$t/9.0+32.0}}

 # abs_temperature --
 #    Return the absolute temperature
 #
 # Arguments:
 #    temp         (Relative) temperature
 #
 # Result:
 #    The temperature in kelvin (no capital, there was only one
 #    lord Kelvin, AFAIK :)
 #
 proc ::physics::abs_temperature {temp} {

    if { [llength $temp] == 2 } {
       set temp [Conv_[lindex $temp 1]_K [lindex $temp 0]]
    }
    list $temp K
 }

 # wien_wavelength --
 #    Return the wave length at which a perfect black body emits most
 #    energy, given the temperature
 #
 # Arguments:
 #    temp         (Relative) temperature
 #
 # Result:
 #    The wave length in micrometers (um)
 #
 proc ::physics::wien_wavelength {temp} {
    list [expr {2898.0/[lindex [abs_temperature $temp] 0]}] um
 }

 # velocity_sound_gas --
 #    Return the (approximate) velocity of sound in gas,
 #    given the kind of gas and the temperature
 #
 # Arguments:
 #    cpcv         Ratio of cp and cv or "mono", "di", "poly" for
 #                 the type of molecules the gas consists of
 #    molarw       Molar weight (g/mol)
 #    temp         Temperature (K)
 # Result:
 #    The sound velocity in m/s
 #
 proc ::physics::velocity_sound_gas {cpcv molarw temp} {
    variable gas_constant

    set abstemp [lindex [abs_temperature $temp] 0]
    set molarw  [expr {$molarw/1000.0}]
    switch -- $cpcv {
    "mono" { set cpcv 1.66 }
    "di"   { set cpcv 1.4 }
    "poly" { set cpcv 1.3 }
    }
    set csound [expr {sqrt($cpcv*$gas_constant*$abstemp/$molarw)}]
    list $csound "m/s"
 }

 # diffusivePoint1d --
 #    Create a function that describes the concentration from a
 #    point source under diffusion in one dimension
 #
 # Arguments:
 #    name         Name of the function to create
 #    diff         Diffusion coefficient (m/s2)
 #    mass         Mass that was released
 #    pos          Position of the release point
 # Result:
 #    Name of the function
 #
 proc ::physics::diffusivePoint1d {name diff mass pos} {
    proc ::physics::$name {x t} \
       [string map [list %diff $diff %mass $mass %pos $pos] {
       expr {%mass/sqrt(4.0*3.1415926*%diff*$t)*exp(-($x-%pos)*($x-%pos)/4.0/$t/%diff)}
    }]
    return ::physics::$name
 }

 # superpose --
 #    Create a function that computes the sum of the individual
 #    functions - the superposition principle
 #
 # Arguments:
 #    name         Name of the function to create
 #    args         List of one or more functions to superpose
 # Result:
 #    Name of the function
 #
 proc ::physics::superpose {name args} {
    set body "expr {"
    set arglist [info args [lindex $args 0]]
    set values  "\$[join $arglist " $"]"

    foreach fnc $args {
       append body "+\[$fnc $values\]"
    }
    append body "}"

    proc ::physics::$name $arglist $body
    return ::physics::$name
 }

 # main --
 #    A simple demonstration
 #
 namespace import ::physics::*

 puts "20 degrees centigrade = [abs_temperature {20 oC}]"
 puts "50 degrees Fahrenheit = [abs_temperature {50 oF}]"
 puts "At 20 degrees centigrade a black body emits most
 energy with a wave length of [wien_wavelength {20 oC}]"

 puts "Velocity of sound in air at 20 oC: \
 [velocity_sound_gas di 28.9 {20 oC}]"

 #
 # Use a fairly large diffusion coefficient:
 # otherwise the concentration will be nearly zero
 # and we risk underflow
 # (I ought to use larger times ;)
 #
 set pointsrc  [diffusivePoint1d point1 1.0 1.0e3 0.0]
 set pointsrc2 [diffusivePoint1d point2 1.0 1.0e3 4.0]
 foreach t {0.1 0.3 1.0 3.0 10.0} {
    puts "t=$t: concentration = [$pointsrc 1.0 $t]"
 }
 foreach x {0.1 0.3 1.0 2.0 3.0} {
    puts "x=$x: concentration = [$pointsrc $x 1.0]"
 }

 puts "Two sources:"
 set all [superpose twopoints $pointsrc $pointsrc2]
 foreach x {0.1 0.3 1.0 2.0 3.0 3.9 4.0 4.1 5.0 6.0} {
    puts "x=$x: concentration = [$all $x 1.0] (= sum of [$pointsrc $x 1.0] and [$pointsrc2 $x 1.0])"
 }