Finding max separation of planets

[Syntax]

I'm looking for some help/guidance on what is turning out to be a surprisingly challenging math problem:

I'm trying to find the maximum possible separation between any two given planets. I think it would be useful, or at least interesting, to know this kind of thing when constructing comms relay networks, etc.

I've got some idea how to proceed, but my efforts keep coming up dry (not least because of my merely intermediate math knowledge). Seems my problem needs new eyes and new ideas. I'm all ears! 

[Terwin]

Seems like the max possible separation would be  Orbital Radius(P1)+Orbital Radius(P2).

This would give the distance between the two planets when they are on opposite sides of the sun.

 

For elliptic orbits, use AP(P1)+AP(P2), this gives a maximum as it is the furthest each planet gets form the sun, even if their orbits make it impossible for them both to be at their furthest from the sun and on opposite sides of the sun at the same time.

[Shpaget]

But what about orbital inclination?

edit:

A quick googling found this:

https://pdfs.semanticscholar.org/0c8c/a8b1fe0eb533c0f12e45f3d689235fbb4957.pdf

Quote

Abstract. We are interested in the distance problem between two objects in three dimensional Euclidean space. There are many distance problems for various types of objects including line segments, boxes, polygons, circles, disks, etc. In this paper we present an iterative algorithm for finding the distance between two given ellipses. Numerical examples are given.

 

Edited November 20, 2018 by Shpaget

[Terwin]

1 hour ago, Shpaget said:

But what about orbital inclination?

Already part of the apoapsis.

apoapsis(Planet 1)+apoapsis(Planet 2) gives you the max possible distance between Planet 1 and Planet 2.

As you are looking at communications, that should be a reasonable estimate.

If you want the max distance for some moons, you just add that in:

apoapsis(Planet 1)+apoapsis(Moon 1)+apoapsis(Planet 2)+apoapsis(Moon 2)

This gives you the distance should the moons, their planets, and the sun all be in a straight-line with each of them being as far as possible form each of the others.

Orbital inclination can only make that distance shorter(Imagine a you have two short lines and one long line.  The short lines both touch the origin with one end, and the long line extends to touch the other ends of both short lines.  The longest you can possibly make the long line is short1+short2 by having both short lines pointing directly away form each other.  Any other angle, will allow the long line to get shorter than short1+short2)

[K^2]

1 hour ago, Terwin said:

Already part of the apoapsis.

apoapsis(Planet 1)+apoapsis(Planet 2) gives you the max possible distance between Planet 1 and Planet 2.

It's certainly an upper bound, but not the least upper bound. Most configurations are never going to reach that, even if many will get close to that at some point or another. As an extreme example, more plausible for comets than planets, imagine two orbits with similar inclination that are both highly elliptic, but share argument of periapsis. Their maximum separation will actually be between apoapsis of one and periapsis of another.

I agree that sum of the apoapsides is a good starting point. If your com system is capable of handling that, it will handle any distance that is going to arise. But if you're trying to cut corners and scrape by on absolute minimum, you'll have to do more math.

And unfortunately, I don't think there is a simple formula. I can write down a general equation for distance between two bodies as function of true anomalies, and I can even give condition for the extreme separation, but actually solving for it? Ugh! In contrast, this is a very simple problem for a computer. You're trying to find argmax of a two variable function which will have just a few local maxima, if any.

@Syntax I don't know how good you are with numerical methods / code. If you want, I can write you a script in language of your choice that takes orbital elements as inputs and spits out max separation and relevant parameters for it. If that does you any good, it shouldn't take me more than half an hour to write.

[Syntax]

1 hour ago, K^2 said:

 @Syntax I don't know how good you are with numerical methods / code. If you want, I can write you a script in language of your choice that takes orbital elements as inputs and spits out max separation and relevant parameters for it. If that does you any good, it shouldn't take me more than half an hour to write.

That would be incredible! I'm very much a beginner with code, but do have some limited experience with javascript, ruby, and (to an even lesser extent) python. If any of those work for you, I'd be "over the moon" just to have something to work from. I will say, though, that I'd need some clear commenting. All the internets to you for such a favor!

[Xd the great]

Use coordinate geometry.

[K^2]

I wrote it in ES6-ish JS, and I figured I might as well wrap it in HTML so that you have an immediate use case.The HTML portion of this is not compliant with anything, but it should work^TM on most browsers. The JS portion can be made a lot cleaner as well, and I didn't bother with variable names that much, but it's functional. Keep in mind that the two optimization functions I've written are not terribly generic. That said, this isn't a terribly generic kind of problem.

The input parameters are the standard orbital elements for each of the bodies. Semimajor axis can be specified in whatever units you want, so long as they're the same. The distance reported will be in the same units. AU would be a good choice for Solar System. Eccentricity is just a number in [0, 1) range. The three angles should be in radians. You can add conversion from degrees to radians inside the compute() function, where the values are extracted from input fields, if this makes your life simpler. I used radians because that's what Math library expects.

Anyways, just paste this mess into the editor and save it as an .html file. No error checking, so make sure all fields are filled in before you hit the button. Otherwise, should work.

Spoiler

 

<html>
<head>
<script language="JavaScript">
// Inverse square and inverse of golden ratio.
const phi1 = 0.38196601125;
const phi2 = 0.61803398875;

// Golden section search for a maximum.
function golden(x1, x2, f) {
	x0 = x2;
	x2 = x2 - x1;
	let v1 = f(x1 + phi1 * x2);
	let v2 = f(x1 + phi2 * x2);
	let ref = f(x1);
	for (let i = 0; i < 40; ++i) {
		if (v1 > v2) {
			v2 = v1;
			x2 *= phi2;
			v1 = f(x1 + phi1 * x2);
		} else {
			v1 = v2;
			x1 += phi1 * x2;
			x2 *= phi2;
			v2 = f(x1 + phi2 * x2);
		}
	}
	return x1 + 0.5 * x2;
}

// Rudimentary 2D maximization. Assumes nothing fancy going on.
function max2d(x1, x2, y1, y2, f) {
	let x = 0.5 * (x1 + x2);
	let y = 0.5 * (y1 + y2);
	let v = f(x, y);
	while (true) {
		y = golden(y1, y2, (z) => f(x, z));
		x = golden(x1, x2, (z) => f(z, y));
		const nv = f(x, y);
		if (Math.abs(nv - v) / (Math.abs(v) + Math.abs(nv)) < 1E-6) {
			return [nv, x, y];
		}
		v = nv;
	}
}

// Transforming X axis from 2D ellipse to 3D.
function transformX(i, O, w) {
	const x = Math.cos(w);
	const r = Math.sin(w);
	const y = Math.cos(i) * r;
	const z =-Math.sin(i) * r;
	const c = Math.cos(O);
	const s = Math.sin(O);
	return [x * c - y * s, x * s + y * c, z];
}

// Transforming Y axis from 2D ellipse to 3D.
function transformY(i, O, w) {
	const x =-Math.sin(w);
	const r = Math.cos(w);
	const y = Math.cos(i) * r;
	const z =-Math.sin(i) * r;
	const c = Math.cos(O);
	const s = Math.sin(O);
	return [x * c - y * s, x * s + y * c, z];
}

// 3D Position at true anomaly T.
function position(tx, ty, a, e, T) {
	const c = Math.cos(T);
	const r = a * (1 - e * e) / (1 + e * c);
	const x = r * c;
	const y = r * Math.sin(T);
	return [tx[0] * x + ty[0] * y, tx[1] * x + ty[1] * y, tx[2] * x + ty[2] * y];
}

// Eucleidan distance.
function distance2(p1, p2) {
	const x = p1[0] - p2[0];
	const y = p1[1] - p2[1];
	const z = p1[2] - p2[2];
	return x * x + y * y + z * z;
}

function plot(tx, ty, a, e, r, color = '#000') {
	const cvs = document.getElementById('canvas');
	const w = cvs.width;
	const h = cvs.height;
	const factor = Math.min(w, h) / (2.1 * r);

	let p = position(tx, ty, a, e, 0);

	const ctx = cvs.getContext('2d');
	ctx.beginPath();
	ctx.moveTo(w/2 + factor * p[0], h/2 - factor * p[1]);
	for (let i = 0; i < 2 * 6.3; i += 0.05) {
		p = position(tx, ty, a, e, i);
		ctx.lineTo(w/2 + factor * p[0], h/2 - factor * p[1]);
	}
	ctx.strokeStyle = color;
	ctx.stroke();
}

function draw(tx, ty, a, e, T, r, color = '#000') {
	const cvs = document.getElementById('canvas');
	const w = cvs.width;
	const h = cvs.height;
	const factor = Math.min(w, h) / (2.1 * r);

	let p = position(tx, ty, a, e, T);

	const ctx = cvs.getContext('2d');
	ctx.beginPath();
	ctx.arc(w/2 + factor * p[0], h/2 - factor * p[1], 5, 0, 2 * Math.PI);
	ctx.fillStyle = color;
	ctx.fill();
}

// Putting it all together.
function compute() {
	// Set up parameters.
	const i1 = parseFloat(document.getElementById("i1").value);
	const O1 = parseFloat(document.getElementById("O1").value);
	const w1 = parseFloat(document.getElementById("w1").value);
	const tx1 = transformX(i1, O1, w1);
	const ty1 = transformY(i1, O1, w1);

	const i2 = parseFloat(document.getElementById("i2").value);
	const O2 = parseFloat(document.getElementById("O2").value);
	const w2 = parseFloat(document.getElementById("w2").value);
	const tx2 = transformX(i2, O2, w2);
	const ty2 = transformY(i2, O2, w2);

	const a1 = parseFloat(document.getElementById("a1").value);
	const e1 = parseFloat(document.getElementById("e1").value);
	const a2 = parseFloat(document.getElementById("a2").value);
	const e2 = parseFloat(document.getElementById("e2").value);

	// Scaling for plotting.
	const r = Math.max(a1 * (1 + e1), a2 * (1 + e2));

	// This is the actual distance function.
	const f = (T1, T2) => distance2(position(tx1, ty1, a1, e1, T1), position(tx2, ty2, a2, e2, T2));

	// Compute distances in 4 quadrants.
	distances = []
	for (let i = 0; i < 6.3; i += 0.2) {
		for (let j = 0; j < 6.3; j += 0.2) {
			distances.push(max2d(i, i + 0.2, j, j + 0.2, f));
		}
	}

	// Final result.
	const res = distances.reduce((a, b) => a[0] > b[0] ? a : b, [0, 0, 0]);
	const d = Math.sqrt(res[0]);
	document.getElementById("max").innerHTML = d;
	document.getElementById("t1").innerHTML = res[1];
	document.getElementById("t2").innerHTML = res[2];

	// Center.
	draw([1, 0, 0], [0, 1, 0], 0, 0, 0, r);

	// Objects at max separation.
	draw(tx1, ty1, a1, e1, res[1], r, 'red');
	draw(tx2, ty2, a2, e2, res[2], r, 'blue');

	// Orbits.
	plot(tx1, ty1, a1, e1, r, 'red');
	plot(tx2, ty2, a2, e2, r, 'blue');
}
</script>
<style>
.grid {
	display: grid;
	grid-template-columns: repeat(4, auto) 1fr;
	grid-template-rows: repeat(10, auto) 1fr;
	justify-content: start;
	grid-gap: 5px;
}
.head { grid-column: any / span 2; }
.label { justify-self: end; }
</style>
</head>
<body>
<div class="grid">
<div style="grid-column: 1/span 2;">Body 1 orbital elements:</div>
<div style="grid-column: 3/span 2;">Body 2 orbital elements:</div>
<div style="grid-column: 5; grid-row: 1/span 11;"><canvas id="canvas" width=600 height=600/></div>
<div class="label">Semimajor axis:</div><div><input size=10 id="a1" value="68773.6"/></div>
<div class="label">Semimajor axis:</div><div><input size=10 id="a2" value="40839.3"/></div>
<div class="label">Eccentricity:</div><div><input size=10 id="e1" value="0.05"/></div>
<div class="label">Eccentricity:</div><div><input size=10 id="e2" value=" 0.145"/></div>
<div class="label">Inclination (rad):</div><div><input size=10 id="i1" value="0.02276"/></div>
<div class="label">Inclination (rad):</div><div><input size=10 id="i2" value="0.08727"/></div>
<div class="label">Ascending node (rad):</div><div><input size=10 id="O1" value="0.9076"/></div>
<div class="label">Ascending node (rad):</div><div><input size=10 id="O2" value="4.887"/></div>
<div class="label">Periapsis (rad):</div><div><input size=10 id="w1" value="0"/></div>
<div class="label">Periapsis (rad):</div><div><input size=10 id="w2" value="1.571"/></div>
<div style="grid-column: 1 / span 4;">
<input type="button" value="Compute" onclick="compute()"/></div>
<div class="label">Greatest separation:</div><div id="max" style="grid-column: 2/span 3;">&nbsp;</div>
<div class="label">True anomaly 1:</div><div id="t1" style="grid-column: 2/span 3;">&nbsp;</div>
<div class="label">True anomaly 2:</div><div id="t2" style="grid-column: 2/span 3;">&nbsp;</div>
</div>
</body>
</html>

 

 

 

Edit: The code's been updated to fix an issue with argument of periapsis, make it more resilient to local maxima, and to add visualization.

Edited November 29, 2018 by K^2

[Syntax]

On 11/21/2018 at 2:41 AM, K^2 said:

Anyways, just paste this mess into the editor and save it as an .html file. No error checking, so make sure all fields are filled in before you hit the button. Otherwise, should work.

This. Is. Incredible! You're my hero! Thank you!

Purely for curiosity's sake, do you think there would be a relatively simple way to determine the points the two bodies would have to be at in their orbits to achieve their max separation? Looking at your code, I see a position function devoted to "3D Position at true anomaly T." That looks like the ticket!

I'm having a hard time getting that value onto the browser page, though. I think it's because I haven't looked at JS in some time and am not at all familiar with the arrow function... what do I need to do here?

[K^2]

@Syntax max2d currently returns just the distance. You can modify that function to return x and y at the same time. These are true anomalies of first and second bodies respextively. You can feed each one in turn to the position function with corresponding parameters if you want the cartesian coordinates. So position(tx1, ty1, a1, e1, x) should do the trick for the first body.

[Syntax]

On 11/23/2018 at 12:32 AM, K^2 said:

You can feed each one in turn to the position function with corresponding parameters if you want the cartesian coordinates.

My cup overflows! Thank you!

I will try this out tonight. 

Am I correct in my understanding that the centre of the ellipse (and not one of the foci) is at the origin?

Edited November 24, 2018 by Syntax

[K^2]

On 11/23/2018 at 5:14 AM, Syntax said:

Am I correct in my understanding that the centre of the ellipse (and not one of the foci) is at the origin?

The focus is at the origin. The polar form used is r(θ) = a (1 - e²) / (1 + e cos(θ)). So r(θ) is minimal at θ = 0, corresponding to periapsis. This is then cast to Cartesian coordinates, matching θ = 0 to the positive x axis, and keeping the ellipse in the XY plane. The transform components tx and ty are then used to rotate the ellipse in 3D to match orbital elements, keeping the focus at the origin.

As an aside, you can quickly convince yourself that r(0) + r(π) = a (1 - e²) / (1 + e) + a (1 - e²) / (1 + e) = a (1 - e) + a (1 + e) = 2a. Or in other words, sum of periapsis and apoapsis is twice the semimajor axis.

[Syntax]

On 11/24/2018 at 7:58 PM, K^2 said:

The focus is at the origin. The polar form used is r(θ) = a (1 - e²) / (1 + e cos(θ)). So r(θ) is minimal at θ = 0, corresponding to periapsis. This is then cast to Cartesian coordinates, matching θ = 0 to the positive x axis, and keeping the ellipse in the XY plane. The transform components tx and ty are then used to rotate the ellipse in 3D to match orbital elements, keeping the focus at the origin.

@K^2 That's fantastic. Thanks for clarifying. I see it now.

Thing is... something's not adding up. I've made this model (linked in my signature) against which I'm checking the angles and distances this script is spitting out. I've been using Dres and Jool... and I'm getting a max distance that I can't seem to replicate no matter how hard I try (off by 3Gm!). Plus, the angles are just plain nonsensical. I don't know what is leading to the inconsistencies. I'm trying to get it to add up, but no luck yet.

If you get a chance to look at what I'm talking about, I'd be interested to hear your thoughts on a diagnosis.

Here's what I've done with the code you provided. I'll warn you: it's not pretty; I'm not good enough at this to make it streamlined.

 

Spoiler

<html>
<head>
<script language="JavaScript">
// Inverse square and inverse of golden ratio.
const phi1 = 0.38196601125;
const phi2 = 0.61803398875;
// One, three, and five quarters of pi.
const oqpi = 0.78539816339;
const tqpi = 2.35619449019;
const fqpi = 3.92699081699;

// Golden section search for a maximum.
function golden(x1, x2, f) {
	x0 = x2;
	x2 = x2 - x1;
	let v1 = f(x1 + phi1 * x2);
	let v2 = f(x1 + phi2 * x2);
	let ref = f(x1);
	for (i = 0; i < 30; ++i) {
		if (v1 > v2) {
			v2 = v1;
			x2 *= phi2;
			v1 = f(x1 + phi1 * x2);
		} else {
			v1 = v2;
			x1 += phi1 * x2;
			x2 *= phi2;
			v2 = f(x1 + phi2 * x2);
		}
	}
	return x1 + 0.5 * x2;
}

// Rudimentary 2D maximization. Assumes nothing fancy going on.
function max2d(x1, x2, y1, y2, f) {
	let x = 0.5 * (x1 + x2);
	let y = 0.5 * (y1 + y2);
	let v = f(x, y);
	while (true) {
		y = golden(y1, y2, (z) => f(x, z));
		x = golden(x1, x2, (z) => f(z, y));
		const nv = f(x, y);
		if (Math.abs(nv - v) / (Math.abs(v) + Math.abs(nv)) < 1E-6) {
          // modify return value
			return [nv, x, y];
		}
		v = nv;
		
	}
}

// Transforming X axis from 2D ellipse to 3D.
function transformX(i, O, w) {
	const x = Math.cos(w - O);
	const r = Math.sin(w - O);
	const y = Math.cos(i) * r;
	const z =-Math.sin(i) * r;
	const c = Math.cos(O);
	const s = Math.sin(O);
	return [x * c - y * s, x * s + y * c, z];
}

// Transforming Y axis from 2D ellipse to 3D.
function transformY(i, O, w) {
	const x =-Math.sin(w - O);
	const r = Math.cos(w - O);
	const y = Math.cos(i) * r;
	const z =-Math.sin(i) * r;
	const c = Math.cos(O);
	const s = Math.sin(O);
	return [x * c - y * s, x * s + y * c, z];
}

// 3D Position at true anomaly T.
function position(tx, ty, a, e, T) {
	const c = Math.cos(T);
	const r = a * (1 - e * e) / (1 + e * c);
	const x = r * c;
	const y = r * Math.sin(T);
	return [tx[0] * x + ty[0] * y, tx[1] * x + ty[1] * y, tx[2] * x + ty[2] * y];
}

// Eucleidan distance.
function distance(p1, p2) {
	const x = p1[0] - p2[0];
	const y = p1[1] - p2[1];
	const z = p1[2] - p2[2];
	return Math.sqrt(x * x + y * y + z * z);
}

// Putting it all together.
function compute() {
	// Set up parameters.
	const i1 = parseFloat(document.getElementById("i1").value);
	const O1 = parseFloat(document.getElementById("O1").value);
	const w1 = parseFloat(document.getElementById("w1").value);
	const tx1 = transformX(i1, O1, w1);
	const ty1 = transformY(i1, O1, w1);

	const i2 = parseFloat(document.getElementById("i2").value);
	const O2 = parseFloat(document.getElementById("O2").value);
	const w2 = parseFloat(document.getElementById("w2").value);
	const tx2 = transformX(i2, O2, w2);
	const ty2 = transformY(i2, O2, w2);

	const a1 = parseFloat(document.getElementById("a1").value);
	const e1 = parseFloat(document.getElementById("e1").value);
	const a2 = parseFloat(document.getElementById("a2").value);
	const e2 = parseFloat(document.getElementById("e2").value);

	// This is the actual distance function.
	const f = (T1, T2) => distance(position(tx1, ty1, a1, e1, T1), position(tx2, ty2, a2, e2, T2));

	// Compute distances in 4 quadrants.
	// modify this section to handle modified max2d output
	const d0 = max2d(oqpi, tqpi, oqpi, tqpi, f);
		const d00 = d0[0];
		
	const d1 = max2d(oqpi, tqpi, tqpi, fqpi, f);
		const d01 = d1[0];
		
	const d2 = max2d(tqpi, fqpi, oqpi, tqpi, f);
		const d10 = d2[0];
		
	const d3 = max2d(tqpi, fqpi, tqpi, fqpi, f);
		const d11 = d3[0];

	// Final result.
	const d = Math.max(d00, d01, d10, d11);
	document.getElementById("max").innerHTML = d;
	
	// Modify to show position result.
	const arrdarr = [d0, d1, d2, d3];
	const darr = [d00, d01, d10, d11];
	const dind = darr.indexOf(Math.max(...darr));
	
	const apearr = arrdarr[dind];
	const tran1 = apearr[1];
	const tran2 = apearr[2];
	
	document.getElementById("T1").innerHTML = tran1;
	document.getElementById("T2").innerHTML = tran2;
	
	// True anomaly to eccentric anomaly
	document.getElementById("ecc1").innerHTML = Math.PI-(Math.asin((Math.sqrt(1-e1*e1)*Math.sin(tran1))/(1+e1*Math.cos(tran1))));
	document.getElementById("ecc2").innerHTML = Math.PI-(Math.asin((Math.sqrt(1-e2*e2)*Math.sin(tran2))/(1+e2*Math.cos(tran2))));
	
	// Modify to show body coordinates.
	const pt1 = position(tx1, ty1, a1, e1, tran1);
	const pt2 = position(tx2, ty2, a2, e2, tran2);
	document.getElementById("pos1").innerHTML = pt1;
	document.getElementById("pos2").innerHTML = pt2;
}
</script>
<style>
.grid { display: grid; grid-template-columns: repeat(4, auto); justify-content: start; grid-gap: 5px; }
.head { grid-column: any / span 2; }
.label { justify-self: end; }
</style>
</head>
<body>
<div class="grid">
<div style="grid-column: 1/span 2;">Body 1 orbital elements:</div>
<div style="grid-column: 3/span 2;">Body 2 orbital elements:</div>
<div class="label">Semimajor axis:</div><div><input size=10 id="a1"/></div>
<div class="label">Semimajor axis:</div><div><input size=10 id="a2"/></div>
<div class="label">Eccentricity:</div><div><input size=10 id="e1"/></div>
<div class="label">Eccentricity:</div><div><input size=10 id="e2"/></div>
<div class="label">Inclination (rad):</div><div><input size=10 id="i1"/></div>
<div class="label">Inclination (rad):</div><div><input size=10 id="i2"/></div>
<div class="label">Ascending node (rad):</div><div><input size=10 id="O1"/></div>
<div class="label">Ascending node (rad):</div><div><input size=10 id="O2"/></div>
<div class="label">Periapsis (rad):</div><div><input size=10 id="w1"/></div>
<div class="label">Periapsis (rad):</div><div><input size=10 id="w2"/></div>
<div style="grid-column: 1 / span 4;">
<input type="button" value="Compute" onclick="compute()"/></div>
</div>
<br>
<div class="label">Greatest separation:</div><div id="max"> </div>
<div class="label">True Anomaly 1:</div><div id="T1"> </div>
<div class="label">True Anomaly 2:</div><div id="T2"> </div>
<div class="label">Ecc. Anomaly 1:</div><div id="ecc1"> </div>
<div class="label">Ecc. Anomaly 2:</div><div id="ecc2"> </div>
<div class="label">Position 1:</div><div id="pos1"> </div>
<div class="label">Position 2:</div><div id="pos2"> </div>
</body>
</html>

 

 

Edited November 26, 2018 by Syntax
typos

[K^2]

@Syntax Hm. I'll take a thorough look in the evening. There's always a chance that I made a mistake somewhere. Can you give me the parameters you've been using, just to make sure I'm reproducing exactly the same thing?

I was kind of considering adding a graphical plot of the orbits and the solution, but it seemed way beyond the scope. With debugging considerations in mind, however, it might be worth the extra effort.

Edited November 26, 2018 by K^2

[Syntax]

42 minutes ago, K^2 said:

Can you give me the parameters you've been using, just to make sure I'm reproducing exactly the same thing?

Sure can!

Jool
SMA: 68773.56032
e: 0.05
Inc: 0.022759093795546
LAN: 0.907571211037052
APe: 0

Dres
SMA: 40839.348203
e: 0.145
Inc: 0.0872664625997166
LAN: 4.88692190558413
APe: 1.57079632679490

Thanks so much for all your help! This is amazing.

[K^2]

@Syntax I'm still going to take a closer look at the code, but I just made a quick plot, and it looks believable, at least.

Jool position: [-34184.36954230997, 61424.45691029616, -1474.0485511224515]

Jool true anomaly: 2.078729696847789

Dres position: [17768.806076568617, -42588.466156726034, -883.9372007234914]

Dres true anomaly: 3.537671968402128

Note that true anomaly is angle from periapsis, not from datum. So for Dres, you need to add 90^o.

[Syntax]

41 minutes ago, K^2 said:

I'm still going to take a closer look at the code, but I just made a quick plot, and it looks believable, at least.

...

Note that true anomaly is angle from periapsis, not from datum. So for Dres, you need to add 90^o.

@K^2That does look more believable. Seems likely I'm having issues on my end. Not sure why, though. I will also continue looking.

One thing I do notice (and this could help explain our discrepancies), is that on the graph you provided, it looks like Dres' periapsis lies on the positive y-axis, and Jool's orbit swings away from it to the upper left. When I look at my model (and indeed at the map screen in KSP) in that orientation, Jool's orbit swings away from Dres' to the upper right when looking straight down at it.

Edited November 27, 2018 by Syntax
tagged

[K^2]

@Syntax Good catch on periapsis. There are several ways that combination of three angles can describe a 3D rotation, and the one I used is not the standard one. So while orbital plane is correct, location of periapsis on it is not. The fix will be updated TransformX and TransformY functions that interpret orbital parameters in the standard way. I'll edit the post with the code with a fix and post any further updates.

Edit: Actually... Try changing (w - O) to just (w) in first two lines of both the transformX and transformY functions. That might actually fix the whole thing, but I'm going to have to spend a little bit more time with it to be absolutely sure that's correct.

Late edit: Last two days have been trapped at work until late. Tonight's looking much better. Hopefully will get around to all of that.

Edited November 28, 2018 by K^2

[Syntax]

29 minutes ago, K^2 said:

Try changing (w - O) to just (w) in first two lines of both the transformX and transformY functions.

@K^2 This gives a max separation a few Gm less than what I've been able to find just through trial and error.

[K^2]

On 11/26/2018 at 6:17 PM, Syntax said:

@K^2 This gives a max separation a few Gm less than what I've been able to find just through trial and error.

Yeah, there was also a local maximum that I did not expect. Like I mentioned before, max2d is a VERY simple algorithm. So I broke the orbits into more chunks, and now it seems to be good. I've updated the code with the fix to transformX/transformY, made it consider more segments of the orbit to avoid local maximum, and added plotting of the results. The updated code took place of the original further up.

One more thing I've changed is that I replaced distance with distance2, which returns square of the distance. Square root is a convex function, there's really no reason to crank it out every step of evaluation just to find the maximum. Square of the distance will always be highest in the same configuration where the distance is highest. So instead, I take the square root at the very end.

Oh, and I pre-filled the form with the numbers above. It just made testing so much quicker. Don't know why I didn't do that from the start.

[Syntax]

@K^2 Thank you SO MUCH for what you've done here. You've single-handedly brought me to within a hair's breadth of the solution to this problem that has been bugging me for months!

I say "within a hair's breadth" because, oddly, when I put the true anomalies output by your algorithm into my model, I get a distance of 113191.131383 Mm, a difference of 169.675468 Mm from your algorithm's output (113360.806852). But, when I play around with fine adjustments to those body positions in my model, I can maximize the distance calculated in the model to within <10m of the algorithm output! So I just need to figure out why the true anomalies I'm finding are around 0.05 rad off from yours. I really believe my model is accurate... but perhaps more checks on that are in order. Whether the issue lies with my model or your code, I will continue looking into this and keep you updated.

I'm sorry I needed so much hand holding on this... it's that I don't have a good grasp on the math used to transform the ellipses from 2D into 3D. Any resources you could recommend for me to read/watch to deepen that skill set?

[mikegarrison]

4 hours ago, Syntax said:

@K^2 Thank you SO MUCH for what you've done here. You've single-handedly brought me to within a hair's breadth of the solution to this problem that has been bugging me for months!

I say "within a hair's breadth" because, oddly, when I put the true anomalies output by your algorithm into my model, I get a distance of 113191.131383 Mm, a difference of 169.675468 Mm from your algorithm's output (113360.806852). But, when I play around with fine adjustments to those body positions in my model, I can maximize the distance calculated in the model to within <10m of the algorithm output! So I just need to figure out why the true anomalies I'm finding are around 0.05 rad off from yours. I really believe my model is accurate... but perhaps more checks on that are in order. Whether the issue lies with my model or your code, I will continue looking into this and keep you updated.

I'm sorry I needed so much hand holding on this... it's that I don't have a good grasp on the math used to transform the ellipses from 2D into 3D. Any resources you could recommend for me to read/watch to deepen that skill set?

Is there any chance you are looking at numerical artifacts? How good is javascript at math precision?

Possibly it needs something like this? http://mathjs.org/docs/datatypes/bignumbers.html

Edited November 30, 2018 by mikegarrison

[K^2]

@Syntax Well, it depends on how hard-core you want to get into it. If you really want to understand coordinate transformations, you have to approach it from perspective of Linear Algebra (or Differential Geometry, if you're planning to dive DEEP, but if you stick to Euclidean spaces, that's overkill.) I don't know how good your background is, and how much time you want to spend on it. Ideally, you want to either take a course, or at least read a textbook on Linear Algebra, but that's a serious time investment. As far as directly relevant materials, I did a quick search, and Texas A&M has decent slides up from a course a few years back. The three most relevant chapters, leading up to the actual coordinate transformations, are as follow:

http://www.math.tamu.edu/~yvorobet/MATH304-2011C/Lect2-03web.pdf

http://www.math.tamu.edu/~yvorobet/MATH304-504/Lect2-06web.pdf

http://www.math.tamu.edu/~yvorobet/MATH304-2011C/Lect2-10web.pdf

If you can follow these slides from start to end, you basically know as much about coordinate transformations as anyone, (again, so long as we stick to neat, linear transformations, such as rotations in Euclidean space). Naturally, you can dig up slides for earlier lectures to fill in gaps if you need to.

The good thing about a course in Linear Algebra is that it requires very little foreknowledge. It's not like being dropped into, say, a Topology or Partial Differential Equations course without pre-reqs. So long as you understand systems of linear equations and matrices, you're basically set to go through it. And while it's not strictly necessary just to get the mechanical skill of doing coordinate transformations, if you want it to MAKE SENSE, this is the best way to get there.

For a much simpler overview, you can look at some practical tutorials for video games. And it might be way easier to get a hang of it in 2D, since you can draw neat little pictures for yourself. Again, just following a quick search, this one looks like it covers the basics:

https://www.alanzucconi.com/2016/02/10/tranfsormation-matrix/

This sort of thing won't give you the deep insight into why it works, and prepare you for more general cases, but if learning Linear Algebra is beyond your time availability, looking for tutorials like this one might be a good start.

@mikegarrison It should be double, but even single precision would be fine, as that's good to one part per million. Big numbers is waaaaay overkill. That said, Math.js is a fine library.

Edited November 30, 2018 by K^2

[Syntax]

On 11/30/2018 at 1:22 AM, K^2 said:

@Syntax Well, it depends on how hard-core you want to get into it.

@K^2Really, just enough to work out what's causing that 0.05 rad discrepancy. I'm happy to research and learn, but it can be hard to know where to look, and my formal mathematics training ended with high school calculus.

Quote

Texas A&M has decent slides up from a course a few years back. The three most relevant chapters, leading up to the actual coordinate transformations, are as follow:

http://www.math.tamu.edu/~yvorobet/MATH304-2011C/Lect2-03web.pdf

http://www.math.tamu.edu/~yvorobet/MATH304-504/Lect2-06web.pdf

http://www.math.tamu.edu/~yvorobet/MATH304-2011C/Lect2-10web.pdfIf you can follow these slides from start to end, you basically know as much about coordinate transformations as anyone, (again, so long as we stick to neat, linear transformations, such as rotations in Euclidean space).

These are definitely well out of my league. Thank you for the resource, though. Who knows, I may end up to coming back to them.

Quote

For a much simpler overview, you can look at some practical tutorials for video games. And it might be way easier to get a hang of it in 2D, since you can draw neat little pictures for yourself. Again, just following a quick search, this one looks like it covers the basics:

https://www.alanzucconi.com/2016/02/10/tranfsormation-matrix/

This sort of thing won't give you the deep insight into why it works, and prepare you for more general cases, but if learning Linear Algebra is beyond your time availability, looking for tutorials like this one might be a good start.

This helped a lot, and I'm going to try to find more stuff like this. Just about the only thing left in your JS that I don't understand is the transformX and transformY functions, and I'm hoping that with that last puzzle piece I can close the book on this.

Edited December 2, 2018 by Syntax

[K^2]

6 hours ago, Syntax said:

Just about the only thing left in your JS that I don't understand is the transformX and transformY functions, and I'm hoping that with that last puzzle piece I can close the book on this.

These are the rows of the rotation matrix used to transform from 2D to 3D. So if you have initial vector v = {x, y, 0}, we can write down the transformed vector in 3D as v' = v R, where R is the rotation matrix and v' is the vector after transformation. We can then alternatively write it out as v' = x X + y Y. Here, X and Y are first two rows of R and exactly what transformX and transformY compute. We don't need transformZ because z = 0 before transformation. Alternatively, you can think of them as directions of the X axis and Y axis of the original, 2D coordinate system, if we apply transformation to them. Since the equation for ellipse I've used has periapsis on positive X axis, transformX can also be viewed as vector in direction of periapsis, and transformY gives a vector orthogonal to it, but also in the orbital plain.