Describe what the following function does as clearly as you can, referring to equations from the equation sheet when appropriate
def function_H(Temp,height):
R_d=287.
g=9.8
num_temps=len(Temp)
inv_scale_height=0
for index in range(num_temps - 1):
delta_z=height[index+1] - height[index]
inv_scale_height=inv_scale_height + \
g/(R_d*Temp[index])*delta_z
avg_inv_scale_height=inv_scale_height/height[-1]
avg_scale_height=1/avg_inv_scale_height
return avg_scale_heightThe function is calculating the average scale height by doing this integral:
$$\frac{1}{H} = \int_0^\infty \frac{g}{R_d T} dz$$It takes as input two sequences (lists or numpy arrays) Temp and Height of the same length, and approximates the integral by summing the product $g/(R_d T) \Delta z$ for every layer $\Delta z$
Rewrite lines 4-8 using sum and diff to remove the loop (this will speed things up by about a factor of 100).
The following code uses np.sum and np.diff to remove the for loop::
def function_H_vec(Temp,height):
R_d=287.
g=9.8
num_temps=len(Temp)
function=g/(R_d*Temp[:-1])
delta_z=np.diff(height)
the_int=np.sum(function*delta_z)
the_avg=the_int/height[-1]
scale_height=1./the_avg
return scale_heightCheck this by calculating the exact hydrostatic pressure by integrating the hydrostatic equation:
$$\int_{p_0}^p d\log(p^\prime) = \int_0^z -\frac{g dz^\prime}{R_d T}$$$$log(p/p_0) = \int_0^z -\frac{g dz^\prime}{R_d T}$$Which in python, looks like this, where we're using np.cumsum to get intermediate pressure values at each height::
def press_int(heights,Temp,p0):
Rd=287.
g=9.8
dz=np.diff(heights)
function= -g/(Rd*Temp[:-1])
log_p_p0=np.cumsum(function*dz)
p_p0=np.exp(log_p_p0)
#
# omsert the surface as first value,
# where p/p0=1
#
p_p0=np.concatenate(([1.],p_p0))
press_vec=p0*p_p0
return press_vec
In [0]:
%matplotlib inline
In [0]:
from matplotlib import pyplot as plt
import numpy as np
def function_H(Temp,height):
R_d=287.
g=9.8
num_temps=len(Temp)
inv_scale_height=0
for index in range(num_temps - 1):
delta_z=height[index+1] - height[index]
inv_scale_height=inv_scale_height + \
g/(R_d*Temp[index])*delta_z
avg_inv_scale_height=inv_scale_height/height[-1]
avg_scale_height=1/avg_inv_scale_height
return avg_scale_height
def function_H_vec(Temp,height):
"""
find the average scale height
input: Temp (K), height (m): numpy vectors (1D)
output: scale_height (m)
"""
R_d=287.
g=9.8
num_temps=len(Temp)
function=g/(R_d*Temp[:-1])
delta_z=np.diff(height)
the_int=np.sum(function*delta_z)
the_avg=the_int/height[-1]
scale_height=1./the_avg
return scale_height
def press_int(heights,Temp,p0):
"""
integrate the hydrostatic equation
to get log pressure, then take exp
to get pressure
"""
Rd=287.
g=9.8
dz=np.diff(heights)
function= -g/(Rd*Temp[:-1])
log_p_p0=np.cumsum(function*dz)
p_p0=np.exp(log_p_p0)
p_p0=np.concatenate(([1.],p_p0))
press_vec=p0*p_p0
return press_vec
if __name__=="__main__":
ztop=2.e4
p0=1.e5
T0=280.
heights=np.linspace(0,ztop,100.)
Temp= T0 -7.e-3*heights
Temp=np.empty_like(heights)
Temp[:]=T0
press_vec=press_int(heights,Temp,p0)
scale_height=function_H(Temp,heights)
scale_height_II=function_H_vec(Temp,heights)
press_scale=p0*np.exp(-heights/scale_height)
plt.close('all')
fig=plt.figure(1,figsize=(9,9))
ax=fig.add_subplot(111)
ax.plot(press_vec*1.e-3,heights*1.e-3,label='exact integration')
ax.plot(press_scale*1.e-3,heights*1.e-3,'r+',label='exp using scale height')
ax.set_xlabel('pressure (kPa)')
ax.set_ylabel('height (km)')
ax.legend(loc='best')
plt.show()
In [0]: