<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def solar_azimuth_analytical(latitude, hourangle, declination, zenith):
 """ Analytical expression of solar azimuth angle based on spherical trigonometry. Parameters latitude : numeric Latitude of location in radians. hourangle : numeric Hour angle in the local solar time in radians. declination : numeric Declination of the sun in radians. zenith : numeric Solar zenith angle in radians. Returns ------- azimuth : numeric Solar azimuth angle in radians. References [1] J. A. Duffie and W. A. Beckman, "Solar Engineering of Thermal Processes, 3rd Edition" pp. 14, J. Wiley and Sons, New York (2006) [2] J. H. Seinfeld and S. N. Pandis, "Atmospheric Chemistry and Physics" p. 132, J. Wiley (1998) [3] `Wikipedia: Solar Azimuth Angle <https://en.wikipedia.org/wiki/Solar_azimuth_angle>`_ [4] `PVCDROM: Azimuth Angle <http://www.pveducation.org/pvcdrom/2- properties-sunlight/azimuth-angle>`_ See Also -------- declination_spencer71 declination_cooper69 hour_angle solar_zenith_analytical """


    numer = (np.cos(zenith) * np.sin(latitude) - np.sin(declination))
    denom = (np.sin(zenith) * np.cos(latitude))

    # cases that would generate new NaN values are safely ignored here
    # since they are dealt with further below
    with np.errstate(invalid='ignore', divide='ignore'):
        cos_azi = numer / denom

    # when zero division occurs, use the limit value of the analytical
    # expression
    cos_azi = \
        np.where(np.isclose(denom,    0.0, rtol=0.0, atol=1e-8),  1.0, cos_azi)

    # when too many round-ups in floating point math take cos_azi beyond
    # 1.0, use 1.0
    cos_azi = \
        np.where(np.isclose(cos_azi,  1.0, rtol=0.0, atol=1e-8),  1.0, cos_azi)
    cos_azi = \
        np.where(np.isclose(cos_azi, -1.0, rtol=0.0, atol=1e-8), -1.0, cos_azi)

    # when NaN values occur in input, ignore and pass to output
    with np.errstate(invalid='ignore'):
        sign_ha = np.sign(hourangle)

    return sign_ha * np.arccos(cos_azi) + np.pi