Unity轨道预测器仅在特定起始位置生效的问题求助

• 仅当卫星从行星右侧（坐标(1200, 0, 0)）起始时，轨道预测准确；从Z/Y轴任一坐标起始时，预测轨道完全错误、反向或无法绘制
• 游戏运行中按下空格键改变轨道后，同样出现预测错误
• 0倾角轨道因升交点赤经方程未做特殊处理，导致轨道无法绘制

环境信息 #

• Unity脚本，挂载于球体游戏对象（场景中唯一激活脚本）
• 引力常数G = 0.005，行星质量1e+09

原代码 #

using NUnit.Framework;
using UnityEngine;
using System.Collections.Generic;
using System;
using Unity.VisualScripting;
using UnityEngine.UIElements;
using JetBrains.Annotations;
using System.Collections;
using Unity.Mathematics.Geometry;

public class Physics : MonoBehaviour
{
    public Material shipTexture;
    public Material trailTexture;
    public Vector3 shipCoords;
    public Vector3 shipThrust;
    GameObject central;
    Rigidbody planetMass;
    float size;
    public float gravity;
    public float mass;

    public class Ship
    {
        public string name;
        public float xPos;
        public float yPos;
        public float zPos;
        public float xVel;
        public float yVel;
        public float zVel;
        public Vector3 spam;
        public Vector3 novec;
        public float std;
        public Vector3 ecvec;
        public float semaj;
        public float smen;
        public float inclination;
        public float rAscensionNode;
        public float argumentOfPeriapsis;
        public float trueAnomaly;
        public float gnd;
        public float dist;
        public Vector3 point;
        public GameObject body;
        public Rigidbody mass;
        public Transform pos;
        public MeshFilter mFilter;
        public Material mats;

        public Ship(string shipName, Vector3 coords)
        {
            name = shipName;
            xPos = coords.x;
            yPos = coords.y;
            zPos = coords.z;
            dist = MathF.Sqrt(MathF.Pow(xPos, 2) + MathF.Pow(yPos, 2) + MathF.Pow(zPos, 2));
            body = new GameObject(name);
            mass = body.AddComponent<Rigidbody>();
            mass.mass = 10;
            mass.useGravity = false;
            mass.angularDamping = 0;
            mass.AddRelativeTorque(new Vector3(-0.3f, 1, 0.5f), ForceMode.Impulse);
            pos = body.GetComponent<Transform>();
            pos.localScale = new Vector3(10, 10, 10);
            mFilter = body.AddComponent<MeshFilter>();
            GameObject go = GameObject.CreatePrimitive(PrimitiveType.Cube);
            mFilter.mesh = go.GetComponent<MeshFilter>().mesh;
            mFilter.gameObject.AddComponent<MeshRenderer>();
            Destroy(go);
            pos.position = new Vector3(xPos, yPos, zPos);
        }
    }

    public class Equations
    {
        public static float Gravity(Transform m2Loc, Rigidbody m1Mass, Rigidbody m2Mass, float gravityConstant)
        {
            return (m1Mass.mass * m2Mass.mass) * gravityConstant / MathF.Pow(m2Loc.position.magnitude, 3);
        }

        public static Vector3 SpecificAngularMomentum(Vector3 r, Vector3 v) // h
        {
            return Vector3.Cross(r, v);
        }

        public static Vector3 NodeVector(Vector3 SAM) // n
        {
            return Vector3.Cross(new Vector3(0, -1, 0), SAM);
        }

        public static float StandardGravParam(float gConst, float earthMass, float shipMass)
        {
            return gConst * (earthMass + shipMass);
        }

        public static Vector3 EccentricityVector(Vector3 r, Vector3 v, Vector3 h, float u) // e
        {
            return ((MathF.Pow(v.magnitude, 2) - u / r.magnitude) * r - (Vector3.Dot(r, v) * v)) / u;
        }

        public static float SpecialMechanicalEnergy(Vector3 v, Vector3 r, float u) // E
        {
            return (MathF.Pow(v.magnitude, 2) / 2) - (u / r.magnitude);            
        }

        public static float Momentum(float a, Vector3 e) // p
        {
            return a * (1 - MathF.Pow(e.magnitude, 2));
        }
        
        public static float Inclination(Vector3 h) // i
        {
            return MathF.Acos(-1 * h.y / h.magnitude);
        }

        public static float RightAscendingNode(Vector3 n) // Ω
        {
            return MathF.Acos(n.x / n.magnitude);
        }

        public static float ArgumentOfPeriapsis(Vector3 n, Vector3 e) // ω
        {
            float equation = MathF.Acos(Vector3.Dot(n, e) / (n.magnitude * e.magnitude));
            if (e.z < 0)
            {
                return 2 * MathF.PI - equation;
            }
            else
            {
                return equation;
            }
        }

        public static float SemiMajorAxis(float u, float SME) // a
        {
            return -1 * (u / (2 * SME));
        }

        public static float TrueAnomaly(Vector3 e, Vector3 r) // ν
        {
            return Mathf.Acos(Vector3.Dot(e, r) / (e.magnitude * r.magnitude));
        }
    }

    public List<Ship> ships = new List<Ship>();
    public List<GameObject> paths = new List<GameObject>();

    IEnumerator GetElements()
    {
        yield return new WaitForSeconds(1.005f);
        Predictor();
    }

    void Start()
    {
        central = GameObject.Find("M1");
        planetMass = central.AddComponent<Rigidbody>();
        planetMass.mass = mass;
        planetMass.angularDamping = 0;
        size = central.GetComponent<Transform>().localScale.magnitude;
        planetMass.angularVelocity = new Vector3(0, 0.1f, 0);
        Ship mil = new Ship("Player", shipCoords);       
        ships.Add(mil);
        mil.mass.AddForce(shipThrust, ForceMode.Impulse);
        foreach (Ship i in ships)
        {
            i.body.GetComponent<Renderer>().material = shipTexture;
        }
        StartCoroutine(GetElements());
    }

    void Trail()
    {
        foreach(Ship i in ships)
        {
            Vector3 tailPos = i.pos.position;
            GameObject tail = new GameObject("Tail");
            tail.GetComponent<Transform>().position = tailPos;
            tail.GetComponent<Transform>().localScale = i.pos.localScale * 0.5f;
            MeshFilter mFilter = tail.AddComponent<MeshFilter>();
            GameObject go = GameObject.CreatePrimitive(PrimitiveType.Cube);
            mFilter.mesh = go.GetComponent<MeshFilter>().mesh;
            mFilter.gameObject.AddComponent<MeshRenderer>();
            tail.GetComponent<Renderer>().material = shipTexture;
            Destroy(go);
        }
    }

    void Predictor()
    {
        foreach(GameObject o in paths)
        {
            Destroy(o);
        }
        paths.Clear();
        foreach (Ship i in ships)
        {
            for(int x = 0; x < 360; x++)
            {
                float r = (i.semaj * (1 - MathF.Pow(i.ecvec.magnitude, 2))) / (1 + i.ecvec.magnitude * MathF.Cos(x * Mathf.Deg2Rad));
                Vector3 flatPos = new Vector3(r * MathF.Cos(x * Mathf.Deg2Rad), r * MathF.Sin(x * Mathf.Deg2Rad), 0);

                Quaternion rotationPeriapsis = Quaternion.Euler(0, 0, i.argumentOfPeriapsis * Mathf.Rad2Deg);
                Quaternion rotationInclination = Quaternion.Euler(i.inclination * Mathf.Rad2Deg, 0, 0);
                Quaternion rotationAscendingNode = Quaternion.Euler(0, 0, i.rAscensionNode * Mathf.Rad2Deg);

                Quaternion totalRotation = rotationAscendingNode * rotationInclination * rotationPeriapsis * Quaternion.Euler(90, 0, 0);
                Vector3 nodePos = totalRotation * flatPos;

                GameObject orbit = new GameObject("Orbit");
                orbit.GetComponent<Transform>().position = nodePos;
                orbit.GetComponent<Transform>().localScale = i.pos.localScale * 0.5f;
                MeshFilter mFilter = orbit.AddComponent<MeshFilter>();
                GameObject go = GameObject.CreatePrimitive(PrimitiveType.Cube);
                mFilter.mesh = go.GetComponent<MeshFilter>().mesh;
                mFilter.gameObject.AddComponent<MeshRenderer>();
                orbit.GetComponent<Renderer>().material = shipTexture;
                Destroy(go);
                paths.Add(orbit);
            }
        }
    }

    void FixedUpdate()
    {
        foreach(Ship i in ships)
        {
            i.gnd = MathF.Abs(i.pos.position.magnitude) - (size / 2);
            i.spam = Equations.SpecificAngularMomentum(i.pos.position, i.mass.linearVelocity);
            i.novec = Equations.NodeVector(i.spam);
            i.std = Equations.StandardGravParam(gravity, central.GetComponent<Rigidbody>().mass, i.mass.mass);
            i.ecvec = Equations.EccentricityVector(i.pos.position, i.mass.linearVelocity, i.spam, i.std);
            i.smen = Equations.SpecialMechanicalEnergy(i.mass.linearVelocity, i.pos.position, i.std);
            i.semaj = Equations.SemiMajorAxis(i.std, i.smen);
            i.inclination = Equations.Inclination(i.spam);
            i.rAscensionNode = Equations.RightAscendingNode(i.novec);
            i.argumentOfPeriapsis = Equations.ArgumentOfPeriapsis(i.novec, i.ecvec);
            i.trueAnomaly = Equations.TrueAnomaly(i.ecvec, i.pos.position);
            i.mass.AddForce(i.pos.position * -1 * Equations.Gravity(i.pos, planetMass, i.mass, gravity), ForceMode.Force);
            i.pos.eulerAngles = i.mass.linearVelocity;

            if (Input.GetKeyDown("space"))
            {
                i.mass.AddForce(0, -50, 0, ForceMode.Impulse);
            }
        }
    }
}

问题根源与修复方案 #

1. 升交点赤经的象限判断缺失 #

原方法仅计算Arccos结果，未考虑升交点向量n的Z分量符号，导致象限判断错误。同时添加0倾角的特殊处理（此时升交点无意义，返回0避免除以0）：

public static float RightAscendingNode(Vector3 n) // Ω
{
    if (n.magnitude < 1e-6)
        return 0f;
        
    float angle = MathF.Acos(n.x / n.magnitude);
    return n.z < 0 ? 2 * MathF.PI - angle : angle;
}

2. 真近点角的象限判断缺失 #

原方法未考虑径向速度（r·v）的符号，导致轨道方向反向。修改方法并传入速度向量：

public static float TrueAnomaly(Vector3 e, Vector3 r, Vector3 v) // ν
{
    float dotEr = Vector3.Dot(e, r);
    float denominator = e.magnitude * r.magnitude;
    if (denominator < 1e-6)
        return 0f;
        
    float angle = MathF.Acos(dotEr / denominator);
    float radialVelocity = Vector3.Dot(r, v);
    return radialVelocity < 0 ? 2 * MathF.PI - angle : angle;
}

同时修改FixedUpdate中的调用：

i.trueAnomaly = Equations.TrueAnomaly(i.ecvec, i.pos.position, i.mass.linearVelocity);

3. 轨道旋转顺序错误 #

原旋转顺序不符合轨道要素定义，正确顺序为：升交点赤经（Y轴）→ 倾角（X轴）→ 近心点幅角（Z轴），并移除多余的90度旋转：

// 修正旋转顺序与轴
Quaternion rotationAscendingNode = Quaternion.Euler(0, i.rAscensionNode * Mathf.Rad2Deg, 0);
Quaternion rotationInclination = Quaternion.Euler(i.inclination * Mathf.Rad2Deg, 0, 0);
Quaternion rotationPeriapsis = Quaternion.Euler(0, 0, i.argumentOfPeriapsis * Mathf.Rad2Deg);

Quaternion totalRotation = rotationAscendingNode * rotationInclination * rotationPeriapsis;
Vector3 nodePos = totalRotation * flatPos;

4. 引力计算的向量与公式错误 #

原方法返回标量且分母为r³（正确应为r²），修改为向量计算：

public static Vector3 Gravity(Vector3 r, Rigidbody m1Mass, Rigidbody m2Mass, float gravityConstant)
{
    float rMagnitude = r.magnitude;
    if (rMagnitude < 1e-6)
        return Vector3.zero;
        
    float forceMagnitude = (gravityConstant * m1Mass.mass * m2Mass.mass) / MathF.Pow(rMagnitude, 2);
    return -r.normalized * forceMagnitude;
}

修改FixedUpdate中的调用：

i.mass.AddForce(Equations.Gravity(i.pos.position, planetMass, i.mass, gravity), ForceMode.Force);

内容的提问来源于stack exchange，提问作者Privvet