def blend_air_and_wave_velocity_cpp(wave_cpp, air_cpp, elevation_cpp,
    direction, eta_eps, air=None, psi_wave_cpp=None, psi_air_cpp=None,
    slope_cpp=None):
    """
    Return C++ code which blends the velocities in the water into the velocities
    in the air in such a way that the C++ code replicates the Python results
    from the blend_air_and_wave_velocities() function
    """
    if air is None:
        return 'x[2] < (%s) + %r ? (%s) : (%s)' % (elevation_cpp, eta_eps,
            wave_cpp, '0.0')
    cpp = (
        """[&]() {{
        const double elev = ({ecpp});
        
        const double val_water = ({uw_cpp});
        if (x[2] < elev + {eps!r}) {{
            // The point is below the free surface
            return val_water;
        }}
        
        const double d = {d!r};
        const double val_air = ({ua_cpp});
        if (x[2] < d) {{
            // The point is in the blending zone
            const double Z = (x[2] - elev) / (d - elev);
            
            // Cubic smoothstep
            //const double f = Z * Z * (3 - 2 * Z);
            //const double dfdZ = Z * (6 - 6 * Z);
            
            // Fift order smootherstep
            const double f = Z * Z * Z * (Z * (Z * 6 - 15) + 10);
            const double dfdZ = Z * Z * (30 + Z * (30 * Z - 60));
            
            // Derivatives needed in the product rules
            const double dZdx = (x[2] - d) / pow(d - elev, 2) * ({slope_cpp});
            const double dZdz = 1.0 / (d - elev);
            const double dfdx = dfdZ * dZdx;
            const double dfdz = dfdZ * dZdz;
            
            return (1 - f) * val_water + f * val_air + %r * (%s);
        }}
        
        // The point is above the blending zone
        return val_air;
    }}()"""
        .format(ecpp=elevation_cpp, uw_cpp=wave_cpp, ua_cpp=air_cpp, eps=
        eta_eps, d=air.blending_height + air.depth_water, slope_cpp=slope_cpp))
    if direction == 'x':
        sign = -1
        prcode = 'dfdz * (%s) - dfdz * (%s)' % (psi_wave_cpp, psi_air_cpp)
    elif direction == 'z':
        sign = +1
        prcode = 'dfdx * (%s) - dfdx * (%s)' % (psi_wave_cpp, psi_air_cpp)
    return cpp % (sign, prcode)
