Plotting a Circle with lat/long

[Eccentric]

Scenario

I was doing the UCGO Europa Base mission, and decided to get creative with placing lights in a circle around the Arrow Freighter as a kind of temporary landing pad.

I thought I was being very clever when I jumped into a spreadsheet and whipped up some calculations with rudimentary trigonometry to figure out the lat/long points to place the lights at with a radius of whatever 0.006 degrees worked out to be.

Problem

You've probably already realized what I ran into. The base is at about 68 N and it appears the longitudinal distance equal to 0.006 degree is about 1/3 or 1/4 by eyeball of what it equals for the latitudinal distance. I got a nice ellipse with a North-South orientation. Suited the Arrow freighter since it was aligned that way, but was not what I wanted.

Attempted Solution

I figured it would be straight forward to go find one of those Great Circle distance plotter, and plug that into my spreadsheet. Get a distance for say 10 deg of lat, another for 10 deg of long, and the ratio should give me what I need to adjust my numbers, right?

Bad Math?

I found many places where formula where offered. One such, the so-called Haversine formula seemed straight forward. I found it here. Happily, they have a nice spot where I could put in my numbers and see how my spreadsheet said the same thing.

Except I don't get the same distance. Not even close. And not wrong by a factor of 10, or 0.1. No, just way too big.

I checked, checked again. Started from scratch and went to great lengths to be sure I was using the right numbers in the right places. All in radians, yes yes yes.

I have a headache now, and it still doesn't work.

Help? Maybe?

 

[Fizyk]

Here is how I would calculate this:

First, instead of using lat/long for denoting positions on the surface of a sphere, it may be (and will be in this case) better to use vectors from the center to given points.

For a lat/long of [math](\theta, \phi)[/math] the vector will be:

[math]\vec{v} = (\cos \theta \cos \phi, \cos \theta \sin \phi, \sin \theta)[/math]

Vectors for points on the circle of radius [math]\alpha[/math] degrees will be vectors that are at angle [math]\alpha[/math] to vector [math]\vec{v}[/math]. To generate such vectors, we will need two unit vectors [math]\vec{a}[/math] and [math]\vec{b}[/math] perpendicular to [math]\vec{v}[/math].

First such vector can be generated from [math]\vec{z} = (0,0,1)[/math], unless [math]\vec{v}[/math] is parallel to it (which it isn't if your base is at 68 deg N). Generation looks like this:

[math]\vec{a_0} = \vec{z} - (\vec{z} \cdot \vec{v}) \vec{v}[/math]
[math]\vec{a} = \frac{\vec{a_0}}{|\vec{a_0}|}[/math]

Explicitly:
[math]\vec{a} = (-\sin \theta \cos \phi, -\sin\theta\sin\phi, \cos\theta[/math]

[math]\vec{b}[/math] can be calculated for example as [math]\vec{a} \times \vec{v}[/math], which gives:

[math]\vec{b} = (-\sin\phi, \cos\phi, 0)[/math]

Now, to get a point displaced by [math]\alpha[/math] in the direction [math]\psi[/math] (it will correspond to heading), you need to calculate:
[math]\vec{u} = \vec{v} + \tan\alpha\cos\psi \vec{a} + \tan\alpha\sin\psi\vec{b}[/math]

This vector won't be a unit vector, but the angle between this vector and [math]\vec{v}[/math] will be exactly [math]\alpha[/math], so this is what we are looking for.

This will equal:
[math]u_x = \cos\theta\cos\phi - \tan\alpha\cos\psi\sin\theta\cos\phi - \tan\alpha\sin\psi\sin\phi[/math]
[math]u_y = \cos\theta\sin\phi - \tan\alpha\cos\psi\sin\theta\sin\phi + \tan\alpha\sin\psi\cos\phi[/math]
[math]u_z = \sin\theta + \tan\alpha\cos\psi\cos\theta[/math]
(It wouldn't fit in one equation, so I wrote the coordinates separately )

To calculate lat/long from this, you can use:
[math]\theta' = \arcsin \frac{u_z}{|\vec{u}|} = \arcsin \frac{\sin\theta + \tan\alpha\cos\psi\cos\theta}{\sqrt{1+\tan^2\alpha}} = \arcsin(\cos\alpha\sin\theta + \sin\alpha\cos\psi\cos\theta)[/math]

[math]\phi' = \arctan\frac{u_y}{u_x} = \arctan\frac{\cos\theta\sin\phi - \tan\alpha\cos\psi\sin\theta\sin\phi + \tan\alpha\sin\psi\cos\phi}{\cos\theta\cos\phi - \tan\alpha\cos\psi\sin\theta\cos\phi - \tan\alpha\sin\psi\sin\phi}[/math]

Your [math]\alpha[/math] here will be 0.006 degrees, change [math]\psi[/math] to get different points on the circle.

 

[Eccentric]

Here is how I would calculate this:

First, instead of using lat/long for denoting positions on the surface of a sphere, it may be (and will be in this case) better to use vectors from the center to given points.

... snip

Before I sequester myself in a dusty room with food and water for about a week while I come to understand what you have given me, let me do this first...

More Information

While I think I can come to use your formulas (I like math, but have never gotten into it that deeply) I should further explain that all I am doing is attempting to place objects (UMMU Lights in this case) in a circle on the surface of a planet/moon (any planet, though this has begun with me on Europa) well enough to fool the naked eye. I am not trying to plot a trajectory or something where small fractions of degree of error mean you miss your target. I would like to go there later, but this is a simple case right here.

Also, the utility is the scenario editor, thus I am constrained to work with lat/long inputs as I position the lights around the landing pad.

Finally

It should lastly be said that now I have tried to do something and it didn't work. Now I have the bone in my teeth and I don't want to let go. I have attempted to figure great circle distance and it did not work for me.

Can Anyone Fix This?

Here is what I have done via spreadsheet:

Using the Haversine formula 
a = sin²(Δlat/2) + cos(lat1).cos(lat2).sin²(Δlong/2)
c = 2.atan2(√a, √(1−a))
d = R.c

R (Earth mean radius) = 6371km

First case I tried was a 10 degree distance along the equator.

lat1 = lat2 = 0 deg = 0 rad
Δlong = 10 deg = 0.174533 rad

then
a = 0.007596
c = 2.967060
d = 18903.137530km

That is clearly an absurd result. The web site I sourced this from says that should be 1112km. What have I done wrong? I note that the website supplying this does not in fact actually use JavaScript as quoted to perform this calculation.

 

[agentgonzo]

You've probably already realized what I ran into. The base is at about 68 N and it appears the longitudinal distance equal to 0.006 degree is about 1/3 or 1/4 by eyeball of what it equals for the latitudinal distance. I got a nice ellipse with a North-South orientation. Suited the Arrow freighter since it was aligned that way, but was not what I wanted.

In your calculations you'll need to add a multiplicative constant of [math]cos(\theta)[/math] (where [math]\theta[/math] is your latitude) to make it a circle rather than an ellipse.

 

[Eccentric]

In your calculations you'll need to add a multiplicative constant of [math]cos(\theta)[/math] (where [math]\theta[/math] is your latitude) to make it a circle rather than an ellipse.

Perfect. Thank you sir. Exactly what I needed.

Of course, I still have that itch to scratch, so I'll be back to Fizyk's equations sometime soon anyway .

 

[Fizyk]

Yeah, my equations are a bit of an overkill and for small circles agentgonzo's method should be perfectly sufficient. I just always have that urge to find an exact solution if it is possible

 

[Notebook]

Its a bit off your needs, but some of this lot might help

http://www.orbiter-forum.com/showthread.php?t=21843

I'll have a go in Mathcad, and see how far I get with your and Fizyk's work.

N.

 

[Eccentric]

Its a bit off your needs, but some of this lot might help

http://www.orbiter-forum.com/showthread.php?t=21843

Turns out Urwumpe's little tidbit here...

most computer implementations take the data in the order (y, x)

...was exactly what my problem was. Reversed the order of the variables in the ATAN2 function and now I get the right numbers spitting out.

Cheers Notebook and Urwumpe.

 

[Notebook]

Glad you got it to work.

Functions taking paramaters in different order is something I wouldn't have thought of, and took me a long time to find out, if I didn't give up.
Seems a cheap trick if you ask me!

N.

 

[agentgonzo]

Turns out Urwumpe's little tidbit here...



...was exactly what my problem was. Reversed the order of the variables in the ATAN2 function and now I get the right numbers spitting out.

Cheers Notebook and Urwumpe.

I've fallen fowl of the atan ordering of parameters problem before. It's a pain.

 

[george7378]

You could always do it by using two formulas for x and y:

x = r cos a
y = r sin a

where x is the x-coordinate (assuming your centre is 0,0) and y is the y-coordinate, and a is the angle between the start point and the radius.

 

[Fizyk]

You could always do it by using two formulas for x and y:

x = r cos a
y = r sin a

where x is the x-coordinate (assuming your centre is 0,0) and y is the y-coordinate, and a is the angle between the start point and the radius.

Not on a sphere with x,y being latitude and longitude (well, you can use it, but it won't give you a circle).

 

[george7378]

Oh yes, but you could use it on a small scale where you can assume it to be flat right?

 

[Fizyk]

Almost. That was the initial problem - 1 degree difference in latitude doesn't give the same distance as 1 degree difference in longitude anywhere else than on the equator. But this is what agentgonzo's proposition fixes: you need to add cos(latitude) somewhere (either multiply latitude by it, or divide longitude) - then this method will give something close enough to a circle.

 

[Eccentric]

Almost. That was the initial problem - 1 degree difference in latitude doesn't give the same distance as 1 degree difference in longitude anywhere else than on the equator. But this is what agentgonzo's proposition fixes: you need to add cos(latitude) somewhere (either multiply latitude by it, or divide longitude) - then this method will give something close enough to a circle.

Indeed, it is close enough to fool the eye, which for this was truly all I wanted to achieve.