Have you considered finding the intersections using an implicit form for the circles, $$\frac{x^2}{r^2} + \frac{y^2}{r^2} + ax + by + c = 0?$$ This representation doesn't have any coefficients that diverge as the circle approaches a straight line. To find intersections, you'll have to solve a quadratic equation whose leading coefficient could be zero or arbitrarily close to it, but the alternative form of the quadratic formula should be able to deal with that robustly.
You'll then have to do some jiggery-pokery to figure out whether the intersection points lie within the arcs. If the arc's bending angle is smaller than $\pi$, a projection onto the line joining the endpoints will suffice.
(Disclaimer: While all of this feels like it should work, I haven't analyzed it in any detail. Also, there could still be a problem when the circle is close to a line and you want the longer arc. But I can't imagine that's a case that would turn up in any practical application.)
Update: For a concrete example, here is the equation for a circular arc passing through the three points $(0,0)$, $(0.5, h)$, and $(1,0)$: $$\kappa^2 x^2 + \kappa^2 y^2 - \kappa^2 x - 2\eta y = 0,$$ where $$\begin{align}\kappa &= \frac{8h}{4h^2 + 1}, \\ \eta &= \frac{8h(4h^2-1)}{(4h^2+1)^2}.\end{align}$$ As you can see, the coefficients remain bounded as $h \to 0$.
Update 2: Wait, that equation becomes trivial if $h = 0$, which is bad. We really want something like $x^2/r + y^2/r + ax + by + c,$ i.e. multiply the previous expression through by $r$. Then for the same example, our equation becomes $$\kappa x^2 + \kappa y^2 - \kappa x - 2\eta' y = 0,$$ where $\eta' = (4h^2-1)/(4h^2+1)$. Here are some explicit values.
$h = 1/2$: $$2 x^2 + 2 y^2 - 2 x = 0,$$ $h = 0.01$: $$0.07997 x^2 + 0.07997 y^2 - 0.07997 x + 1.998 y = 0,$$ $h = 0$: $$2 y = 0.$$
By the way, in this format, the linear terms will always be simply $-2(x_0/r)x$ and $-2(y_0/r)y$, where the center of the circle is at $(x_0,y_0)$. As the center goes to infinity but the endpoints remain fixed, these coefficients remain bounded and nonzero (i.e. not both zero).
Your code is almost completely correct (and avoids the use of atan2, which is nice). The one problem you have is that when the vector from the disk-center to your point happens to point perpendicular to the disk, i.e., when the test point happens to be exactly above the disk's center, it falls apart: projecting that vector to the $xy$-plane gives the zero vector, whose tip does not lie on the disk's boundary circle.
Here's replacement code (mostly yours):
vec1 = vector between disk center and particle
vec1 = [(5 - 0), (5 - 0), (5 - 0)]
vec1 = [5,5,5]
pvec1 = vec1 with its z-coordinate set to 0.
if (norm (pvec1) = 0 ):
answer = sqrt( norm(vec1) * norm(vec1) + radius * radius )
return
unitvec1 = unit vector in direction of pvec1
unitvec1 = pvec1/norm(pvec1)
unitvec1 = [0.7071, 0.7071, 0]
vec2 = vector between disk center and point on the perimeter closest to the particle
vec2 = disk radius * unitvec1
vec2 = 3 * [0.7071, 0.7071, 0]
vec2 = [2.1..., 2.1..., 0]
vec3 = vector between particle and point on the perimeter closest to the particle
vec3 = vec1 - vec2
So the min distance is
norm(vec3) = 6.8087
That'll handle the special case on which your current code is failing. Note that you might want to replace the test "norm (pvec1) = 0" with "norm (pvec1) < .001" to make the code a bit more numerically stable, and introduce almost no error.
Best Answer
There are four cases to consider:
1) Endpoints of both arcs
2) An endpoint of one and an interior point of the other, which is on the line through that endpoint and the centre of the other arc.
3) Interior points of both arcs, which are on the line through the centres of the two arcs.
4) Intersections of the two arcs (thanks for pointing that out, Lopsy)