' demo_us76.bas

' this program demonstrates how to interact with the
' us76.bas subroutine which computes the properties
' of the U.S. 1976 standard atmosphere between zero
' and 1000 kilometers altitude

''''''''''''''''''''''''''''''''''''''

' dimension atmospheric data arrays

dim zs(49) as float, ps(49) as float, tms(49) as float, wms(49) as float

' read atmospheric data

for i% = 1 to 49

read zs(i%)

next i%

data 0.0, 11.019, 20.063, 32.162, 47.35, 51.413
data 71.802, 86.0, 91.0, 94.0, 97.0, 100.0, 103.0
data 106.0, 108.0, 110.0, 112.0, 115.0, 120.0, 125.0
data 130.0, 135.0, 140.0, 145.0, 150.0, 155.0, 160.0
data 165.0, 170.0, 180.0, 190.0, 210.0, 230.0, 265.0
data 300.0, 350.0, 400.0, 450.0, 500.0, 550.0, 600.0
data 650.0, 700.0, 750.0, 800.0, 850.0, 900.0, 950.0, 1000.0

for i% = 1 to 49

read tms(i%)

next i%

data 288.15, 216.65, 216.65, 228.65, 270.65
data 270.65, 214.65, 186.95, 186.87, 187.74
data 190.40, 195.08, 202.23, 212.89, 223.29
data 240.00, 264.00, 300.00, 360.00, 417.23
data 469.27, 516.59, 559.63, 598.78, 634.39
data 666.80, 696.29, 723.13, 747.57, 790.07
data 825.31, 878.84, 915.78, 955.20, 976.01
data 990.06, 995.83, 998.22, 999.24, 999.67
data 999.85, 999.93, 999.97, 999.99, 999.99
data 1000.00, 1000.00, 1000.00, 1000.00

for i% = 1 to 49

read wms(i%)

next i%

data 28.9644, 28.9644, 28.9644, 28.9644, 28.9644
data 28.9644, 28.9644, 28.9522, 28.889, 28.783
data 28.620, 28.395, 28.104, 27.765, 27.521
data 27.268, 27.020, 26.680, 26.205, 25.803
data 25.436, 25.087, 24.749, 24.422, 24.103
data 23.792, 23.488, 23.192, 22.902, 22.342
data 21.809, 20.825, 19.952, 18.688, 17.726
data 16.735, 15.984, 15.247, 14.330, 13.092
data 11.505, 9.718, 7.998, 6.579, 5.543
data 4.849, 4.404, 4.122, 3.940

for i% = 1 to 49

read ps(i%)

next i%

data 1013.25, 226.32, 54.7487, 8.68014, 1.10905
data 0.66938, 0.039564, 3.7338e-3, 1.5381e-3
data 9.0560e-4, 5.3571e-4, 3.2011e-4, 1.9742e-4
data 1.2454e-4, 9.3188e-5, 7.1042e-5, 5.5547e-5
data 4.0096e-5, 2.5382e-5, 1.7354e-5, 1.25054e-5
data 9.3568e-6, 7.2028e-6, 5.6691e-6, 4.5422e-6
data 3.6930e-6, 3.0395e-6, 2.5278e-6, 2.1210e-6
data 1.5271e-6, 1.1266e-6, 6.4756e-7, 3.9276e-7
data 1.7874e-7, 8.7704e-8, 3.4498e-8, 1.4518e-8
data 6.4468e-9, 3.0236e-9, 1.5137e-9, 8.2130e-10
data 4.8865e-10, 3.1908e-10, 2.2599e-10, 1.7036e-10
data 1.3415e-10, 1.0873e-10, 8.9816e-11, 7.5138e-11

' define altitude and compute properties

alt = 100.0

us76(alt, atmtmp, atmdns, atmprs, atmsos)

' display results

print " "

print "properties of the U.S. 1976 standard atmosphere"

print "-----------------------------------------------"

print " "

print "altitude (kilometers) = ", alt

print " "

print "temperature (degrees Kelvin) = ", atmtmp

print " "

print "density (kilograms/cubic meter) = ", atmdns

print " "

print "pressure (newton/square meter) = ", atmprs

print " "

print "speed of sound (meters/second) = ", atmsos

print " "

end

'''''''''''''''''''''''''''''''''''''''''''''
'''''''''''''''''''''''''''''''''''''''''''''

sub us76(alt, atmtmp, atmdns, atmprs, atmsos)

' U.S. 1976 standard atmosphere subroutine

' input

' alt = altitude (kilometers)

' output

' atmtmp = temperature (degrees Kelvin)
' atmdns = density (kilograms/cubic meter)
' atmprs = pressure (newton/square meter)
' atmsos = speed of sound (meters/second)

''''''''''''''''''''''''''''''''''''''''''

local i%, j%

z = alt

' altitude must be less than 1000 kilometers

if (z > 1000.0) then

z = 999.0

end if

' determine positive altitude values

if (z >= 0.0) then

' set initial values

ro = 6356.766

go = 9.80665

wmo = 28.9644

rs = 8314.32

zm = z * 1000.0

rom = ro * 1000.0

' ------------------------------------
' altitudes greater than 86 kilometers
' ------------------------------------

if (z >= 86.0) then

' find pointer location

for i% = 8 to 48

if (zs(i%) <= z and z < zs(i% + 1)) then

exit for

end if

next i%

' calculate temperature

if (i% = 8) then

t = tms(9)

elseif (not (i% < 16 or i% >= 19)) then

t = 240.0 + 12.0 * (z - 110.0)

elseif (i% >= 19) then

xi = (z - 120.0) * (ro + 120.0) / (ro + z)

t = 1000.0 - 640.0 * exp(-0.01875 * xi)

else

t = 263.1905 - 76.3232 * sqr(1.0 - ((z - 91.0) / 19.9429)^2)

end if

' set pointers

j% = i%

if (i% = 48) then

j% = i% - 1

end if

' calculate altitude indices values for molecular weight calculation

z0 = zs(j%)

z1 = zs(j% + 1)

z2 = zs(j% + 2)

' calculate molecular weight

wma = wms(j%) * (z - z1) * (z - z2) / ((z0 - z1) * (z0 - z2))

wma = wma + wms(j% + 1) * (z - z0) * (z - z2) / ((z1 - z0) * (z1 - z2))

wma = wma + wms(j% + 2) * (z - z0) * (z - z1) / ((z2 - z0) * (z2 - z1))

' calculate atmospheric pressure indices

alp0 = log(ps(j%))

alp1 = log(ps(j% + 1))

alp2 = log(ps(j% + 2))

' calculate pressure

alpa = alp0 * (z - z1) * (z - z2) / ((z0 - z1) * (z0 - z2))

alpa = alpa + alp1 * (z - z0) * (z - z2) / ((z1 - z0) * (z1 - z2))

alpa = alpa + alp2 * (z - z0) * (z - z1) / ((z2 - z0) * (z2 - z1))

' set second side of molecular weight and pressure value

alpb = alpa

wmb = wma

' perform same calculations for far side values for smoothing

if (not (i% = 8 or i% = 48)) then

j% = j% - 1

' calculate altitude indices values for molecular weight calculation

z0 = zs(j%)

z1 = zs(j% + 1)

z2 = zs(j% + 2)

' calculate pressure indices values for pressure calculation

alp0 = log(ps(j%))

alp1 = log(ps(j% + 1))

alp2 = log(ps(j% + 2))

' calculate atmospheric pressure

alpb = alp0 * (z - z1) * (z - z2) / ((z0 - z1) * (z0 - z2))

alpb = alpb + alp1 * (z - z0) * (z - z2) / ((z1 - z0) * (z1 - z2))

alpb = alpb + alp2 * (z - z0) * (z - z1) / ((z2 - z0) * (z2 - z1))

' calculate molecular weight

wmb = wms(j%) * (z - z1) * (z - z2) / ((z0 - z1) * (z0 - z2))

wmb = wmb + wms(j% + 1) * (z - z0) * (z - z2) / ((z1 - z0) * (z1 - z2))

wmb = wmb + wms(j% + 2) * (z - z0) * (z - z1) / ((z2 - z0) * (z2 - z1))

end if

' calculate pressure and molecular weight from derived constants

p = 100.0 * exp((alpa + alpb) / 2.0)

wm = (wma + wmb) / 2.0

else

' find pointer location

for i% = 1 to 7

if (zs(i%) <= z and z < zs(i% + 1)) then

exit for

end if

next i%

' calculate altitude breakpoints

zl = ro * zs(i%) / (ro + zs(i%))

zu = ro * zs(i% + 1) / (ro + zs(i% + 1))

' convert to meters

zlm = zl * 1000.0

zum = zu * 1000.0

wm = wmo

ht = (ro * z) / (ro + z)

hm = ht * 1000.0

g = (tms(i% + 1) - tms(i%)) / (zu - zl)

gm = g * 0.001

' calculate pressure and temperature

if (g < 0.0 or g > 0.0) then

p = ps(i%) * ((tms(i%) / (tms(i%) + g * (ht - zl)))^((go * wmo) / (rs * gm))) * 100.0

else

p = ps(i%) * exp(-(go * wmo * (hm - zlm)) / (rs * tms(i%))) * 100.0

end if

t = tms(i%) + g * (ht - zl)

end if

' calculate density

d = (wm * p) / (rs * t)

else

' set default values for negative altitudes

t = 288.15

p = 101325.0

d = 1.224999155888

end if

' calculate speed of sound

if (z < 91.0) then

ss = 20.046707 * sqr(t)

else

ss = 274.0392912

end if

atmdns = d

atmtmp = t

atmsos = ss

atmprs = p

end sub