Solve this particular system of equations

Okay so I have 2 equations, for the latitude and longitude, of a sphere when looking at it from a distance. The latitude and longitude can both be found by giving the equation a value for theta (the rotation angle from the vertical) and phi (the azimuthal angle) - just regular spherical polar coordinates.

i.e. theta = 70, phi = 90 ==> b (latitude) = 34, l (longitude) = 0. For example.

My question is now: For any given value of l and b, I would like to find the path length through the sphere. So any given (l,b) pair will give me two solutions of both theta and phi, which when the radius of the sphere is known (which I do), will give me a way of measuring the distance from one intersection point of the surface to the other.

How exactly do I go about doing this? I tried using Matlab's symbolic toolbox solve function, but it seems to me that by doing that, it will be extremely tedious, as I will want this to be as detailed as possible, so I will need to know (theta,phi) pairs for many (l,b) pairs. So ideally I'd think I would want to write a function yeah? but by using the symbolic toolbox solve, the equation entered has to be a string (i.e. S = solve('34 = f(theta,phi), '10 = g(theta,phi), theta, phi) where here, I'm trying to solve for l = 10, b = 34). Ideally, I would want to be able to enter a vector of b values and l values, and get a vector solution out for both theta and phi - basically do it in one hit if I can.

Please let me know if this is possible, and if you need any further clarification or any of my formulae then just let me know!

Thanks,

Joshua

Answer 1

[John D'Errico]

A confusing question. Sigh. As Mohammad said, a bit unclear. I'd call that a bit of an understatement.

To me, my interpretation is that you have a point on the surface of a sphere, given in terms of latitude and longitude. Then given another point in space that is in general external to the sphere, you wish to find the two intersections of the line that passes through these two points and the sphere. You say that you want the path length through the sphere, so given those two intersection points, the Euclidean distance between them is trivial to compute.

I would suppose that if the point on the sphere is on the back side of the sphere, from the point you are viewing from, then you want to know the distance through the sphere along that line.

Is this your question? If so, then be more clear, because there are three people here who each have chosen totally different interpretations of your question. I'll bet if we looked, we could find another interpretation too.

So what is your question, really? BE CLEAR! While I think I have interpreted it correctly, I'm a bit unsure, since two others have gotten completely different conclusions about what you need. I know, they are wrong. Or maybe its just me.

If your question is what I have said it is, then the answer is simple, doable with pencil and paper, though I'd probably use symbolic tools because, well, I'm lazy.

The point on the surface of a sphere can be converted to cartesian coordinates. Assume it is centered at the origin. If not, then it is trivial to translate the entire problem. After all, in the end, what you want is a distance. So a simple translation has no affect on distance. In fact, translating the origin of the sphere to (0,0,0) only makes the problem both simpler and better posed numerically. Then write the equation of the sphere as a basic quadratic form.

x^2 + y^2 + z^2 = r^2

The line through two points is simply written in a parametric form like this:

P(t) = P0*(1-t) + P1*t
     = P0 + t*(P1 - P0)

where P0 and P1 are vectors, defined by the two points in question. So any point on the line is given for some value of t. When t = 0, we get the point P0. When t = 1, we get P1 out. A nice thing about this linear parametric form is that it is easy to know where a point lies along that line, simply by knowing the value of t. Sometimes this parametric form is known as a convex combination. So if t is between 0 and 1, then the point located is between P0 and P1. Negative values of t, or values of t greater than 1 also have easily interpreted meanings.

Now, while I could go through the algebra, I'd only bother to do so IF it were of interest. It is not that hard to be honest. (On the other hand, I've too often invested a serious amount of effort, only to find the question was not as I interpreted it. I won't do so here. At least, I won't do so until I know what the real question is, and if I have interpreted what you said correctly.)

[Joshua D'Agostino]

Of course you've interpreted it correctly, I'm sorry but I fail to see another way of interpreting it... but in any case, do you mind elaborating please? Would really help! Thanks

[John D'Errico]

As I said, without loss of generality, the sphere can be centered at (0,0,0) by translating the entire problem. So the sphere has equation

x^2 + y^2 + z^2 = r^2

And we have a line with equation

P(t) = P0 + t*(P1 - P0)

Here I'll call P0 the given point on the sphere, and P1 is the given point external to the sphere.

P0 = [x0,y0,z0]
P1 = [x1,y1,z1]

So as a function of the line parameter t, any point on the line can be written as

P_t = P0 + t*(P1 - P0)
P_t = [x0 - t*(x0 - x1), y0 - t*(y0 - y1), z0 - t*(z0 - z1)]

Substitute the values for x,y,z into the sphere equation.

subs(x^2+y^2+z^2-r^2,{x,y,z},{x0 - t*(x0 - x1), y0 - t*(y0 - y1), z0 - t*(z0 - z1)})

(x0 - t*(x0 - x1))^2 + (y0 - t*(y0 - y1))^2 + (z0 - t*(z0 - z1))^2 - r^2

Again, you should see that when t == 0, this just reduces to the sphere equation. Solve for the two values of t.

solve(ans,t)
ans =
((r^2*x0^2 - 2*r^2*x0*x1 + r^2*x1^2 + r^2*y0^2 - 2*r^2*y0*y1 + r^2*y1^2 + r^2*z0^2 - 2*r^2*z0*z1 + r^2*z1^2 - x0^2*y1^2 - x0^2*z1^2 + 2*x0*x1*y0*y1 + 2*x0*x1*z0*z1 - x1^2*y0^2 - x1^2*z0^2 - y0^2*z1^2 + 2*y0*y1*z0*z1 - y1^2*z0^2)^(½) - y0*y1 - z0*z1 - x0*x1 + x0^2 + y0^2 + z0^2)/(x0^2 - 2*x0*x1 + x1^2 + y0^2 - 2*y0*y1 + y1^2 + z0^2 - 2*z0*z1 + z1^2)

 -(x0*x1 + y0*y1 + z0*z1 + (r^2*x0^2 - 2*r^2*x0*x1 + r^2*x1^2 + r^2*y0^2 - 2*r^2*y0*y1 + r^2*y1^2 + r^2*z0^2 - 2*r^2*z0*z1 + r^2*z1^2 - x0^2*y1^2 - x0^2*z1^2 + 2*x0*x1*y0*y1 + 2*x0*x1*z0*z1 - x1^2*y0^2 - x1^2*z0^2 - y0^2*z1^2 + 2*y0*y1*z0*z1 - y1^2*z0^2)^(1/2) - x0^2 - y0^2 - z0^2)/(x0^2 - 2*x0*x1 + x1^2 + y0^2 - 2*y0*y1 + y1^2 + z0^2 - 2*z0*z1 + z1^2)

Yes, it is a bit messy. I have a funny feeling that pencil and paper would yield a less messy solution. But who cares? I'm still feeling lazy.

Lets try it out. I'll pick a sphere of radius 3, and a point on the sphere at [2 2 -1]. Take the external point as, say [1 1 5]. I've chosen these points so the line must pass through the sphere.

I'll first write a pair of anonymous functions to do the computations for us.

T = @(r,x0,y0,z0,x1,y1,z1) [(sqrt(r^2*x0^2 - 2*r^2*x0*x1 + r^2*x1^2 + r^2*y0^2 - 2*r^2*y0*y1 + r^2*y1^2 + r^2*z0^2 - 2*r^2*z0*z1 + r^2*z1^2 - x0^2*y1^2 - x0^2*z1^2 + 2*x0*x1*y0*y1 + 2*x0*x1*z0*z1 - x1^2*y0^2 - x1^2*z0^2 - y0^2*z1^2 + 2*y0*y1*z0*z1 - y1^2*z0^2) - y0*y1 - z0*z1 - x0*x1 + x0^2 + y0^2 + z0^2)/(x0^2 - 2*x0*x1 + x1^2 + y0^2 - 2*y0*y1 + y1^2 + z0^2 - 2*z0*z1 + z1^2);
 -(x0*x1 + y0*y1 + z0*z1 + sqrt(r^2*x0^2 - 2*r^2*x0*x1 + r^2*x1^2 + r^2*y0^2 - 2*r^2*y0*y1 + r^2*y1^2 + r^2*z0^2 - 2*r^2*z0*z1 + r^2*z1^2 - x0^2*y1^2 - x0^2*z1^2 + 2*x0*x1*y0*y1 + 2*x0*x1*z0*z1 - x1^2*y0^2 - x1^2*z0^2 - y0^2*z1^2 + 2*y0*y1*z0*z1 - y1^2*z0^2) - x0^2 - y0^2 - z0^2)/(x0^2 - 2*x0*x1 + x1^2 + y0^2 - 2*y0*y1 + y1^2 + z0^2 - 2*z0*z1 + z1^2)];

P_t = @(t,x0,y0,z0,x1,y1,z1) [x0 - t*(x0 - x1), y0 - t*(y0 - y1), z0 - t*(z0 - z1)];

Tinter = T(3, 2, 2, -1, 1, 1, 5)
Tinter =
       0.526315789473684
                       0

As I said before, when t=0, we get back P0, so that works. The second point of intersection is at t=10/19=0.526315789473684...

We can actually compute the intersection point as...

P_t(10/19, 2, 2, -1, 1, 1, 5)
ans =
    1.47368421052632  1.47368421052632  2.15789473684211

and verify that it lies on the sphere.

norm(ans)
ans =
    3

That does not really matter, since it was only the distance between the points you care about. We can get that as

norm(diff(P_t(Tinter, 2, 2, -1, 1, 1, 5)))
ans =
        3.24442842261525

So you can do the entire computation in essentially two lines for any pair of points.

[Joshua D'Agostino]

Oh wow, thank you!

Answer 2

[Chad Greene]

For path length through a sphere, transform your spherical coordinates to 3D cartesian coordinates x,y,z with sph2cart. Then the distance from a to b is

sqrt((xb-xa)^2 + (yb-ya)^2 + (zb-za)^2)