Adding to @whuber's answer here is python code to reproduce the ideas he explained, as well to validate the results.

    from scipy.stats import beta
    import numpy as np
    
    def rsphere(n,n_dim,rad=1.0):
        X = np.random.normal(size=(n,n_dim+1))
        X_norm=np.divide(np.linalg.norm(X,axis=1),rad)
        X = np.divide(X.T,X_norm).T
        return X
    
    def reflect(X, target):
        n_dim = X.shape[1]
        pole = np.array([0]*(n_dim-1) + [1])
        v = pole - target
        X_new = X - np.outer(2.0/np.dot(v,v) * np.matmul(X,v).T, v)
        return X_new
    
    
    r_max=0.5 # max radius to sample
    n_dim=200 # number of dimension of the space
    n_sphere=n_dim-1 # dimension of the n-sphere
    n = int(1e3) # numer of points to sample
    
    target = np.zeros(n_dim) # target
    target[-2]=1 
    pole = np.array([0]*(n_dim-1) + [1]) # north pole

    # generate points
    X = rsphere(n,n_sphere,rad=1.)
    rads = np.random.uniform(0.0,r_max,n)
    q = beta.cdf(np.square(rads)/4.0, n_sphere/2.0, n_sphere/2.0)
    unif = np.random.uniform(0,q,n)
    z = 1 - 2*beta.ppf(unif, n_sphere/2.0, n_sphere/2.0)
    rho = np.sqrt(1.0 - np.square(z))
    X_pole = np.hstack((rsphere(n, n_sphere-1, rad=rho),z.reshape(-1,1)))
    X_nei = reflect(X_pole, target)

To validate results we can look at a histogram of distances of the generated points vs the samples points:

    import matplotlib.pyplot as plt
    X_dist = np.linalg.norm(X_nei-target,axis=1)
    plt.hist(X_dist)
    plt.xlabel('$\|z_i-z_{target} \|$')
    plt.show()

[![Distribution of distances][1]][1]
And to visualize the 3D case:

    from mpl_toolkits.mplot3d import Axes3D
    import matplotlib.pyplot as plt
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.set_aspect('equal')
    # Create a sphere
    r = 1
    pi = np.pi
    cos = np.cos
    sin = np.sin
    phi, theta = np.mgrid[0.0:pi:40j, 0.0:2.0*pi:40j]
    x = r*sin(phi)*cos(theta)
    y = r*sin(phi)*sin(theta)
    z = r*cos(phi)
    # sphere
    ax.plot_surface(
        x, y, z,  rstride=1, cstride=1, alpha=0.25, linewidth=0)
    # pole
    ax.scatter(X_pole[:,0],X_pole[:,1],X_pole[:,2],c='b',alpha=0.5)
    ax.scatter(pole[0],pole[1],pole[2],s=40,c='k')
    # target
    ax.scatter(X_nei[:,0],X_nei[:,1],X_nei[:,2],c='r',alpha=0.5)
    ax.scatter(target[0],target[1],target[2],s=40,c='k')
    plt.show()

[![Sphere sampling][2]][2]


  [1]: https://i.sstatic.net/trSIm.png
  [2]: https://i.sstatic.net/bFhua.png