typescript 如何根据初始状态值(r0、v0、deltat和mu)正确计算说明空间中物体轨道的椭圆

hc8w905p  于 2023-05-01  发布在  TypeScript
关注(0)|答案(1)|浏览(78)

我正在开发一款游戏,它的主用户界面是自上而下显示太阳系的(地球的北极大致是“向上”的)。我需要展示一个物体在行星之间的半现实的旅行路径。我一直在用兰伯特问题来确定物体的初速度矢量,我从中得到的结果似乎对我的游戏是合理的。
我已经试着解决这个问题几个星期了,甚至试着用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.17335743a
我的findSecondFoci方法返回{x: -82494511.2804895, y: -126466117.64851947}
这将导致下面的屏幕截图。抱歉,有点乱,我刚开始在绘制的路径上添加点以进行调试。屏幕截图具有在游戏空间中可视化的每个点的描述。

正如您所看到的,椭圆不仅没有接近于截取r 0和r1位置,而且似乎在错误的方向上。
有趣的是,在看到这个之后,我又截取了另一个屏幕截图,显示了路径的X值反转而不是Y,椭圆的尖端非常接近目标目的地。我正在做一些测试,看看这是否只是一个巧合。如果这不是巧合,我只需要弄清楚为什么我的短半轴这么大。

nx7onnlm

nx7onnlm1#

终于找到了我一直在寻找的解决方案!
我偶然发现了这个图形,显示了椭圆的可视化。

从That中我生成了一些代码来找到第二个焦点的两个可能的位置,通过生成一个圆圈来表示2a - r2a - r0,围绕它们各自的点。它们相交的地方,是我可能的焦点。
代码如下:

private static circleCircleIntersection(
    c1: Vector,
    r1: number,
    c2: Vector,
    r2: number
  ): [Vector, Vector] | null {
    const d = c2.subtract(c1).magnitude();
    if (d > r1 + r2 || d < Math.abs(r1 - r2)) {
      return null;
    }

    const a = (r1 * r1 - r2 * r2 + d * d) / (2 * d);
    const h = Math.sqrt(r1 * r1 - a * a);
    const p = c1.add(c2.subtract(c1).normalize().multiply(a));

    const x1 = p.x + (h * (c2.y - c1.y)) / d;
    const y1 = p.y - (h * (c2.x - c1.x)) / d;

    const x2 = p.x - (h * (c2.y - c1.y)) / d;
    const y2 = p.y + (h * (c2.x - c1.x)) / d;

    return [new Vector(x1, y1), new Vector(x2, y2)];
  }

导致了这些省略号包括我需要的那个它也和初始速度矢量完美匹配,进一步证实了这是我需要的结果。

相关问题