def METHOD_NAME( z, r ):
    dens = np.ones_like( z )
    # Make the density smooth at rmax
    dens = np.where( r > rmax-smooth_r,
        np.cos(0.5*np.pi*(r-smooth_r)/smooth_r)**2, dens)
    # Make the density 0 below p_zmin
    dens = np.where( z < p_zmin, 0., dens )
    # Make a linear ramp
    dens = np.where( (z>=p_zmin) & (z<p_zmin+ramp),
                     (z-p_zmin)/ramp*dens, dens )
    return( dens )