Package aerocalc :: Module std_atm
[frames] | no frames]

Source Code for Module aerocalc.std_atm

   1  #!/usr/bin/python 
   2  # -*- coding: utf-8 -*- 
   3   
   4  # ############################################################################# 
   5  # Copyright (c) 2008, Kevin Horton 
   6  # All rights reserved. 
   7  # Redistribution and use in source and binary forms, with or without 
   8  # modification, are permitted provided that the following conditions are met: 
   9  # * 
  10  #     * Redistributions of source code must retain the above copyright 
  11  #       notice, this list of conditions and the following disclaimer. 
  12  #     * Redistributions in binary form must reproduce the above copyright 
  13  #       notice, this list of conditions and the following disclaimer in the 
  14  #       documentation and/or other materials provided with the distribution. 
  15  #     * The name of Kevin Horton may not be used to endorse or promote products 
  16  #       derived from this software without specific prior written permission. 
  17  # * 
  18  # THIS SOFTWARE IS PROVIDED BY KEVIN HORTON ``AS IS'' AND ANY EXPRESS OR 
  19  # IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF 
  20  # MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO 
  21  # EVENT SHALL KEVIN HORTON BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 
  22  # SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 
  23  # PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; 
  24  # OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
  25  # WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR 
  26  # OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF 
  27  # ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
  28  # ############################################################################# 
  29  # 
  30  # version 0.16, 06 May 2007 
  31  # 
  32  # Version History: 
  33  # vers     date     Notes 
  34  #  0.1   14 May 06  First release. 
  35  # 
  36  # 0.11   17 May 06  Cleaned up documentation. 
  37  # 
  38  # 0.12   24 May 06  Added temp2speed_of_sound function. 
  39  # 
  40  # 0.13   02 Jun 06  Added temp2isa and isa2temp functions. 
  41  # 
  42  # 0.14   22 Apr 07  Added density altitude vs temperature functions 
  43  # 
  44  # 0.15   29 Apr 07  Broke out sat_press and dry_press as public functions. 
  45  # 
  46  # 0.16   05 May 07  Reworked to use default units from default_units module. 
  47  # ############################################################################# 
  48  # 
  49  # To Do: 1. Done. 
  50  # 
  51  #        2. Done. 
  52  # 
  53  #        3. Done. 
  54  # 
  55  #        4. Won't do. 
  56  # 
  57  #        5. Add temp2temp_ratio, press2press_ratio and density2density_ratio. 
  58  # 
  59  #        6. Done. 
  60  # 
  61  #        7. Add tests for all functions to test/test_std_atm.py 
  62  #           Tests to add: 
  63  #           dry_press 
  64  #           density_alt_table ? (probably won't make a test for this) 
  65  # 
  66  #        8. Review code to catch duplicated code blocks in different functions. 
  67  #           Move these blocks to internal functions. 
  68  # 
  69  #        9. Review API for all public functions for consistency of units, etc. 
  70  # 
  71  # Done:  1. consider replacing calculations by constants where possible.  See 
  72  #           http://jcoppens.com/globo/teoria/atm_est.en.php 
  73  # 
  74  #           Tested replacing calculations by constants in press2alt.  The 
  75  #           perf improvement was only about 15%.  Probably not worth the trouble. 
  76  #           Better to keep the pedigree of the equations visible. 
  77  # 
  78  #        2. Validate against published data.  Created unittests using data: 
  79  #           http://www.sworld.com.au/steven/space/atmosphere/ 
  80  # 
  81  #        3. Added relative humidity to density altitude calc.  See: 
  82  #           http://wahiduddin.net/calc/density_altitude.htm 
  83  # 
  84  #        4. Change formulae to use pressure in pa, not in HG.  Won't do. 
  85  #           Instead, changed to use default units specified in default_units.py 
  86  # 
  87  #        6. Added functions: 
  88  #           isa2temp 
  89  #           temp2isa 
  90  # 
  91  #        7. Added tests for functions: 
  92  #           pressure_alt 
  93  #           sat_press 
  94  #           density_alt2temp 
  95  # ############################################################################# 
  96   
  97  """Calculate standard atmosphere parametres. 
  98   
  99  Calculates standard atmosphere parametres, using the 1976 International 
 100  Standard Atmosphere.  The default units for the input and output are defined 
 101  in default_units.py 
 102       
 103  All altitudes are geopotential altitudes (i.e. it is assumed that there is 
 104  no variation with altitude of the acceleration due to gravity). 
 105       
 106  Works up to 84.852 km (278,386 ft) altitude. 
 107   
 108  """ 
 109   
 110  import unit_conversion as U 
 111  import math as M 
 112  import locale as L 
 113  try: 
 114      from default_units import * 
 115  except ImportError: 
 116      default_area_units = 'ft**2' 
 117      default_power_units = 'hp' 
 118      default_speed_units = 'kt' 
 119      default_temp_units = 'C' 
 120      default_weight_units = 'lb' 
 121      default_press_units = 'in HG' 
 122      default_density_units = 'lb/ft**3' 
 123      default_length_units = 'ft' 
 124      default_alt_units = default_length_units 
 125      default_avgas_units = 'lb' 
 126   
 127  try: 
 128      L.setlocale(L.LC_ALL, 'en_US') 
 129  except: 
 130      pass 
 131   
 132  g = 9.80665  # Acceleration of gravity at 45.542 deg latitude, m/s**s 
 133  Rd = 287.05307  # Gas constant for dry air, J/kg K 
 134   
 135  # conditions starting at sea level, in a region with temperature gradient 
 136   
 137  T0 = 288.15  # Temperature at sea level, degrees K 
 138  L0 = -6.5  # Temperature lapse rate, at sea level deg K/km 
 139  P0 = 29.9213  # Pressure at sea level, in HG 
 140  Rho0 = 1.2250  # Density at sea level, kg/m**3 
 141   
 142  # conditions starting at 11 km, in an isothermal region 
 143   
 144  T11 = T0 + 11 * L0  # Temperature at 11,000 m, degrees K 
 145  PR11 = (T11 / T0) ** ((-1000 * g) / (Rd * L0))  # pressure ratio at 11,000 m 
 146  P11 = PR11 * P0 
 147  Rho11 = (Rho0 * PR11) * (T0 / T11) 
 148   
 149  # conditions starting at 20 km, in a region with temperature gradient 
 150   
 151  T20 = T11 
 152  PR20 = PR11 * M.exp(((-1000 * g) * (20 - 11)) / (Rd * T11)) 
 153  L20 = 1  # temperature lapse rate, starting at 20,000 m, deg K/km 
 154  P20 = PR20 * P0 
 155  Rho20 = (Rho0 * PR20) * (T0 / T20) 
 156   
 157  # conditions starting at 32 km, in a region with temperature gradient 
 158   
 159  T32 = 228.65  # Temperature at 32 km, degrees K 
 160  PR32 = PR20 * (T32 / T20) ** ((-1000 * g) / (Rd * L20)) 
 161   
 162  # PR32 = PR20 * M.exp((-1000 * g) * (32 - 20)/(R * T20)) 
 163   
 164  L32 = 2.8  # temperature lapse rate, starting at 32,000 m, deg K/km 
 165  P32 = PR32 * P0 
 166  Rho32 = (Rho0 * PR32) * (T0 / T32) 
 167   
 168  # conditions starting at 47 km, in an isothermal region 
 169   
 170  T47 = 270.65 
 171  PR47 = PR32 * (T47 / T32) ** ((-1000 * g) / (Rd * L32)) 
 172  P47 = PR47 * P0 
 173  Rho47 = (Rho0 * PR47) * (T0 / T47) 
 174   
 175  # conditions starting at 51 km, in a region with temperature gradient 
 176   
 177  T51 = 270.65  # Temperature at 51 km, degrees K 
 178  PR51 = PR47 * M.exp(((-1000 * g) * (51 - 47)) / (Rd * T47)) 
 179  L51 = -2.8  # temperature lapse rate, starting at 51,000 m, deg K/km 
 180  P51 = PR51 * P0 
 181  Rho51 = (Rho0 * PR51) * (T0 / T51) 
 182   
 183  # conditions starting at 71 km, in a region with temperature gradient 
 184   
 185  T71 = 214.65  # Temperature at 71 km, degrees K 
 186  PR71 = PR51 * (T71 / T51) ** ((-1000 * g) / (Rd * L51)) 
 187  L71 = -2.  # temperature lapse rate, starting at 71,000 m, deg K/km 
 188  P71 = PR71 * P0 
 189  Rho71 = (Rho0 * PR71) * (T0 / T71) 
 190   
 191  # temp_units_list = ['C', 'F', 'K', 'R'] 
 192   
 193  # ############################################################################# 
 194  # 
 195  # Altitude to temperature 
 196  # 
 197  # ############################################################################# 
 198   
 199   
200 -def alt2temp(H, alt_units=default_alt_units, 201 temp_units=default_temp_units):
202 """Return the standard temperature for the specified altitude. Altitude 203 units may be feet ('ft'), metres ('m'), statute miles, ('sm') or 204 nautical miles ('nm'). Temperature units may be degrees C, F, K or R 205 ('C', 'F', 'K' or 'R') 206 207 If the units are not specified, the units in default_units.py are used. 208 209 Examples: 210 211 Calculate the standard temperature (in default temperature units) at 212 5,000 (default altitude units): 213 >>> alt2temp(5000) 214 5.0939999999999941 215 216 Calculate the standard temperature in deg F at sea level: 217 >>> alt2temp(0, temp_units = 'F') 218 59.0 219 220 Calculate the standard temperature in deg K at 11,000 m: 221 >>> alt2temp(11000, alt_units = 'm', temp_units = 'K') 222 216.64999999999998 223 224 Calculate the standard temperature at 11 statute miles in deg R: 225 >>> alt2temp(11, alt_units = 'sm', temp_units = 'R') 226 389.96999999999997 227 228 The input value may be an expression: 229 >>> alt2temp(11 * 5280, temp_units = 'R') 230 389.96999999999997 231 232 """ 233 234 # Validated to 84000 m 235 # uses meters and degrees K for the internal calculations 236 237 # function tested in tests/test_std_atm.py 238 239 H = U.len_conv(H, from_units=alt_units, to_units='km') 240 241 if H <= 11: 242 temp = T0 + H * L0 243 elif H <= 20: 244 temp = T11 245 elif H <= 32: 246 temp = T20 + (H - 20) * L20 247 elif H <= 47: 248 temp = T32 + (H - 32) * L32 249 elif H <= 51: 250 temp = T47 251 elif H <= 71: 252 temp = T51 + (H - 51) * L51 253 elif H <= 84.852: 254 temp = T71 + (H - 71) * L71 255 else: 256 raise ValueError, \ 257 'This function is only implemented for altitudes of 84.852 km and below.' 258 259 return U.temp_conv(temp, to_units=temp_units, from_units='K')
260 261
262 -def alt2temp_ratio(H, alt_units=default_alt_units):
263 """ 264 Return the temperature ratio (temperature / standard temperature for 265 sea level). The altitude is specified in feet ('ft'), metres ('m'), 266 statute miles, ('sm') or nautical miles ('nm'). 267 268 If the units are not specified, the units in default_units.py are used. 269 270 Examples: 271 272 Calculate the temperature ratio at 8,000 (default altitude units) 273 >>> alt2temp_ratio(8000) 274 0.94499531494013533 275 276 Calculate the temperature ratio at 8,000 m. 277 >>> alt2temp_ratio(8000, alt_units = 'm') 278 0.81953843484296374 279 """ 280 281 # function tested in tests/test_std_atm.py 282 283 return alt2temp(H, alt_units, temp_units='K') / T0
284 285 286 # ############################################################################# 287 # 288 # ISA deviation to temperature 289 # 290 # ############################################################################# 291 292
293 -def isa2temp( 294 ISA_dev, 295 altitude, 296 temp_units=default_temp_units, 297 alt_units=default_alt_units, 298 ):
299 """ 300 Return the temperature that is a specified amount warmer or cooler 301 than the standard temperature for the altitude. 302 303 The temperature may be in deg C, F, K or R. 304 305 The altitude may be in feet ('ft'), metres ('m'), kilometres ('km'), 306 statute miles, ('sm') or nautical miles ('nm'). 307 308 If the units are not specified, the units in default_units.py are used. 309 310 Examples: 311 312 Determine the temperature that is 10 deg (default temperature units) warmer 313 than the standard temperature at 8,000 (default altitude units): 314 >>> isa2temp(10, 8000) 315 9.1503999999999905 316 317 Determine the temperature that is 25 degrees K cooler than the standard 318 temperature at 2000 m. 319 >>> isa2temp(-25, 2000, temp_units = 'K', alt_units = 'm') 320 250.14999999999998 321 """ 322 323 # function tested in tests/test_std_atm.py 324 325 temp = ISA_dev + alt2temp(altitude, alt_units, temp_units) 326 327 return temp
328 329 330 # ############################################################################# 331 # 332 # temperature to ISA deviation 333 # 334 # ############################################################################# 335 336
337 -def temp2isa( 338 temp, 339 altitude, 340 temp_units=default_temp_units, 341 alt_units=default_alt_units, 342 ):
343 """ 344 Return the amount that the specified temperature is warmer or cooler 345 than the standard temperature for the altitude. 346 347 The temperature may be in deg C, F, K or R. 348 349 The altitude may be in feet ('ft'), metres ('m'), kilometres ('km'), 350 statute miles, ('sm') or nautical miles ('nm'). 351 352 If the units are not specified, the units in default_units.py are used. 353 354 Examples: 355 356 Determine the ISA deviation for a temperature of 30 deg (default 357 temperature units) at an altitude of 2000 (default altitude units): 358 >>> temp2isa(30, 2000) 359 18.962400000000002 360 361 Determine the ISA deviation in degrees F for a temperature of 45 deg F 362 at an altitude of 1000 m: 363 >>> temp2isa(45, 1000, temp_units = 'F', alt_units = 'm') 364 -2.2999999999999972 365 """ 366 367 # function tested in tests/test_std_atm.py 368 369 std_temp = alt2temp(altitude, alt_units, temp_units) 370 ISA_dev = temp - std_temp 371 372 return ISA_dev
373 374 375 # ############################################################################# 376 # 377 # Altitude to pressure and pressure ratio 378 # 379 # ############################################################################# 380 381
382 -def _alt2press_ratio_gradient( 383 H, 384 Hb, 385 Pb, 386 Tb, 387 L, 388 ):
389 390 # eqn from USAF TPS PEC binder, page PS1-31 391 392 return (Pb / P0) * (1 + (L / Tb) * (H - Hb)) ** ((-1000 * g) / (Rd 393 * L))
394 395
396 -def _alt2press_ratio_isothermal( 397 H, 398 Hb, 399 Pb, 400 Tb, 401 ):
402 403 # eqn from USAF TPS PEC binder, page PS1-26 404 405 return (Pb / P0) * M.exp((-1 * (H - Hb)) * ((1000 * g) / (Rd * Tb)))
406 407
408 -def alt2press_ratio(H, alt_units=default_alt_units):
409 """ 410 Return the pressure ratio (atmospheric pressure / standard pressure 411 for sea level). The altitude is specified in feet ('ft'), metres ('m'), 412 statute miles, ('sm') or nautical miles ('nm'). 413 414 If the units are not specified, the units in default_units.py are used. 415 416 Examples: 417 418 Calculate the pressure ratio at 5000 (default altitude units): 419 >>> alt2press_ratio(5000) 420 0.8320481158727735 421 422 Calculate the pressure ratio at 1000 m: 423 >>> alt2press_ratio(1000, alt_units = 'm') 424 0.88699304638887044 425 426 The functions are only implemented at altitudes of 84.852 km and lower. 427 >>> alt2press_ratio(90, alt_units = 'km') 428 Traceback (most recent call last): 429 File '<stdin>', line 1, in ? 430 File './std_atm.py', line 189, in alt2press_ratio 431 if H <= 20: 432 ValueError: This function is only implemented for altitudes of 84.852 km and below. 433 """ 434 435 # uses meters and degrees K for the internal calculations 436 437 # function tested in tests/test_std_atm.py 438 439 H = U.len_conv(H, from_units=alt_units, to_units='km') 440 441 if H <= 11: 442 return _alt2press_ratio_gradient(H, 0, P0, T0, L0) 443 if H <= 20: 444 return _alt2press_ratio_isothermal(H, 11, P11, T11) 445 if H <= 32: 446 return _alt2press_ratio_gradient(H, 20, P20, T20, L20) 447 if H <= 47: 448 return _alt2press_ratio_gradient(H, 32, P32, T32, L32) 449 if H <= 51: 450 return _alt2press_ratio_isothermal(H, 47, P47, T47) 451 if H <= 71: 452 return _alt2press_ratio_gradient(H, 51, P51, T51, L51) 453 if H <= 84.852: 454 return _alt2press_ratio_gradient(H, 71, P71, T71, L71) 455 else: 456 raise ValueError, \ 457 'This function is only implemented for altitudes of 84.852 km and below.'
458 459
460 -def alt2press(H, alt_units=default_alt_units, 461 press_units=default_press_units):
462 """ 463 Return the atmospheric pressure for a given altitude, with the 464 altitude in feet ('ft'), metres ('m'), statute miles, ('sm') or nautical 465 miles ('nm'), and the pressure in inches of HG ('in HG'), mm of HG 466 ('mm HG'), psi, lb per sq. ft ('psf'), pa, hpa or mb. 467 468 If the units are not specified, the units in default_units.py are used. 469 470 Examples: 471 472 Calculate the pressure in inches of mercury at 5,000 (default altitude 473 units): 474 >>> alt2press(5000) 475 24.895961289464015 476 477 Calculate the pressure in pounds per square foot at 10,000 (default 478 altitude units): 479 >>> alt2press(10000, press_units = 'psf') 480 1455.3301392981359 481 482 Calculate the pressure in pascal at 20 km: 483 >>> alt2press(20, press_units = 'pa', alt_units = 'km') 484 5474.8827144576408 485 """ 486 487 # uses meters, inches of HG and degrees K for the internal calculations 488 489 # function tested in tests/test_std_atm.py 490 491 H = U.len_conv(H, from_units=alt_units, to_units='m') 492 493 press = P0 * alt2press_ratio(H, alt_units='m') 494 press = U.press_conv(press, from_units='in HG', 495 to_units=press_units) 496 497 return press
498 499 500 # ############################################################################# 501 # 502 # Pressure altitude from barometric altitude and altimeter setting 503 # 504 # ############################################################################# 505 506
507 -def pressure_alt(H, alt_setting, alt_units=default_alt_units):
508 """ 509 Return the pressure altitude, given the barometric altitude and the 510 altimeter setting. 511 512 Altimeter setting may have units of inches of HG, or hpa or mb. If the 513 altimeter setting value is less than 35, the units are assumed to be 514 in HG, otherwise they are assumed to be hpa. The altimeter setting must 515 be in the range of 25 to 35 inches of mercury. 516 517 The altitude may have units of feet ('ft'), metres ('m'), statute miles, 518 ('sm') or nautical miles ('nm'). 519 520 If the units are not specified, the units in default_units.py are used. 521 522 Examples: 523 524 Calculate the pressure altitude for 1,000 (default altitude units) 525 barometric altitude with altimeter setting of 30.92 in HG: 526 >>> pressure_alt(1000, 30.92) 527 88.612734282205338 528 529 Calculate the pressure altitude for 1,000 (default altitude units) 530 barometric altitude with altimeter setting of 1008 mb: 531 >>> pressure_alt(1000, 1008) 532 1143.6503495627171 533 534 Calculate the pressure altitude in metres for 304.8 m barometric 535 altitude with altimeter setting of 1008 mb: 536 >>> pressure_alt(304.8, 1008, alt_units = 'm') 537 348.58462654671621 538 """ 539 540 H = U.len_conv(H, from_units=alt_units, to_units='ft') 541 if alt_setting > 35: 542 alt_setting = U.press_conv(alt_setting, from_units='hpa', 543 to_units='in HG') 544 if alt_setting < 25 or alt_setting > 35: 545 raise ValueError, 'Altimeter setting out of range.' 546 HP = H + 145442.2 * (1 - (alt_setting / P0) ** 0.190261) 547 HP = U.len_conv(HP, from_units='ft', to_units=alt_units) 548 return HP
549 550
551 -def QNH( 552 HP, 553 H, 554 alt_units=default_alt_units, 555 alt_setting_units='in HG', 556 ):
557 """ 558 Return the altimeter setting, given the pressure altitude (HP) and the 559 barometric altitude (H). 560 """ 561 562 HP = U.len_conv(HP, from_units=alt_units, to_units='ft') 563 H = U.len_conv(H, from_units=alt_units, to_units='ft') 564 QNH = P0 * (1 - (HP - H) / 145442.2) ** 5.255594 565 QNH = U.press_conv(QNH, from_units='in HG', 566 to_units=alt_setting_units) 567 568 return QNH
569 570 571 # ############################################################################# 572 # 573 # Altitude to density and density ratio 574 # 575 # ############################################################################# 576 577
578 -def alt2density_ratio(H, alt_units=default_alt_units):
579 """ 580 Return the density ratio (atmospheric density / standard density 581 for sea level). The altitude is specified in feet ('ft'), metres ('m'), 582 statute miles, ('sm') or nautical miles ('nm'). 583 584 If the units are not specified, the units in default_units.py are used. 585 586 Examples: 587 588 Calculate the density ratio at 7,500 (default altitude units): 589 >>> alt2density_ratio(7500) 590 0.79825819881753035 591 592 Calculate the density ratio at 2 km: 593 >>> alt2density_ratio(2, alt_units = 'km') 594 0.8216246960994622 595 """ 596 597 # function tested in tests/test_std_atm.py 598 599 return alt2press_ratio(H, alt_units) / alt2temp_ratio(H, alt_units)
600 601
602 -def alt2density(H, alt_units=default_alt_units, 603 density_units=default_density_units):
604 """ 605 Return the density given the pressure altitude. The altitude is 606 specified in feet ('ft'), metres ('m'), statute miles, ('sm') or 607 nautical miles ('nm'). 608 609 The desired density units are specified as 'lb/ft**3', 'slug/ft**3' or 610 'kg/m**3'. 611 612 If the units are not specified, the units in default_units.py are used. 613 614 Examples: 615 616 Calculate the density in lb / ft cubed at 7,500 (default altitude units): 617 >>> alt2density(7500) 618 0.061046199847730374 619 620 Calculate the density in slugs / ft cubed at 5,000 (default altitude units): 621 >>> alt2density(5000, density_units = 'slug/ft**3') 622 0.0020480982157718704 623 624 Calculate the density in kg / m cubed at 0 (default altitude units: 625 >>> alt2density(0, density_units = 'kg/m**3') 626 1.2250000000000001 627 628 Calculate the density in kg / m cubed at 81,000 m: 629 >>> alt2density(81000, density_units = 'kg/m**3', alt_units = 'm') 630 1.3320480184052337e-05 631 """ 632 633 # function tested in tests/test_std_atm.py 634 635 # get density in kg/m**3 636 637 density = Rho0 * alt2density_ratio(H, alt_units) 638 return U.density_conv(density, from_units='kg/m**3', 639 to_units=density_units)
640 641 642 # ############################################################################# 643 # 644 # Density to altitude and density ratio to altitude 645 # 646 # ############################################################################# 647 648
649 -def _density2alt_gradient( 650 Rho, 651 Rhob, 652 Hb, 653 Tb, 654 L, 655 ):
656 657 return Hb + (Tb / L) * ((Rho / Rhob) ** (-1 / ((1000 * g) / (Rd * L) 658 + 1)) - 1)
659 660
661 -def _density2alt_isothermal( 662 Rho, 663 Rhob, 664 Hb, 665 Tb, 666 ):
667 668 return Hb - ((Rd * Tb) * M.log(Rho / Rhob)) / (1000 * g)
669 670
671 -def density2alt(Rho, density_units=default_density_units, 672 alt_units=default_alt_units):
673 """ 674 Return the altitude corresponding to the specified density, with 675 density in 'lb/ft**3', 'slug/ft**3' or 'kg/m**3'. 676 677 The altitude is specified in feet ('ft'), metres ('m'), statute miles, 678 ('sm') or nautical miles ('nm'). 679 680 If the units are not specified, the units in default_units.py are used. 681 682 Examples: 683 684 Calculate the altitude in default altitude units where the density is 685 0.056475 in default density units: 686 >>> density2alt(.056475) 687 9999.8040934937271 688 689 Calculate the altitude in metres where the density is 0.018012 kg / m 690 cubed: 691 >>> density2alt(.018012, alt_units = 'm', density_units = 'kg/m**3') 692 29999.978688508152 693 """ 694 695 # function tested in tests/test_std_atm.py 696 697 Rho = U.density_conv(Rho, from_units=density_units, 698 to_units='kg/m**3') 699 700 if Rho > Rho11: 701 H = _density2alt_gradient(Rho, Rho0, 0, T0, L0) 702 elif Rho > Rho20: 703 H = _density2alt_isothermal(Rho, Rho11, 11, T11) 704 elif Rho > Rho32: 705 H = _density2alt_gradient(Rho, Rho20, 20, T20, L20) 706 elif Rho > Rho47: 707 H = _density2alt_gradient(Rho, Rho32, 32, T32, L32) 708 elif Rho > Rho51: 709 H = _density2alt_isothermal(Rho, Rho47, 47, T47) 710 elif Rho > Rho71: 711 H = _density2alt_gradient(Rho, Rho51, 51, T51, L51) 712 else: 713 H = _density2alt_gradient(Rho, Rho71, 71, T71, L71) 714 715 if H > 84.852: 716 raise ValueError, \ 717 'This function is only implemented for altitudes of 84.852 km and below.' 718 719 return U.len_conv(H, from_units='km', to_units=alt_units)
720 721
722 -def density_ratio2alt(DR, alt_units=default_alt_units):
723 """ 724 Return the altitude for the specified density ratio. The altitude is in 725 feet ('ft'), metres ('m'), statute miles, ('sm') or nautical miles 726 ('nm'). 727 728 If the units are not specified, the units in default_units.py are used. 729 730 Examples: 731 732 Calculate the altitude in default altitude units where the density ratio is 733 1: 734 >>> density_ratio2alt(1) 735 0.0 736 737 Calculate the altitude in feet where the density ratio is 0.5: 738 >>> density_ratio2alt(.5) 739 21859.50324995652 740 741 Calculate the altitude in km where the density ratio is 0.1 742 >>> density_ratio2alt(.1, alt_units = 'km') 743 17.9048674520646 744 """ 745 746 # function tested in tests/test_std_atm.py 747 748 D = DR * Rho0 749 return density2alt(D, alt_units=alt_units, density_units='kg/m**3')
750 751 752 # ############################################################################# 753 # 754 # Density Altitude 755 # 756 # ############################################################################# 757 758
759 -def density_alt( 760 H, 761 T, 762 alt_setting=P0, 763 DP='FALSE', 764 RH=0.0, 765 alt_units=default_alt_units, 766 temp_units=default_temp_units, 767 ):
768 """ 769 Return density altitude, given the pressure altitude and the 770 temperature with altitudes in units of feet ('ft'), metres ('m'), 771 statute miles, ('sm') or nautical miles ('nm'), and temperature in units 772 of deg C, F, K or R ('C', 'F', 'K' or 'R'). 773 774 Mandatory parametres: 775 H = altitude 776 T = temperature 777 778 Optional parametres: 779 alt_setting = altimeter setting (defaults to 29.9213 if not provided 780 DP = dew point 781 RH = relative humidity 782 alt_units = units for the altitude. 'ft', 'm', or 'km'. 783 temp_units = units for the temperature and dew point. 'C', 'F', 'K' 784 or 'R'. 785 786 The altimeter setting units are assumed to be inches of HG, unless the 787 value is greater than 35. In this case the units are assumed to be mb. 788 789 If the dew point or relative humidity are not specified, the air is 790 assumed to be completely dry. If both the dew point and relative humidity 791 are specified, the relative humidity value is ignored. 792 793 If the units are not specified, the units in default_units.py are used. 794 795 The method is from: http://wahiduddin.net/calc/density_altitude.htm 796 797 Examples: 798 799 Calculate the density altitude in default altitude units for a pressure 800 altitude of 7000 default altitude units and a temperature of 15 deg 801 (default temperature units). The altimeter setting is not specified, so it 802 defaults to standard pressure of 29.9213 in HG or 1013.25 mb: 803 >>> density_alt(7000, 15) 804 8595.3465863232504 805 806 Calculate the density altitude in default altitude units for a pressure 807 altitude of 7000 default altitude units and a temperature of 85 deg F. 808 The altimeter setting is not specified, so it defaults to standard pressure 809 of 29.9213 in HG or 1013.25 mb. The dew point and relative humidity are 810 not specified, so the air is assumed to be dry: 811 >>> density_alt(7000, 85, temp_units = 'F') 812 10159.10696106757 813 814 Calculate the density altitude in default altitude units for a pressure 815 altitude of 7000 default altitude units, an altimeter setting of 29.80 and 816 a temperature of 85 deg F and a dew point of 55 deg F: 817 >>> density_alt(7000, 85, 29.80, 55, temp_units = 'F') 818 10522.776013011618 819 820 Calculate the density altitude in metres for a pressure altitude of 821 2000 m, an altimeter setting of 1010 mb, a temperature of 15 deg (default 822 temperature units) and a relative humidity of 50%: 823 >>> density_alt(2000, 15, 1010, alt_units = 'm', RH = 0.5) 824 2529.8230634449737 825 826 The dew point may be specified in one of two ways: as the fourth 827 argument on the command line, or via the keyword argument DP. 828 >>> density_alt(2000, 15, 1010, alt_units = 'm', DP = 5) 829 2530.7528237990618 830 831 The relative humidity must be in the range of 0 to 1: 832 >>> density_alt(2000, 15, 1010, alt_units = 'm', RH = 1.1) 833 Traceback (most recent call last): 834 File '<stdin>', line 1, in ? 835 File 'std_atm.py', line 533, in density_alt 836 raise ValueError, 'The relative humidity must be in the range of 0 to 1.' 837 ValueError: The relative humidity must be in the range of 0 to 1. 838 """ 839 840 Rv = 461.495 # gas constant for water vapour 841 842 # saturated vapour pressure 843 844 if DP == 'FALSE' and RH == 0: 845 Pv = 0 846 else: 847 Pv = sat_press(T, DP, RH, temp_units, press_units='pa') 848 849 # dry air pressure 850 851 Pd = dry_press(H, Pv, alt_setting=alt_setting, alt_units=alt_units, 852 press_units='pa') 853 854 T = U.temp_conv(T, from_units=temp_units, to_units='K') 855 D = Pd / (Rd * T) + Pv / (Rv * T) 856 857 DR = D / Rho0 858 859 return density_ratio2alt(DR, alt_units)
860 861
862 -def _sat_press(T):
863 """ 864 Return the saturation pressure in mb of the water vapour, given 865 temperature in deg C. Equation from: 866 http://wahiduddin.net/calc/density_altitude.htm 867 """ 868 869 eso = 6.1078 870 c0 = 0.99999683 871 c1 = -0.90826951e-2 872 c2 = 0.78736169e-4 873 c3 = -0.61117958e-6 874 c4 = 0.43884187e-8 875 c5 = -0.29883885e-10 876 c6 = 0.21874425e-12 877 c7 = -0.17892321e-14 878 c8 = 0.11112018e-16 879 c9 = -0.30994571e-19 880 881 p = c0 + T * (c1 + T * (c2 + T * (c3 + T * (c4 + T * (c5 + T * (c6 882 + T * (c7 + T * (c8 + T * c9)))))))) 883 sat_press = eso / p ** 8 884 return sat_press
885 886
887 -def sat_press( 888 T='FALSE', 889 DP='FALSE', 890 RH=0.0, 891 temp_units=default_temp_units, 892 press_units=default_press_units, 893 ):
894 """ 895 Return the saturated vapour pressure of water. Either the dew point, or 896 the temperature and the relative humidity must be specified. If both the 897 dew point and relative humidity are specified, the relative humidity value 898 is ignored. 899 900 If the temperature and dew point are both specified, the dew point cannot 901 be greater than the temperature: 902 903 If the units are not specified, the units in default_units.py are used. 904 905 >>> sat_press(T=10, DP=11) 906 Traceback (most recent call last): 907 File '<stdin>', line 1, in <module> 908 File 'std_atm.py', line 795, in sat_press 909 raise ValueError, 'The dew point cannot be greater than the temperature.' 910 ValueError: The dew point cannot be greater than the temperature. 911 912 Dew point is 11 deg (default temperature units). Find the water vapour 913 pressure in default pressure units: 914 >>> sat_press(DP=11) 915 0.38741015927568667 916 917 Dew point is 65 deg F. Find the water vapour pressure in default pressure units: 918 >>> sat_press(DP=65, temp_units = 'F') 919 0.62207710701956165 920 921 Dew point is 212 deg F (the boiling point of water at sea level). 922 Find the water vapour pressure in lb per sq. inch: 923 >>> sat_press(DP=212, temp_units = 'F', press_units = 'psi') 924 14.696764873564959 925 926 Temperature is 30 deg C. Find the water vapour pressure in default pressure units: 927 for 50% relative humidity: 928 >>> sat_press(T=30, RH = 0.5) 929 0.62647666996057927 930 """ 931 932 if DP != 'FALSE': 933 934 # use dew point method 935 936 if T != 'FALSE': 937 if DP > T: 938 raise ValueError, \ 939 'The dew point cannot be greater than the temperature.' 940 941 DP = U.temp_conv(DP, from_units=temp_units, to_units='C') 942 943 # calculate vapour pressure 944 945 Pv = _sat_press(DP) * 100 946 else: 947 948 if RH == 'FALSE': 949 raise ValueError, \ 950 'Either DP (dew point) or RH (relative humidity) must be specified.' 951 952 # relative humidity is specified 953 # confirm relative humidity is in range 954 955 if RH < 0 or RH > 1: 956 raise ValueError, \ 957 'The relative humidity must be in the range of 0 to 1.' 958 959 if T == 'FALSE': 960 raise ValueError, \ 961 'If the relative humidity is specified, the temperature must also be specified.' 962 963 T = U.temp_conv(T, from_units=temp_units, to_units='C') 964 965 Pv = _sat_press(T) * 100 966 Pv *= RH 967 968 Pv = U.press_conv(Pv, from_units='pa', to_units=press_units) 969 970 return Pv
971 972
973 -def dry_press( 974 H, 975 Pv, 976 alt_setting=P0, 977 alt_units=default_alt_units, 978 press_units=default_press_units, 979 ):
980 """ 981 Returns dry air pressure, i.e. the total air pressure, less the water 982 vapour pressure. 983 """ 984 985 HP = pressure_alt(H, alt_setting, alt_units=alt_units) 986 P = alt2press(HP, press_units=press_units, alt_units=alt_units) 987 Pd = P - Pv 988 989 return Pd
990 991
992 -def density_alt2temp( 993 density_alt_seek, 994 press_alt, 995 alt_units=default_alt_units, 996 temp_units=default_temp_units, 997 ):
998 """ 999 Return temperature to achieve a desired density altitude. 1000 1001 If the units are not specified, the units in default_units.py are used. 1002 """ 1003 1004 low = -100 # initial lower guess 1005 high = 100 # initial upper guess 1006 1007 # confirm initial low and high are OK: 1008 1009 da_low = density_alt(press_alt, low, alt_units=alt_units) 1010 if da_low > density_alt_seek: 1011 raise ValueError, 'Initial low guess too high.' 1012 1013 da_high = density_alt(press_alt, high, alt_units=alt_units) 1014 if da_high < density_alt_seek: 1015 raise ValueError, 'Initial high guess too low.' 1016 1017 guess = (low + high) / 2. 1018 da_guess = density_alt(press_alt, guess, alt_units=alt_units) 1019 1020 # keep iterating until da is within 1 ft of desired value 1021 1022 while M.fabs(da_guess - density_alt_seek) > 1: 1023 if da_guess > density_alt_seek: 1024 high = guess 1025 else: 1026 low = guess 1027 1028 guess = (low + high) / 2. 1029 da_guess = density_alt(press_alt, guess, alt_units=alt_units) 1030 1031 guess = U.temp_conv(guess, from_units='C', to_units=temp_units) 1032 1033 return guess
1034 1035
1036 -def density_alt_table( 1037 density_alt_seek, 1038 alt_range=2000, 1039 alt_inc=100, 1040 alt_units=default_alt_units, 1041 temp_units=default_temp_units, 1042 multi_units=False, 1043 file='', 1044 format='text', 1045 ):
1046 """ 1047 Return a text or html table of required temperature vs pressure altitude. 1048 1049 If the units are not specified, the units in default_units.py are used. 1050 """ 1051 1052 line_buffer = [] 1053 if format == 'text': 1054 line_buffer.append('Pressure altitudes and temperatures for a density ' 1055 ) 1056 line_buffer.append('altitude of ' + str(density_alt_seek) + ' ' 1057 + alt_units) 1058 line_buffer.append('(assuming dry air)\n') 1059 if multi_units: 1060 line_buffer.append(' Pressure Temp Temp') 1061 line_buffer.append(' Altitude') 1062 line_buffer.append(' (' + alt_units 1063 + ') (deg C) (deg F)') 1064 else: 1065 line_buffer.append(' Pressure Temp') 1066 line_buffer.append(' Altitude') 1067 line_buffer.append(' (' + alt_units + ') (deg ' 1068 + temp_units + ')') 1069 elif format == 'html': 1070 print 'creating html' 1071 else: 1072 raise ValueError, \ 1073 'Invalid format. Must be either "text" or "html"' 1074 1075 if multi_units: 1076 for alt in range(max(density_alt_seek - alt_range / 2., 0), 1077 density_alt_seek + alt_range / 2. + alt_inc, 1078 alt_inc): 1079 temp_c = density_alt2temp(density_alt_seek, alt, 1080 alt_units=alt_units) 1081 temp_f = U.temp_conv(temp_c, from_units='C', to_units='F') 1082 alt_str = L.format('%.*f', (0, alt), grouping=True) 1083 temp_c_str = '%.1f' % temp_c 1084 temp_f_str = '%.1f' % temp_f 1085 line_buffer.append(alt_str.rjust(6) + temp_c_str.rjust(11) 1086 + temp_f_str.rjust(10)) 1087 else: 1088 for alt in range(max(density_alt_seek - alt_range / 2., 0), 1089 density_alt_seek + alt_range / 2. + alt_inc, 1090 alt_inc): 1091 alt_str = L.format('%.*f', (0, alt), grouping=True) 1092 temp_str = '%.1f' % density_alt2temp(density_alt_seek, alt, 1093 temp_units=temp_units, alt_units=alt_units) 1094 line_buffer.append(alt_str.rjust(6) + temp_str.rjust(11)) 1095 1096 if file != '': 1097 OUT = open(file, 'w') 1098 for line in line_buffer: 1099 OUT.write(line + '\n') 1100 1101 print 'file selected' 1102 else: 1103 return '\n'.join(line_buffer)
1104 1105 1106 # ############################################################################# 1107 # 1108 # Pressure to altitude and pressure ratio to altitude 1109 # 1110 # ############################################################################# 1111 1112
1113 -def _press2alt_gradient( 1114 P, 1115 Pb, 1116 Hb, 1117 Tb, 1118 L, 1119 ):
1120 1121 return Hb + (Tb / L) * ((P / Pb) ** (((-1 * Rd) * L) / (1000 * g)) 1122 - 1)
1123 1124
1125 -def _press2alt_isothermal( 1126 P, 1127 Pb, 1128 Hb, 1129 Tb, 1130 ):
1131 1132 return Hb - ((Rd * Tb) * M.log(P / Pb)) / (1000 * g)
1133 1134
1135 -def press2alt(P, press_units=default_press_units, 1136 alt_units=default_alt_units):
1137 """ 1138 Return the altitude corresponding to the specified pressure, with 1139 pressure in inches of HG, mm of HG, psi, psf (lb per sq. ft), pa, hpa or 1140 mb. 1141 1142 The altitude is in units of feet ('ft'), metres ('m'), statute miles, 1143 ('sm') or nautical miles ('nm') 1144 1145 If the units are not specified, the units in default_units.py are used. 1146 1147 Examples: 1148 1149 Calculate the pressure altitude in feet for a pressure of 31.0185 inches 1150 of HG: 1151 >>> press2alt(31.0185) 1152 -999.98992888235091 1153 1154 Calculate the pressure altitude in feet for a pressure of 1155 1455.33 lb sq. ft: 1156 >>> press2alt(1455.33, press_units = 'psf') 1157 10000.002466564831 1158 1159 Calculate the pressure altitude in metres for a pressure of 1160 90.3415 mm HG: 1161 >>> press2alt(90.3415, press_units = 'mm HG', alt_units = 'm') 1162 15000.025465320754 1163 1164 Calculate the pressure altitude in metres for a pressure of 1165 1171.86 pascal: 1166 >>> press2alt(1171.86, press_units = 'pa', alt_units = 'm') 1167 30000.029510365184 1168 """ 1169 1170 # function tested in tests/test_std_atm.py 1171 1172 P = U.press_conv(P, from_units=press_units, to_units='in HG') 1173 1174 if P > P11: 1175 H = _press2alt_gradient(P, P0, 0, T0, L0) 1176 elif P > P20: 1177 H = _press2alt_isothermal(P, P11, 11, T11) 1178 elif P > P32: 1179 H = _press2alt_gradient(P, P20, 20, T20, L20) 1180 elif P > P47: 1181 H = _press2alt_gradient(P, P32, 32, T32, L32) 1182 elif P > P51: 1183 H = _press2alt_isothermal(P, P47, 47, T47) 1184 elif P > P71: 1185 H = _press2alt_gradient(P, P51, 51, T51, L51) 1186 else: 1187 H = _press2alt_gradient(P, P71, 71, T71, L71) 1188 1189 if H > 84.852: 1190 raise ValueError, \ 1191 'This function is only implemented for altitudes of 84.852 km and below.' 1192 1193 return U.len_conv(H, from_units='km', to_units=alt_units)
1194 1195
1196 -def press_ratio2alt(PR, alt_units=default_alt_units):
1197 """ 1198 Return the pressure ratio for the specified altitude. The altitude is 1199 specified in feet ('ft'), metres ('m'), statute miles, ('sm') or 1200 nautical miles ('nm'). 1201 1202 If the units are not specified, the units in default_units.py are used. 1203 1204 Examples: 1205 1206 Calculate the altitude in feet where the pressure ratio is 0.5: 1207 >>> press_ratio2alt(.5) 1208 17969.990746028907 1209 1210 Calculate the altitude in metres where the pressure ratio is 0.1: 1211 >>> press_ratio2alt(.1, alt_units = 'm') 1212 16096.249927559489 1213 """ 1214 1215 # function tested in tests/test_std_atm.py 1216 1217 P = PR * P0 1218 return press2alt(P, alt_units=alt_units)
1219 1220 1221 # ############################################################################# 1222 # 1223 # Temperature to speed of sound 1224 # 1225 # ############################################################################# 1226 1227
1228 -def temp2speed_of_sound(temp, temp_units=default_temp_units, 1229 speed_units=default_speed_units):
1230 """ 1231 Return the speed of sound, given the air temperature. 1232 1233 The temperature units may be deg C, F, K or R ('C', 'F', 'K' or 'R'). 1234 1235 The speed units may be 'kt', 'mph', 'km/h', 'm/s' and 'ft/s'. 1236 1237 If the units are not specified, the units in default_units.py are used. 1238 1239 Examples: 1240 1241 Determine speed of sound in knots at 15 deg (default temperature units): 1242 >>> temp2speed_of_sound(15) 1243 661.47882487301808 1244 1245 Determine speed of sound in mph at 120 deg F: 1246 >>> temp2speed_of_sound(120, speed_units = 'mph', temp_units = 'F') 1247 804.73500154991291 1248 """ 1249 1250 # function tested in tests/test_std_atm.py 1251 1252 temp = U.temp_conv(temp, from_units=temp_units, to_units='K') 1253 1254 speed_of_sound = M.sqrt((1.4 * Rd) * temp) 1255 speed_of_sound = U.speed_conv(speed_of_sound, from_units='m/s', 1256 to_units=speed_units) 1257 1258 return speed_of_sound
1259 1260 1261 if __name__ == '__main__': 1262 1263 # run doctest to check the validity of the examples in the doc strings. 1264 1265 import doctest 1266 import sys 1267 doctest.testmod(sys.modules[__name__]) 1268