I'm not entirely sure if you're asking for a closed expression that gives the points, or just a description of these points and am going to assume the latter.
The first step in solving this problem is to change it into an easier one. We may rotate the sphere in such a way the point $P$ becomes the north pole $N=(0,0,1)$. If we call this rotation (which may be expressed as a $3\times3$ matrix) $R$, a point $Q$ is on the circle with radius $r$ around $P$ iff $R(Q)$ is on the circle with radius $r$ around $N$. This is due to the fact that a rotation of the sphere does not change distances on the sphere.
The next step is to consider these particular circles of radius $r$ around $N$. The best way to see this is to draw a sketch of the intersection of the sphere with the $yz$-plane, i.e. a circle of radius 1 around the origin. We can also see to points $P_1, P_2$ on the circle at distance $r$ from $N$. (I'll assume $r$ is small enough for this to make sense) These points should be such that $P_1$ is obtained from $P_2$ by reflecting around the $z$-axis. The spherical coordinate $\theta$ corresponding to these points can now easily be obtained: this angle is given as an angle wrt the $z$-axis, while $r$ is precisely the arc length corresponding to this angle on a circle with radius 1. This leads to the conclusion that $r=\theta$ in this special case.
The next thing to notice, is that all points at distance $r$ from $N$ have the same $z$-coordinate. By projecting either $P_1$ or $P_2$ at this $z$-coordinate (which again can be done insightful in the same picture), one finds this height to be $z_0=\cos(\theta)=\cos(r)$. Moreover, yet again from considering the same picture, the required circle is a circle parallel to the $xy$-plane with center $(0,0,z_0)$ and radius $\sin(\theta)=\sin(r)$.
We may thus conclude that the points at distance $r$ from $N$ are of the form $(\sin(r)\cos(\phi),\sin(r)\sin(\phi),cos(r))$ with $\phi\in[0,2\pi)$. Using the rotation matrix $R$ from before, we thus get the set of points satisfying your original question:
$$\{R^{-1}(\sin(r)\cos(\phi),\sin(r)\sin(\phi),cos(r))\mid \phi\in[0,2\pi)\}.$$
The tangent lines to the sphere form a cone with apex at the given point $\mathbf p$ and axis the line through $\mathbf p$ and the center $\mathbf c$ of the sphere. (To keep things clearer, I’ll mark inhomogeneous Cartesian coordinate vectors with a tilde; if there’s no tilde, this indicates a homogeneous coordinate vector.) The points of tangency lie on the intersection of the sphere with the polar plane of $\mathbf p$. This will be a circle in this polar plane, which is perpendicular to the cone’s axis, with center $\mathbf c'$ at the intersection of the plane and axis. To parameterize this circle, find the radius $r'$ via the Pythagorean theorem, and find a pair of unit vectors $\mathbf u$ and $\mathbf v$ that are orthogonal to the axis and each other. The circle is then $$\mathbf q = \mathbf c'+r'\mathbf u\cos\theta+r'\mathbf v\sin\theta.$$ The line through this point be found in any of the usual ways, such as the parametrization $(1-\lambda)\mathbf p+\lambda\mathbf q$.
To keep things clearer in the following, I’ll mark variables that represent inhomogeneous Cartesian coordinate vectors with a tilde; no tilde indicates a homogeneous coordinate vector. Representing the sphere by the matrix $$S = \left[\begin{array}{c|c} I_3 & -\tilde{\mathbf c} \\ \hline -\tilde{\mathbf c}^T & \|\tilde{\mathbf c}\|^2-r^2\end{array}\right],$$ the polar plane to $\mathbf p$ is $\mathbf\pi=S\mathbf p$. (The components of this vectors correspond to the coefficients of the implicit Cartesian equation of the plane.) Use your favorite method to find its intersection with the line through $\mathbf p$ and $\mathbf c$. Here, I’ll use the Plücker matrix of the line: $$\mathbf c' = (\mathbf p\mathbf c^T-\mathbf c\mathbf p^T)\mathbf\pi = (\mathbf c^T\mathbf\pi)\mathbf p - (\mathbf p^T\mathbf\pi)\mathbf c.$$ (The parenthesized quantities on the r.h.s. are just the dot products of $\mathbf p$ and $\mathbf c$ with $\mathbf\pi$.) The radius of the circle is found from $$r'^2 = r^2-\|\tilde{\mathbf c}-\tilde{\mathbf c}'\|^2.$$ For the unit vectors $\mathbf u$ and $\mathbf v$ of the parametrization, choose one of the cross products of $\tilde{\mathbf c}-\tilde{\mathbf p}$ with the unit basis vectors (at least two of them are guaranteed to be nonzero) and normalize it for $\mathbf u$, and then take $\mathbf u\times(\tilde{\mathbf c}-\tilde{\mathbf p})$, normalized, for $\mathbf v$. You can instead use $\tilde{\mathbf c}-\tilde{\mathbf c}'$ or $\tilde{\mathbf c}'-\tilde{\mathbf p}$ for this calculation; I’m not sure which will provide the best numerical stability.
For your example, I get $$S = \begin{bmatrix} 1&0&0&-100\\0&1&0&50\\0&0&1&-75\\-100&50&-75&17500\end{bmatrix}$$ from which $\mathbf\pi=[400:150:-75:-27500]$. Using the above formula for the axis-plane intersection gives $$\mathbf c' = [400:150:-75:-27500],$$ which corresponds to the inhomogeneous Cartesian coordinates $$\tilde{\mathbf c}' = \left({30500\over301},{14900\over301},{22500\over301}\right)$$ and the circle’s radius works out to be $r'=250\sqrt{\frac3{301}}$.
For the unit vector $\tilde{\mathbf u}$, I'll take the cross product that has the greatest norm, $(-150,400,0)^T$, so $\tilde{\mathbf u} = \frac1{\sqrt{73}}(-3,8,0)$ and (according to Mathematica; the radicals are starting to get a bit ugly) $\tilde{\mathbf v}=\left({24\over\sqrt{21973}},{9\over\sqrt{21973}},2\sqrt{{73\over301}}\right)^T$. The parametrization of the circle finally comes out to be approximately $$x = 101.329 - 8.76349 \cos\theta + 4.04095 \sin\theta \\ y= -49.5017 +
23.3693 \cos\theta + 1.51536 \sin\theta \\ z = 74.7508 +
24.5825 \sin\theta.$$
Depending on what it is you want to do with this, it might be more convenient to work with the tangent cone as a whole. It turns out that its matrix is easily computed from $\mathbf p$ and $C$: $$C_{cone} = (\mathbf p^TC\mathbf p)C - (C\mathbf p)(C\mathbf p)^T.$$ Note that the second term isn’t a dot product; it’s the tensor product of $C\mathbf p$ with itself.
Best Answer
I'm not sure to well understand your question, so I add a figure.
The figure is a plane section of your sphere passing from $P$. If I well understand the distance that you want is the length of the arc $PB$.( If this is wrong than my answer is wrong)
You know the radius of the cap, that is the arc $AB=\beta$. In this case the arc $PB$ is simply $\dfrac{\pi}{2}-\alpha-\beta$ , where $\alpha$ is the arc that fix the position of $P$ with respect to the equatorial plane of the sphere(its latitude).
If you know as radius of the cap the distance $CB$ than you can find $\beta=\arcsin (CB)$.