# gas.tcl
#
# Library of functions for Gas Dynamic calculations
#
# PJ
#
# Revisions...
# 15-Dec-97 : Initial code
#
#--------------------------------------------------

# Steady isentropic flow: temperature ratio
#
# Input:
# M : Mach number
# g : ratio of specific heats
#
# Output:
# T0_T : ratio of total T over static T
#
proc T0_T { M {g 1.4} } {
   if {$M < 0.0} {
      return -code error "T0_T: invalid Mach number $M"
   }
   expr 1.0 + ($g - 1.0)/2.0 * $M * $M
}; # end of procedure T0_T


# Steady isentropic flow: pressure ratio
#
# Input:
# M : Mach number
# g : ratio of specific heats
#
# Output:
# p0_p : ratio of total pressure over static pressure
#
proc p0_p { M {g 1.4} } {
   if {$M < 0.0} {
      return -code error "p0_p: invalid Mach number $M"
   }
   set T0_T [T0_T $M $g]
   set expon [expr $g / ($g - 1.0)]
   expr pow($T0_T, $expon)
}; # end of procedure p0_p


# Steady isentropic flow: density ratio
#
# Input:
# M : Mach number
# g : ratio of specific heats
#
# Output:
# rho0_rho : ratio of stagnation density over local density
#
proc rho0_rho { M {g 1.4} } {
   if {$M < 0.0} {
      return -code error "rho0_rho: invalid Mach number $M"
   }
   set T0_T [T0_T $M $g]
   set expon [expr 1.0 / ($g - 1.0)]
   expr pow($T0_T, $expon)
}; # end of procedure rho0_rho


# Normal shock: pressure ratio
#
# Input:
# M1 : upstream Mach number
# g  : ratio of specific heats
#
# Output:
# p2_p1 : ratio of downstream pressure over upstream pressure
#
proc p2_p1 { M1 {g 1.4} } {
   if {$M1 < 1.0} {
      return -code error "p2_p1: invalid Mach number $M1"
   }
   expr 1.0 + 2.0 * $g * ($M1 * $M1 - 1.0) / ($g + 1.0)
}; # end of procedure p2_p1


# Normal shock: density ratio
#
# Input:
# M1 : upstream Mach number
# g  : ratio of specific heats
#
# Output:
# rho2_rho1 : ratio of downstream density over upstream density
#
proc rho2_rho1 { M1 {g 1.4} } {
   if {$M1 < 1.0} {
      return -code error "rho2_rho1: invalid Mach number $M1"
   }
   set numer [expr ($g + 1.0) * $M1 * $M1]
   set denom [expr 2.0 + ($g - 1.0) * $M1 * $M1]
   expr $numer / $denom
}; # end of procedure rho2_rho1


# Normal shock: velocity ratio
#
# Input:
# M1 : upstream Mach number
# g  : ratio of specific heats
#
# Output:
# u2_u1 : ratio of downstream velocity over upstream velocity
#         The velocities are in the shock-stationary frame.
#
proc u2_u1 { M1 {g 1.4} } {
   if {$M1 < 1.0} {
      return -code error "u2_u1: invalid Mach number $M1"
   }
   set rho2_rho1 [rho2_rho1 $M1 $g]
   expr 1.0 / $rho2_rho1
}; # end of procedure u2_u1


# Normal shock: temperature ratio
#
# Input:
# M1 : upstream Mach number
# g  : ratio of specific heats
#
# Output:
# T2_T1 : ratio of downstream temperature over upstream temperature
#
proc T2_T1 { M1 {g 1.4} } {
   if {$M1 < 1.0} {
      return -code error "T2_T1: invalid Mach number $M1"
   }
   set rho2_rho1 [rho2_rho1 $M1 $g]
   set p2_p1 [p2_p1 $M1 $g]
   expr $p2_p1 / $rho2_rho1
}; # end of procedure T2_T1


# Normal shock: downstream Mach number
#
# Input:
# M1 : upstream Mach number
# g  : ratio of specific heats
#
# Output:
# Mach2 : downstream Mach number in shock-stationary frame
#
proc Mach2 { M1 {g 1.4} } {
   if {$M1 < 1.0} {
      return -code error "Mach2: invalid Mach number $M1"
   }
   set numer [expr 1.0 + 0.5 * ($g - 1.0) * $M1 * $M1]
   set denom [expr $g * $M1 * $M1 - 0.5 * ($g - 1.0)]
   expr sqrt($numer / $denom)
}; # end of procedure Mach2


# Rayleigh Pitot pressure formula
#
# Input:
# M1 : upstream Mach number
# g  : ratio of specific heats
#
# Output:
# p02_p1 : ratio of downstream total pressure over
#          upstream static pressure
#
proc p02_p1 { M1 {g 1.4} } {
   if {$M1 < 0.0} {
      return -code error "p02_p1: invalid Mach number $M1"
   }
   if {$M1 > 1.0} {
      # Go across the normal shock then bring to rest isentropically.
      set p2_p1  [p2_p1 $M1 $g]
      set m2     [Mach2 $M1 $g]
      set p02_p2 [p0_p $m2 $g]
      set p02_p1 [expr $p02_p2 * $p2_p1]
   } else {
      # The whole process is isentropic.
      set p02_p1 [p0_p $M1 $g]
   }
}; # end of procedure p02_p1

# gascalc.tcl
# Calculator for gas-dynamic relations.
#
# P.Jacobs
# Department of Mechanical Enginering
# The University of Queensland
#
# Version...
# 28-Jan-98 : initial code to try out wish and browser plugin
#
#-------------------------------------------------------------

# source gas.tcl

proc initGasCalc {} {
    global mach
    global g
    global isentropic
    global shock

    set mach 1.5
    set g    1.4
    updateRelations

    frame .fTitle
    label .fTitle.lTitle -text "Gas-Dynamic Relations"
    pack .fTitle.lTitle
    pack .fTitle -side top

    # Set up two entry fields for the Mach number and
    # ratio-of-specific-heats.

    frame .fInput
    label .fInput.lTitle -text Inputs:
    label .fInput.lGamma -text Gamma -anchor w
    entry .fInput.eGamma -textvariable g -relief sunken -width 12
    label .fInput.lMach  -text Mach -anchor w
    entry .fInput.eMach  -textvariable mach -relief sunken -width 12
    pack  .fInput.lTitle -side left
    pack  .fInput.lGamma .fInput.eGamma -side left
    pack  .fInput.lMach .fInput.eMach -side left
    pack  .fInput -side top -pady 2

    # Set bindings so that the gas-dynamic relations are
    # updated as soon as a new value is entered.

    bind  .fInput.eGamma  updateRelations
    bind  .fInput.eMach   updateRelations

    # A table of isentropic results

    frame .fIsentropic -relief sunken -borderwidth 2
    label .fIsentropic.lTitle -text "Isentropic Flow"
    frame .fIsentropic.fResults
    foreach name [array names isentropic] {
       label .fIsentropic.fResults.l$name -text $name
       label .fIsentropic.fResults.r$name -textvariable isentropic($name)
       grid  .fIsentropic.fResults.l$name .fIsentropic.fResults.r$name -sticky w
    }
    pack  .fIsentropic.lTitle -side left
    pack  .fIsentropic.fResults
    pack  .fIsentropic -fill x -pady 2

    # A table of shock results.

    frame .fShock -relief sunken -borderwidth 2
    label .fShock.lTitle -text "Normal Shock"
    frame .fShock.fResults
    foreach name [array names shock] {
       label .fShock.fResults.l$name -text $name
       label .fShock.fResults.r$name -textvariable shock($name)
       grid  .fShock.fResults.l$name .fShock.fResults.r$name -sticky w
    }
    pack  .fShock.lTitle -side left
    pack  .fShock.fResults
    pack  .fShock -fill x -pady 2

    # And, at the bottom of the screen, a Dismiss button
    # to shut everything down.

    button .bDismiss -text "Dismiss" -command {closeGasCalc}
    pack   .bDismiss
}; # end of initGasCalc


proc closeGasCalc {} {
    global isentropic
    global shock

    destroy .fTitle.lTitle
    destroy .fTitle

    destroy .fInput.lTitle
    destroy .fInput.lGamma
    destroy .fInput.eGamma
    destroy .fInput.lMach
    destroy .fInput.eMach
    destroy .fInput

    foreach name [array names isentropic] {
       destroy .fIsentropic.fResults.l$name
       destroy .fIsentropic.fResults.r$name
    }
    destroy .fIsentropic.fResults
    destroy .fIsentropic.lTitle
    destroy .fIsentropic

    foreach name [array names shock] {
       destroy .fShock.fResults.l$name
       destroy .fShock.fResults.r$name
    }
    destroy .fShock.fResults
    destroy .fShock.lTitle
    destroy .fShock

    destroy .bDismiss

    unset isentropic
    unset shock
}; # end of closeGasCalc

proc updateRelations {} {
    global mach
    global g
    global isentropic
    global shock

    set isentropic(p0/p) [p0_p $mach $g]
    set isentropic(T0/T) [T0_T $mach $g]
    set isentropic(rho0/rho) [rho0_rho $mach $g]

    if {$mach > 1.0} {
       set shock(p2/p1) [p2_p1 $mach $g]
       set shock(T2/T1) [T2_T1 $mach $g]
       set shock(rho2/rho1) [rho2_rho1 $mach $g]
       set shock(Mach2) [Mach2 $mach $g]
       set shock(p02/p1) [p02_p1 $mach $g]
    } else {
       set shock(p2/p1) "---"
       set shock(T2/T1) "---"
       set shock(rho2/rho1) "---"
       set shock(Mach2) "---"
       set shock(p02/p1) "---"
    }; # end if
}; # end of updateRelations


# -----------------------------------------------------
# Now that we have set everything up, start building...
# -----------------------------------------------------

initGasCalc