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()

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()