```Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: how to judge if an arc of a great circle crosses a quadrangle and calculate the positions of the intersects of  on a sphere surface
Date: Fri, 12 Mar 2010 22:34:06 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 35
Message-ID: <hnefgu\$rq1\$1@fred.mathworks.com>
References: <hnd5rr\$4eq\$1@fred.mathworks.com>
NNTP-Posting-Host: webapp-02-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1268433246 28481 172.30.248.37 (12 Mar 2010 22:34:06 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Fri, 12 Mar 2010 22:34:06 +0000 (UTC)
Xref: news.mathworks.com comp.soft-sys.matlab:616393

"Ha " <scifiles@126.com> wrote in message <hnd5rr\$4eq\$1@fred.mathworks.com>...
> The story below happens on Earth surface.
>
> There is an arc of a great circle (AB), the two endpoints of which are placed at (latA,lonA) and (latB, lonB). Beside, there is a quadrangle. The latitude and longitude of the four vertices (V1,V2,V3,V4) of the quadrangle are also known. The edges (e.g. V1V2,V2V4,...) of the quadrangle all lie on great cirles. Then, how can I know if the arc AB  crosses the quadrangle and calculate their intersects?
-----------------
Your quadrangle question can be broken down into finding whether a pair of great circle arcs, each defined by a pair of points, will intersect, and if so what that point of intersection is.

I would suggest transforming your longitude and latitude coordinates (in radians) over to cartesian coordinates using the equations:

x = cos(theta)*cos(phi);
y = sin(theta)*cos(phi);
z = sin(phi);

where theta is longitude measured positive going east of greenwich and negative west, and where phi is latitude measured positive going north from the equator and negative south.  We should have -pi<=theta<=pi and -pi/2<=phi<=pi/2.  These cartesian coordinates correspond to a hypothetical spherical "earth" of unit radius, but that does not interfere in the following computations.

The following is a sketch of how you might then proceed.  Let a, b, c, and d be vectors of the cartesian coordinate endpoints for the two arcs a-b and c-d obtained in this way.  Carry out the following computations:

p = cross(a,b); % p is along the normal to the plane of arc a-to-b
q = cross(c,d); % Similarly for q and arc c-to-d
t = cross(p,q); % t is along the line of intersection of the planes
s1 = dot(cross(p,a),t);
s2 = dot(cross(b,p),t);
s3 = dot(cross(q,c),t);
s4 = dot(cross(d,q),t);

You will find that if s1, s2, s3, and s4 are all of the same sign, the arcs will then and only then intersect.  In that case they intersect along +t if they are all positive and along -t if all are negative.

If they do intersect, you can transform the corresponding vector, t or -t, back into longitude and latitude (without worrying about its length.)  Letting x,y,z be its cartesian coordinates this reverse transformation can be accomplished this way:

theta = atan2(y,x);
phi = atan2(z,sqrt(x^2+y^2));

The above procedure ignores the indeterminate cases when 1) two ends of an arc are coincident or at antipodes (in which case p or q is of zero length) or 2) when the arcs lie in the same plane (in which case t is of zero length).  I will let you deal with these.

Roger Stafford
```