<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 ISA
Step1:  Get layer data
Step2:  Get $p, T, \rho$
Step3:  Note that the geopotential height can be provided as a numpy array
Step4:  If height is provided as a
Step5:  Height computation
The geopotential height can be recalculated given a temperature or pressure, by using isa.height_from_temperature_isa(T) or  isa.height_from_pressure_isa.
Note that different heights may produce the same temperature. For example
Step6:  On the other hand, pressure-altitude yields a single value of the geopotential height
Step7:  Array-numpy can be used to input more than one temperature or pressure, wich yields a geopoential height numpy array. Note that if more than one height value can be retrieved, then an array of arrays is returned. I.e
Step8:  Finally, suppose $h=7500m, \Delta T_{ISA}=0K$ and $h=9000m, \Delta T_{ISA}=+10K$
Step9:  Now retrieve geopotential height, indicating the ISA deviation
<ASSISTANT_TASK:>
Python Code:
# Import isa library
from pyturb.gas_models import isa
import numpy as np
from matplotlib import pyplot as plt
Explanation: ISA: Example of usage
This Notebook serves as an example on how to use the Standard Atmosphere model provided with pyTurb. The isa functions can be found in the gas_models folder:
pyturb
gas_models
isa:
python
    from pyturb.gas_models import isa
The International Standard Atmosphere stablishes a model for the temperature, 
pressure, density and Mach number as a function of the geopotential heights.
The ISA model refers to the US Atmosphere model of 1975 (COESA 1975), although 
from 50 km on, ISA model stablishes the 1975 US Atmosphere Model as interim.
This notebook is intended to serve as use-case for pressure, temperature and density as a function of the geopotential height, and for geopontial height as a function of pressure and temperature.
Imports:
pyTurb libraries:
- pyturb.gas_models.isa
Others (plotting results):
- numpy
- matplotlib.pyplot
End of explanation
height = [0, 11000, 20000, 32000, 47000, 51000, 71000, 84852]
for i_layer, h in enumerate(height):
    lapse_rate, Tbase, pbase, dbase, heightbase, layer_name = isa.get_atmosdata(h)
    print('{5} -{4:>14s} - zbase={3:9.3f}m, alpha={0:9.2e}K/m, Tbase={1:6.2f}K, pbase={2:10.3f}Pa, dbase={6:6.3e}kg/m^3'.format(lapse_rate, Tbase, pbase, heightbase, layer_name, i_layer+1, dbase))
Explanation: Get layer data:
The layer information of a given h geopotential height ($h [m]$) can be retrived with isa.get_atmosdata(h):
End of explanation
# For geopotential height of 1000m
T0 = isa.temperature_isa(1000)
p0 = isa.pressure_isa(1000)
rho0 = isa.density_state_eq(1000)
print('T=', T0, 'K; p=', p0, 'Pa; rho=', rho0, 'kg/m^3')
Explanation: Get $p, T, \rho$:
Given a geopotential height in meters, $h [m]$, the pressure, temperature and density can be calculated as:
- T = isa.temperature_isa(h, isa_dev)
- p = isa.pressure_isa(h, isa_dev)
- $\rho$ = isa.density_isa(h, isa_dev)
The density can algo be calculated with the state equation, provided a temperature and pressure:
$$\rho =\frac{p}{R_{air}T}$$
- $\rho$ = isa.density_state_eq(h, isa_dev)
Note that isa_dev is the ISA temperature deviation, assuming 15ºC for a standard day. By default, isa_dev is zero.
End of explanation
# Create a vector from troposphere to mesosphere3 (0m to 84852m).
h = np.linspace(0,84851, 100)
# Temperature:
T = isa.temperature_isa(h)
# Pressure
p = isa.pressure_isa(h)
# Density
d = isa.density_state_eq(h)
plt.figure()
plot_temp = plt.plot(h,T, '-x')
plt.ylabel("Temperature [K]")
plt.xlabel("geopotential height [m]")
fig, eje1 = plt.subplots()
eje2 = eje1.twinx()
plot1, = eje1.plot(h, p, "b-+", label="pressure")
plot2, = eje2.plot(h, d, "r-x", label="density_steq")
eje1.set_xlabel("geopotential height [m]")
eje1.set_ylabel("pressure [Pa]")
eje2.set_ylabel("density kg/m3")
eje1.yaxis.label.set_color(plot1.get_color())
eje2.yaxis.label.set_color(plot2.get_color())
tkw = dict(size=4, width=1.5)
eje1.tick_params(axis='y', colors=plot1.get_color(), **tkw)
eje2.tick_params(axis='y', colors=plot2.get_color(), **tkw)
eje1.tick_params(axis='x', **tkw)
lines = [plot1, plot2]
eje1.legend(lines, [l.get_label() for l in lines])
Explanation: Note that the geopotential height can be provided as a numpy array:
End of explanation
print('Temperature:')
print('isa_dev by default: T=', isa.temperature_isa([0, 1000, 2000, 3000]))
print('isa_dev=10, T=', isa.temperature_isa([0, 1000, 2000, 3000], 10))
print('isa_dev as array, T=', isa.temperature_isa([0, 1000, 2000, 3000], [5, 10, 15, 20]))
print('Pressure:')
print('isa_dev by default: ', isa.pressure_isa([5000, 8000, 10000, 0]))
print('isa_dev as array, p=', isa.pressure_isa([5000, 8000, 10000, 0], [5, 10, 15, 20]))
Explanation: If height is provided as a
End of explanation
isa.height_from_temperature_isa(250)
plt.figure()
plot_temp = plt.plot(h,T, '-x')
plt.plot([0,84000], [250, 250])
plt.ylabel("Temperature [K]")
plt.xlabel("geopotential height [m]")
Explanation: Height computation
The geopotential height can be recalculated given a temperature or pressure, by using isa.height_from_temperature_isa(T) or  isa.height_from_pressure_isa.
Note that different heights may produce the same temperature. For example:
$$T=250K \rightarrow h=\left[5869.2m, 39625.0m, 58375.0m\right]$$
End of explanation
plt.figure()
plot_temp = plt.plot(h,p, '-x')
plt.ylabel("Pressure [Pa]")
plt.xlabel("geopotential height [m]")
isa.height_from_pressure_isa([90000, 50000, 2500])
Explanation: On the other hand, pressure-altitude yields a single value of the geopotential height:
End of explanation
isa.height_from_temperature_isa([200, 220, 250, 275])
Explanation: Array-numpy can be used to input more than one temperature or pressure, wich yields a geopoential height numpy array. Note that if more than one height value can be retrieved, then an array of arrays is returned. I.e:
- $ T=200K \rightarrow h=78325m$
- $ T=220K \rightarrow \left[h=10484.6, 23350, 69089.3\right]$
- $T=250K \rightarrow \left[h=5869.2m, 39625.0m, 58375m\right]$
- $ T=275K \rightarrow h=2023.1m$
End of explanation
T75 = isa.temperature_isa(7500)
print(T75)
T90 = isa.temperature_isa(9000, 10)
print(T90)
p75 = isa.pressure_isa(7500)
p90 = isa.pressure_isa(9000,10)
print(p75)
print(p90)
Explanation: Finally, suppose $h=7500m, \Delta T_{ISA}=0K$ and $h=9000m, \Delta T_{ISA}=+10K$
End of explanation
hp75 = isa.height_from_pressure_isa(p75)
hp90 = isa.height_from_pressure_isa(p90, 10)
print(hp75)
print(hp90)
Explanation: Now retrieve geopotential height, indicating the ISA deviation:
End of explanation
<END_TASK>