Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

How to find all points along a great circle given two points on circle

  1. Jan 28, 2012 #1
    Hi,

    I've researched this problem all across the web and most answers involve finding the distance between two points on a great circle, mostly in nautical terms, using latitude and longitude, etc. but that isn't the answer I'm looking for. It's generally agreed if you have two points on a great circle that are not antipodal then there is only one great circle through those points. But nowhere can i find how to generate parametric formulas for a unique great circle or how to find other points on the circle *analytically* (without graphs.) I can specify two points that lie on the unit circle: (0,.5,-.866), (.5,0,-.866). These are separated by 90 degrees on the circle. I would like to be able to compute all the points on the circle (in rectangular coordinates) starting from 0 degrees (the first point) and in increments of tenths or hundreths of a degree to 360 degrees. Can anyone help me with this? My math level includes college trig and a little differential calculus.
     
  2. jcsd
  3. Jan 28, 2012 #2

    Simon Bridge

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member
    2016 Award

    A circle, being a conic section, usually needs three points to be specified.
    However, a great circle has the extra information that it is the biggest on a sphere.
    The radius of the sphere is normally given, as is the center, but I don't think both are needed.

    Presumably you have seen:
    http://mathworld.wolfram.com/GreatCircle.html
    Lat and Long are just spherical coordinates - presumably you can convert to cartesian?

    Also see hints:
    https://www.physicsforums.com/showthread.php?t=291195

    If the sphere is centered on the origin and you have two points on the surface A and B then the great circle through those two points is the intersection between the plane through A and B and the origin and the sphere. Presumably you know how to do that?
     
  4. Jan 28, 2012 #3
    I don't have any experience working with conic sections. My college trig only covered 2D not 3D objects. What I need are formulas in the form of x=f(xt),y=f(xt),z=f(xt), where xt is the angle on the circle and f is the formula with whatever functions are needed to map the circle's orbit. It's a unit sphere and the great circle has an origin of (0,0,0). I tried solving three equations with three variables using three known points on the sphere (another point on the circle plane is at 120 degrees: .866,-.5, 0), as suggested in the thread

    https://www.physicsforums.com/showthread.php?t=291195

    I came up with real values for the variables: x=sqrt(2), y=sqrt(6), etc. so that approach didn't work. The other approach, rotating from 2D to 3D, is hard for me to follow since I don't know what the angle of the plane is to the z-axis. I do know the circle plane is rotating at 45 degrees to the x/y plane, and the z plane is perpendicular to the x/y plane which lies on the plane of the monitor screen. Most examples of spherical and rectangular coordinates have the x and z planes reversed from this, so that makes conversions from spherical to rectangular coordinates and viceversa problematic, since I don't know if the parametric equations for the conversions are still valid for this setup. They don't seem to work as I would expect, but then my experience with spherical coordinates is limited.
     
  5. Jan 28, 2012 #4
    Are your two points always spaced by 90 degrees? If so, what I'd do is to imagine them as a pair of X-Y "axes", and try to construct a circle on this new "system of coordinates".

    Pick one of the vectors that looks to you as pointing to the "right" (define that to your convenience), and call that vector "X". In your example, I'd pick X=(0.5, 0, -0.866); the other one would point to the "front", and we call that one Y=(0, 0.5, -0.866). A cross-product X x Y would give you a third vector pointing "up" but, as it will turn out, you won't need it; just assume that there is one "Z" unit vector in your new system, and call it Z=(z1, z2, z3).

    If you write these three vectors as columns of a matrix,[tex]\left(\begin{array}{ccc}
    0.5 & 0 & z_1 \\
    0 & 0.5 & z_2 \\
    -0.866 & -0.866 & z_3
    \end{array}\right)[/tex]
    you have a matrix to convert from "your" X,Y,Z system into the ordinary X,Y,Z (1,0,0), (0,1,0), (0,0,1) (the system in which you expressed your vectors from the beginning).

    Now, over that slanted X-Y plane of "your" system, the parametric equation of a unit circle looks like x(t) = cos(t), y(t) = sin(t), z(t) = 0. In other words, the vector (cos(t), sin(t), 0) points somewhere in the unit circle you want, as the parameter "t" moves... except that the vector has this simple expression because it is in your "slanted" system of coordinates, which is rather useless to you. So, transform this vector by your matrix, in order to bring it into the ordinary, standard X,Y,Z system. You do that by multiplying the matrix and the vector:[tex]\left(\begin{array}{ccc}
    0.5 & 0 & z_1 \\
    0 & 0.5 & z_2 \\
    -0.866 & -0.866 & z_3
    \end{array}\right) \left(\begin{array}{c}
    \cos(t) \\
    \sin(t) \\
    0
    \end{array}\right) = \left(\begin{array}{c}
    0.5 \cos(t) \\
    0.5 \sin(t) \\
    -0.866 \cos(t) - 0.866 \sin(t)
    \end{array}\right)[/tex]
    and, if I didn't make any gross mistake, the components of that resulting vector in the right-hand side give you the parametric equations of your circle: the first component is the new "x(t)" function, the second the "y(t)", the third the "z(t)".

    Actually, at the start it doesn't really matter which vector you pick as "X" and which as "Y"... the circle will move in one direction or the other, that's all. As "t" increases, the point in the circle will move from the vector you designate as "X" towards the vector you name as "Y". The requisite, however, is that the points be 90-degrees apart on your circle.
     
  6. Jan 28, 2012 #5

    LCKurtz

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    Here's a suggestion I would try. You are given two points u and v on the unit sphere. Think of them as position vectors ##\vec u## and ##\vec v##. It is easy enough to calculate cross products. Calculate ##\vec w = (\vec u \times \vec v)\times \vec u##. Then ##\vec w## and ##\vec u## are unit vectors perpendicular to each other and in the plane of the circle. So a parameterization of the circle is$$
    \vec R(t) = \vec u\cos t + \vec w\sin t $$
     
  7. Jan 28, 2012 #6
    No, the two points aren't always spaced by 90 degrees. The points are easier to calculate visually when one component is zero, that's all.

    I have some new (revised) points, that appear to be accurate:
    0 degrees = (0x, .707y, -.707z)
    30 degrees = (.354x, .354y, -.866z)
    60 degrees = (.707x, 0y, -.707z)
    90 degrees = (.866x, -.354y, -.354z)

    When I plug 0 and 90 into your matrix, I don't get the same coordinates for 30 degrees, though. So I need a more general way of computing the points. Given that the circle plane is 45 degrees offset from the x/y plane, with the z plane perpendicular to the x/y plane and these four points should be enough to calculate all other points on the circle.

     
  8. Jan 29, 2012 #7

    Simon Bridge

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member
    2016 Award

    Conic sections are 2D objects. They are the circle, parabola, ellipse and hyperbola (also a line and a point but we don't normally count those).

    And you have no 3D computation lessons?
    Cross product? Matrix? Nothing? Anything?

    Where can I start?

    I'll start with vectors and the cross-product, LCKurtz suggestion, see if I can't put you on some sort of useful track to follow....

    Your two arbitrary points, P, Q, and the origin, O, all lie in the plane of the great circle, in 3D. I'm going to use the property that the cross product of two vectors is perpendicular to both of them.

    Put [itex]\vec{p}=\overrightarrow{OP}[/itex] and [itex]\vec{q}=\overrightarrow{OQ}[/itex] then the axis of the circle is [itex]\vec{v}=\vec{p}\times\vec{q}[/itex]

    If θ (=xt) is the angle in the plane of the great circle, from vector [itex]\vec{p}[/itex] [1] then the vector pointing at θ=90° must be [itex]\vec{u}=\vec{v}\times\vec{p}[/itex]

    The components of the vectors are the coordinates of points on the circle.
    eg. if p=(a,b,c)t then P=(a,b,c).

    Vectors p,u, and v play the roles of the axis unit vectors for the circle. Like the i,j, and k vectors you are used to.

    If the circle were in the x-y plane then a vector w pointing to an arbitrary point W on the circle would have a parametric equation [itex]\vec{w}=\vec{i}\cos\theta+\vec{j}\sin\theta[/itex] which gives you the usual x=cosθ and y=sinθ. But this is for the special case that the great circle is the equator.

    For an arbitrary great circle, a vector w pointing to point W on the circle, angle θ from p is like that but with p and u in place of i and j:
    [itex]\vec{w}=\vec{p}\cos\theta + \vec{u}\sin\theta[/itex]
    ... you should be able to take it from there.

    ---------------------
    [1] you'll probably have to change yours to whatever reference you are told to use.
     
  9. Jan 29, 2012 #8

    Simon Bridge

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member
    2016 Award

    That simplifies things immensely! Do you know the direction of the tilt too?!

    Let the circle plane be x',y' so the axis perpendicular to the circle-plane is z' - so z' is tilted to make an angle of 45 degrees to z (need direction).
    The equation of the circle in the plane of the circle is [itex]x^{\prime 2}+y^{\prime 2} = 1[/itex] ... you know how to parametrize this - all you need then is x' and y' in terms of x,y,z and θ.
     
  10. Jan 29, 2012 #9
    I think the tilt is -45 degrees (the circle plane is first tilted back 45 degrees then to the left toward the -x axis 45 degrees) Note: the z axis is perpendicular to the x/y plane which is the plane of the page, x axis is horizontal, y axis is vertical, not exactly a standard orientation which has the x or y axis perpendicular to the page. (See attachment.)

     

    Attached Files:

  11. Jan 29, 2012 #10

    Simon Bridge

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member
    2016 Award

    OK - cool. The x-z plane is "horizontal".

    So you really just want the equation of a circle radius 1 centered on the origin and tilted 45deg clockwise viewed down the +z axis.

    I don't see the parameter angle θ - is the angle in the plane of the circle or in the x-z plane?
    (Actually - I don't think it matters!)

    notice: y only depends on x; also the projection in the x-z plane is an ellipse - so x and z only depend on the angle ... which will be the parametric equation for the ellipse.
     
  12. Jan 29, 2012 #11
    Angle 0 is the angle in the circle, starting from the 12 o'clock position. The rest of angles 0-360 follow around the circle clockwise. The projection is actually a circle. It's only drawn as an ellipse to emphasize the angle it is tilted. The circle plane is tilted both 45 degrees from the x/y plane (the page) and 45 degrees from the z axis, so to get a true view of the circle plane you first have to rotate the circle plane (from the page) 45 degrees counter-clockwise toward the -z axis and then 45 degrees counter-clockwise toward the -x axis. I tried to do this using the standard rotation matrix. First I set xt=sin(t),yt=cos(t),zt=0. Then I rotated xt,yt,zt using rotation angle 0,-45, 0 and then again using -45, 0, 0. This gave me 0,.7071,-.7071 for x(0),y(0),z(0), which matched my first angle, but all other angles except 180 were different from what I expected. So what did I do wrong? I would really like someone to fully develop this into the correct parametric formulas for (x,y,z) if that is possible, and then I could study the derivation. The assumption has been that I understand vectors and matrixes well enough to finish the derivation, but I get confused by this stuff very easily...

     
  13. Jan 29, 2012 #12

    LCKurtz

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    @fracs49er: Have you tried my suggestion? For your initial two points you get $$
    \vec R(t) = \langle \frac 1 2 sin t,\frac 1 2 \cos t,\frac{\sqrt 3} 2 (\cos t - \sin t)\rangle$$for the circle. It's really pretty easy and it works whether your original vectors are perpendicular or not.
     
  14. Jan 29, 2012 #13
    Thanks for following this thread. There are problems with your solution: z is positive when angle = 0 (should be negative) and with angle 120 I get a z and radius greater than one. It's a unit circle with origin at 0,0,0! So maybe you made a mistake in your solution.
    Note: I posted a new set of points which are more accurate than the first two:
    0 degrees = (0x, .707y, -.707z)
    30 degrees = (.354x, .354y, -.866z)
    60 degrees = (.707x, 0y, -.707z)
    90 degrees = (.866x, -.354y, -.354z)

    Maybe you could redo your solution using these points and check that it works with them. Unfortunately I can't decipher how you are using the vectors and cross products to solve for x,y,z. But thanks for the effort!


     
  15. Jan 29, 2012 #14

    Simon Bridge

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member
    2016 Award

    But which is the "12O'Clock" position in terms of the axis?
    Nope - a tilted circle projects an ellipse on the plane it is tilted against. Draw the view looking down the y axis and you'll see.
    Trickier - the axis of rotation then is a line in the x-z plane, either x=z or x=-z ... this means that the axes of the x-z ellipse are rotated wrt the x and z axis. The projection onto the x-y plane will also be an ellipse.

    There are several ways to deal with this - the most direct way is to use a rotation matrix (or two of them) like you tried. If you have recently covered matrix rotation - then this would be the way you are expected to do it (I get the impression that this is off an exercise set in class that you'd like to explore more?) - the circle is rotated in the z axis followed by a rotation in the y axis. sin(45)=cos(45)=1/√2 ... easy matrix. eg. for consecutive +45deg rotations, 1st about k and then about j:
    [tex]\mathbf{R}=\mathbf{R}_y \mathbf{R}_z = \frac{1}{2}\left (
    \begin{array}{ccc}
    1 & 0 & 1\\ 0 & \sqrt{2} & 0\\ -1 & 0 & 1
    \end{array}
    \right ) \left (
    \begin{array}{ccc}
    1 & -1 & 0\\ 1 & 1 & 0\\0 & 0 & \sqrt{2}
    \end{array}
    \right )[/tex]... but you have to check that these are the correct rotations.
    I'm sure you would however you will benefit more in terms of your understanding if you do this yourself. I know it's painful and confusing but everyone has to go through this. There is enough information in the above and your course-notes for you to work out the general derivation for yourself.

    The general vector method for an arbitrary pair of points on the unit sphere has been shown to you several times - all that stuff about cross products.

    If you prefer the matrix example which needs the points to be 90 degrees to each other then you have been shown how to use the cross product to find two points 90 degrees from each other.

    You could also use the rotation matrix method (above) in general by working out the axis of rotation and the angle of rotation and putting those into the general rotation matrix.

    So you are spoiled for choice.
     
  16. Jan 29, 2012 #15

    Simon Bridge

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member
    2016 Award

    Thinking about it:
    you need to make the problem clear.

    you want to find the parametric equations for a unit circle that has been tilted with respect to the z-x plane.

    you need to specify the axis and angle of the tilt - specify the axis as a vector and give the tilt direction by the right-hand-screw rule wrt this vector.

    if we imagine the axial vector to the circle starts out pointing in the +y direction, then the angle it ends up forming wrt the y axis is the tilt angle.

    From your description I am imagining a tilt angle of +45 degrees about (1,0,1).
    [see attachment - green=x-axis, red=z-axis, magenta=y-axis]

    that should clear up your thinking somewhat.
     

    Attached Files:

    Last edited: Jan 29, 2012
  17. Jan 29, 2012 #16
    Did you see the attachment I posted? It shows the circle plane clearly, and it's not just tilted with respect to the z-axis. Actually it is rotating at a 45 degree angle to the x/y axis too, in addition to being skewed with respect to the z-axis.

    Your assumption is I'm a math student, but I'm not, though I like the symmetry of mathematics in action, even if I don't understand how a lot of higher-level math works. My college days are long gone, and I really don't want to wrack my brains anymore dredging up a math solution that I'll probably use only once. I didn't think the Physics Forum was limited to students only, but maybe that's who mostly are members and who it's only suppose to help, so perhaps I'm barking up the wrong tree.

    Actually I'm a programmer trying to help a self-made quantum physicist illustrate his life-long thesis, and some of his mathematics has got me going bonkers trying to program it. Even he can't remember (twenty years back) how he derived the initial angles for his circle plane, which he calls a "frequency vector."

    So if you don't have a pain-free solution to the problem of finding points along the circluar vector then we'll have to do without it for now. Then I won't need to waste any more of yours or my time... I thought the solution to this would be trivial to guys of your math ability, and wouldn't need a long drawn-out consideration. Thanks for your time, anyway.

     
  18. Jan 30, 2012 #17

    Simon Bridge

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member
    2016 Award

    I saw that - did you see mine?

    Yours looks like -45deg about (0,0,1) followed by -45deg about (0,1,0)
    Can you confirm?
    [x and z both > 0 means y <0?]

    It got drawn out because you started out with one description of the problem and then changed it... twice now.

    It was trivial - and the complete general solutions are above.
    You should be able to get a computer to crunch the numbers for you - just use a math library which can handle matrixes and vectors.

    I assumed you were interested in learning. You don't have to be enrolled in an education program to be a student. But you are correct, if you don't want to learn, then you are barking up the wrong tree. PF is not a free math consultancy - the service you are asking for is expensive.

    Oh dear - the "frequency vector" is a common name for the angular velocity vector and it is the normal vector to the circle not the circle itself. The circle should have an arrow on it representing rotation. In QM we'd usually use the angular momentum instead.

    ... then perhaps you will be more used to rotations as quaternions? Can you translate vector notation into source code? Do you know how to code for rotations? For cross product?


    This is what I have so far:
    [tex]\left (
    \begin{array}{c}
    x\\y\\z
    \end{array} \right ) = \frac{1}{2}

    \left (
    \begin{array}{ccc}
    0 & \sqrt{2} & 0\\ 1 & 0 & 1\\ -1 & 0 & 1
    \end{array} \right )
    \left (
    \begin{array}{ccc}
    1 & 1 & 0\\ -1 & 1 & 0\\ 0 & 0 & \sqrt{2}
    \end{array} \right )

    \left (
    \begin{array}{c}
    \cos t\\0\\ \sin t
    \end{array} \right )

    [/tex]... but for t=0, x,y,z is not the same as yours - this is because my parameterization starts at a different place.

    You can also do this by parameterizing the ellipses in the x-z and x-y projections and solving the simultaneous equations.

    I do actually get a different result from LCKurtz from his method....
    you should be able to get your computer to do his cross products quite happily.
    The parameterization should even start in the right place if you pick u for t=0.
     
  19. Jan 30, 2012 #18
    Life is a continuous process of learning. So you could say I am a student of life, but I don't always have time to do my homework. I never took calculus in college so vectors and cross-products are strange sounding terms to me. Matrices I have a better understanding of, but if I don't work with them often, I might need a refresher course. I have books on trig and calculus that sometime come in handy to review and understand things like matrices, derivatives and integals, etc. There's also a section on cross products which I could probably program into a 3D function. Velocity vectors and quantum physics may remain beyond my grasp but we all can't be rocket scientists. One question I have about LCKurtz's cross-product method: R⃗(t)=u⃗cost+w⃗sint. To get w you need to multiply (u*w)*u (two consecative cross-products). Then you multiply the u vector against cost and the w vector against sint, both of which are scalar operations, I think they would be termed, such that R(t)= (1,2,3)cost = (1cost,2cost,3cost) + (3,4,5)sint = (3sint,4sint,5sint) = (1cost+3sint,2cost+4sint,3cost+5sint). Did I get that right? If so, I hope that solves my problem. It is certainly the most straight-forward looking solution.

     
  20. Jan 30, 2012 #19

    LCKurtz

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    Hello fracs49er. I just got back to this thread since we have company. Your steps are a bit sketchy but your final answer is correct as to how you simplify things. I did have a mistake in the formula I posted. You need to remember that ##\vec w## must be a unit vector. If the unit vectors ##\vec u## and ##\vec v## are perpendicular to each other then ##\vec w## will automatically be a unit vector. Otherwise it needs to be divided by its length. For your original two points the correct parameterization is $$
    \vec R(t) = \langle -\frac 2 {\sqrt 7}\sin(t),\frac 1 2 \cos t +\frac{3\sqrt 7}{14}\sin t,
    -\frac{\sqrt 3} 2 \cos t + \frac{\sqrt{21}}{14}\sin t\rangle$$I hope I don't have any typos. Here's a picture:
    greatcircle.jpg
     
    Last edited: Jan 30, 2012
  21. Jan 30, 2012 #20
    Your picture certainly looks like how the velocity vector appears on my diagram (Good job!), but then those two original points were supposedly in error, and the points I later posted supposedly correct. I tried setting up the cross products as in your formula, using the updated points for 0 and 90 degrees (0, .707, -.707 and .866,-.354,-.354.) That returns the same identical points when I enter 0 and 90 for t, and all angles for t that are in 90 increments from 90 also appear with the expected values. But when I enter 30 or any other angle that isn't a multiple of 90, then I get points that don't agree with what my associate had calculated, but they are in the approximate range that would be expected. For example, for 30 degrees, instead of .354,.354, -.866, your formula returns .433,.436,-.789. So I don't know what that means. Either my associate was wrong in his original calculations, or your formula only works with angles of 90 degree multiples, or my code is buggy in some way. (Note: using 30 degrees for vector v doesn't work for any angle except 0.) Did you plot your picture using continuous angles 0-360 outputted by your formula or did you simulate the curve in Mathematica? If the former, I have to think that my associate's earlier data for the frequency vector points was flawed in some way. Given that the two input points get handled by your formula correctly, I don't think I made any mistake in my coding, but you might check that:

    void fvector(hcomplex *u, hcomplex *v, hcomplex *ft, double t){
    hcomplex a,b;

    a.x=u->y*v->z-u->z*v->y; // u*v
    a.y=u->z*v->x-u->x*v->z;
    a.z=u->x*v->y-u->y*v->x;
    b.x=(a.y*u->z-a.z*u->y)*sin(t); // wsint
    b.y=(a.z*u->x-a.x*u->z)*sin(t);
    b.z=(a.x*u->y-a.y*u->x)*sin(t);
    ft->x=h1->x*cos(t)+b.x; // ucost + wsint
    ft->y=h1->y*cos(t)+b.y;
    ft->z=h1->z*cos(t)+b.z;
    }

    u.x=0;
    u.y=sqrt(.5); // .707
    u.z=-sqrt(.5); //-.707
    v.x=sqrt(.75); // .866
    v.y=-sqrt(.125); // -.354
    v.z=-sqrt(.125); // -.354
    t=angle*PI/180.0; // convert to radians
    fvector(&u,&v,&ft,t);
    sprintf(s,"%0.3f,%0.3f,%0.3f",ft.x,ft.y,ft.z);
    AfxMessageBox(s);


     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: How to find all points along a great circle given two points on circle
Loading...