我正在开发一款游戏,它的主用户界面是自上而下显示太阳系的(地球的北极大致是“向上”的)。我需要展示一个物体在行星之间的半现实的旅行路径。我一直在用兰伯特问题来确定物体的初速度矢量,我从中得到的结果似乎对我的游戏是合理的。
我已经试着解决这个问题几个星期了,甚至试着用ChatGPT-4来帮助解决我不理解的数学问题。我已经接近了几次,但似乎从来没有能够降落在适当的解决方案。
我对a
(椭圆的半长轴)的计算有时会出现负值,我知道这是不正确的。
我考虑的一个原因是,我的游戏(像大多数游戏一样)使用左手坐标系,但大多数数学都期望右手坐标系。(在我的游戏中,+y是下,在大多数数学中+y是上)。我已经处理过这个问题,每当我对向量进行计算时,我都会反转我的y值,并确保在我显示到屏幕上的时候,结果都是反转的。我也试过反转整个坐标系,但结果是相似的。
下面是TS实现的逻辑生成我的当前椭圆。
import { Constants } from './constants';
import Vector from './vector';
import * as math from 'mathjs';
const solver = require('lambert-orbit');
export default class LambertSolver {
public static async generatePath(
r0: Vector, // Initial position vector (km)
r1: Vector, // Final position vector (km)
deltat: number, // Time of flight (s)
solarMass: number, // Solar mass (kg)
pathResolution: number = 200
): Promise<Vector[]> {
const mu = solarMass * Constants.GRAVITY_KM;
const prograde = true;
const lowPath = true;
const vals = await solver(
mu,
r0.invertY().toArray(),
r1.invertY().toArray(),
deltat,
0,
prograde,
lowPath
);
const v1 = Vector.fromArray(vals[0]);
const [a] = this.orbitalParameters(r0.invertY(), v1, mu);
// const v1Norm = v1.magnitude();
// const r1Norm = r1.magnitude();
// const a = -mu / (v1Norm * v1Norm - (2 * mu) / r1Norm);
const secondFoci = this.findSecondFoci(r0.invertY(), v1, a, mu);
// const [secondFoci2, secondFoci3] = this.findSecondFociAlternative(
// r0.invertY(),
// r1.invertY(),
// a
// );
const path = this.ellipsePath(
new Vector(0, 0),
secondFoci,
a / 0.9,
pathResolution
);
return [
...path.map((x) => x.invertY()),
// ...this.ellipsePath(new Vector(0, 0), secondFoci2, a, pathResolution).map(
// (x) => x.invertY()
// ),
// ...this.ellipsePath(new Vector(0, 0), secondFoci3, a, pathResolution).map(
// (x) => x.invertY()
// ),
r0,
r1,
secondFoci.invertY(),
// secondFoci2,
// secondFoci3,
new Vector(0, 0),
r0.add(v1.invertY().multiply(2000000)),
r0,
];
}
private static orbitalParameters(
r1: Vector,
v1: Vector,
mu: number
): [number, number, number, number, number, number] {
const R = r1.magnitude();
const V = v1.magnitude();
const a = 1 / (2 / R - (V * V) / mu);
const Hv = r1.cross(v1);
const H = Hv.magnitude();
const p = (H * H) / mu;
const q = r1.dot(v1);
const e = Math.sqrt(1 - p / a);
const i = Math.acos(Hv.z / H);
let omega = 0;
if (i != 0) {
omega = Math.atan2(Hv.x, -Hv.y);
}
const Cw = (r1.x * Math.cos(omega) + r1.y * Math.sin(omega)) / R;
let Sw = 0;
if (i == 0 || i == Math.PI) {
Sw = (r1.y * Math.cos(omega) - r1.x * Math.sin(omega)) / R;
} else {
Sw = r1.z / (R * Math.sin(i));
}
const TAx = (H * H) / (R * mu) - 1;
const TAy = (H * q) / (R * mu);
const v = Math.atan2(TAy, TAx);
var w = Math.atan2(Sw, Cw) - v;
return [a, e, i, omega, w, v];
}
private static findSecondFoci(
r1: Vector,
v1: Vector,
a: number,
mu: number
): Vector {
// Compute the specific angular momentum vector
const h = r1.cross(v1);
// Compute the eccentricity vector
const epsilon = v1.magnitude() ** 2 / 2 - mu / r1.magnitude();
const e_vec = r1
.multiply(epsilon)
.subtract(v1.multiply(r1.dot(v1)))
.multiply(1 / mu);
// Compute the eccentricity (magnitude of the eccentricity vector)
const e = e_vec.magnitude();
// Compute the distance from the center of the ellipse to each focus
const c = a * e;
// Compute the center of the ellipse
const firstFocus = r1;
const center = firstFocus.multiply(1 - e);
// Find the direction vector from the first focus to the second focus
const direction = e_vec.normalize();
// Compute the position of the second focus
const secondFocus = center.add(direction.multiply(c));
return secondFocus;
}
private static ellipsePath(
focalPoint1: Vector,
focalPoint2: Vector,
a: number,
numberOfPoints: number = 50
): Vector[] {
// Calculate the distance between the two focal points
const focalDistance = focalPoint1.distanceTo(focalPoint2);
// Calculate the semi-minor axis (b)
const c = focalDistance / 2;
const b = Math.sqrt(a * a - c * c);
// Calculate the center of the ellipse
const centerX = (focalPoint1.x + focalPoint2.x) / 2;
const centerY = (focalPoint1.y + focalPoint2.y) / 2;
const center = new Vector(centerX, centerY);
// Calculate the angle of the line connecting the focal points with the positive x-axis
const angle = Math.atan2(
focalPoint2.y - focalPoint1.y,
focalPoint2.x - focalPoint1.x
);
// Generate ellipse points in polar coordinates and convert them to Cartesian coordinates
const path: Vector[] = [];
for (let i = 0; i < numberOfPoints; i++) {
const theta = (2 * Math.PI * i) / numberOfPoints;
// Ellipse equation in polar coordinates: r(theta) = (a * b) / sqrt((b * cos(theta))^2 + (a * sin(theta))^2)
const r =
(a * b) /
Math.sqrt(
Math.pow(b * Math.cos(theta), 2) + Math.pow(a * Math.sin(theta), 2)
);
// Convert polar coordinates to Cartesian coordinates
const x = r * Math.cos(theta);
const y = r * Math.sin(theta);
// Rotate the point by the angle and translate it to the center
const rotatedX = x * Math.cos(angle) - y * Math.sin(angle);
const rotatedY = x * Math.sin(angle) + y * Math.cos(angle);
const translatedX = rotatedX + centerX;
const translatedY = rotatedY + centerY;
path.push(new Vector(translatedX, translatedY));
}
return path;
}
}
在我自己的代码中有了几个版本的Lambert问题求解器(AI辅助)之后,我决定使用lambert-orbit包,它似乎正是我所需要的,并添加了一些检查和其他逻辑,我稍后会需要。
我的vector类实现如下:
export default class Vector {
public x: number;
public y: number;
public z: number;
constructor(x: number, y: number, z: number = 0) {
this.x = x;
this.y = y;
this.z = z;
}
// Add another vector to this one
public add(v: Vector): Vector {
return new Vector(this.x + v.x, this.y + v.y, this.z + v.z);
}
// Subtract another vector from this one
public subtract(v: Vector): Vector {
return new Vector(this.x - v.x, this.y - v.y, this.z - v.z);
}
// Multiply the vector by a scalar
public multiply(scalar: number): Vector {
return new Vector(this.x * scalar, this.y * scalar, this.z * scalar);
}
// Divide the vector by a scalar
public divide(scalar: number): Vector {
return new Vector(this.x / scalar, this.y / scalar, this.z / scalar);
}
// Calculate the magnitude of the vector
public magnitude(): number {
return Math.sqrt(this.x * this.x + this.y * this.y + this.z * this.z);
}
// Calculate the dot product of two vectors
public dot(v: Vector): number {
return this.x * v.x + this.y * v.y + this.z * v.z;
}
public distanceTo(other: Vector): number {
const dx = this.x - other.x;
const dy = this.y - other.y;
const dz = this.z - other.z;
return Math.sqrt(dx * dx + dy * dy + dz * dz);
}
// Copy the components of another vector into this one
public copy(v: Vector): Vector {
this.x = v.x;
this.y = v.y;
this.z = v.z;
return this;
}
public scale(scalar: number): Vector {
return new Vector(this.x * scalar, this.y * scalar, this.z * scalar);
}
public descale(scalar: number): Vector {
return new Vector(this.x / scalar, this.y / scalar, this.z / scalar);
}
public unit(): Vector {
const magnitude = this.magnitude();
if (magnitude === 0) {
throw new Error('Cannot compute the unit vector of a zero vector.');
}
return new Vector(
this.x / magnitude,
this.y / magnitude,
this.z / magnitude
);
}
public absolute(): Vector {
return new Vector(Math.abs(this.x), Math.abs(this.y), Math.abs(this.z));
}
// Create a new vector from this one
public clone(): Vector {
return new Vector(this.x, this.y, this.z);
}
public toArray(): number[] {
return [this.x, this.y, this.z];
}
public static fromArray(array: number[]): Vector {
return new Vector(array[0], array[1], array[2]);
}
public angleTo(other: Vector): number {
const dx = other.x - this.x;
const dy = other.y - this.y;
const dz = other.z - this.z;
const d = Math.sqrt(dx * dx + dy * dy + dz * dz);
return Math.acos(this.dot(other) / (this.magnitude() * d));
}
public invert(): Vector {
return new Vector(-this.x, -this.y, -this.z);
}
public invertY(): Vector {
return new Vector(this.x, -this.y, this.z);
}
public invertX(): Vector {
return new Vector(-this.x, this.y, this.z);
}
public invertZ(): Vector {
return new Vector(this.x, this.y, -this.z);
}
public angle(): number {
return Math.atan2(this.y, this.x);
}
public normalize(): Vector {
const norm = this.magnitude();
return new Vector(this.x / norm, this.y / norm, this.z / norm);
}
public cross(other: Vector): Vector {
const x = this.y * other.z - this.z * other.y;
const y = this.z * other.x - this.x * other.z;
const z = this.x * other.y - this.y * other.x;
return new Vector(x, y, z);
}
}
举个例子,如果我将下面的内容传递给generatePath
方法:r0:{x: 126619585.3943946, y: -79668083.39934838}
r1:{x: 78468331.12907715, y: -214058085.17912537}
增量:7776000
太阳能质量:1.989e+30
我得到一个速度矢量:{x: 4.544347796414998, y: 28.198168671041273}
我的orbitalParameters
方法返回138428143.17335743
的a
值
我的findSecondFoci
方法返回{x: -82494511.2804895, y: -126466117.64851947}
这将导致下面的屏幕截图。抱歉,有点乱,我刚开始在绘制的路径上添加点以进行调试。屏幕截图具有在游戏空间中可视化的每个点的描述。
正如您所看到的,椭圆不仅没有接近于截取r 0和r1位置,而且似乎在错误的方向上。
有趣的是,在看到这个之后,我又截取了另一个屏幕截图,显示了路径的X值反转而不是Y,椭圆的尖端非常接近目标目的地。我正在做一些测试,看看这是否只是一个巧合。如果这不是巧合,我只需要弄清楚为什么我的短半轴这么大。
1条答案
按热度按时间nx7onnlm1#
终于找到了我一直在寻找的解决方案!
我偶然发现了这个图形,显示了椭圆的可视化。
从That中我生成了一些代码来找到第二个焦点的两个可能的位置,通过生成一个圆圈来表示
2a - r
和2a - r0
,围绕它们各自的点。它们相交的地方,是我可能的焦点。代码如下:
导致了这些省略号包括我需要的那个它也和初始速度矢量完美匹配,进一步证实了这是我需要的结果。