导弹追踪问题

相关资料见
Matlab建模—导弹追踪问题
追及问题真的有那么简单吗？–追及问题的本质
geogebra追及问题
以B站视频那道题的数据为例。

simucpp目前只支持绘制波形，不支持绘制坐标图，因此打印出导弹的横纵坐标后将打印结果复制到单独的data.csv文件中，并用python画图。绘图代码如下：

import numpy as np
import csv
import zhnplot
filename = 'data.csv'
with open(filename) as f:
    reader = csv.reader(f)
    ans = list(reader)
x, y = [], []
for adata in ans:
    x.append(adata[0])
    y.append(adata[1])
x=np.array(x, dtype=np.float32)
y=np.array(y, dtype=np.float32)
zhnplot.plot(x, y)

其中zhnplot见 封装matplotlib绘制单个子图。

方法一

$$\begin{gathered} \frac{dx}{dt}=v_2\cos \theta \\ \frac{dy}{dt}=v_2\sin \theta \\ \tan \theta =\frac{v_{1}t-y}{-x} \end{gathered}$$

方法二

$$\begin{gathered} \frac{dy}{dx}=\frac{v_{1}t-y}{-x} \\ {\int }_{-L}^x\sqrt{1+y^{\prime 2}}dx=v_{2}t \end{gathered}$$
上式代入下式得
$${\int }_{-L}^x\sqrt{1+y^{\prime 2}}dx=(y-x{y}^{\prime })\frac{v_2}{v_1}$$
两边同时求导得
$$\sqrt{1+y^{\prime 2}}=-\frac{v_2}{v_1}x{y}^{\prime \prime }$$
令$z=\frac{dy}{dx}$，则
$$\frac{dz}{dx}=-\frac{v_1}{v_{2}x}\sqrt{1+z^2}$$

方法三

也是来自B站视频的方法，可以算出时间，但轨迹不知道咋算，而且也不理解为什么可以这样算时间。
$$\begin{gathered} v_{2}t-v_{1}t\sin \theta =L \\ {v}_{1}t-v_{2}t\sin \theta =0 \end{gathered}$$
解得
$$t=\frac{L}{v_2-v_1\sin \theta }=\frac{L}{v_2-\frac{v_1^2}{v_2}}$$
simucpp代码：

// 导弹追踪问题
#include "simucpp.hpp"
#include <cmath>
int main()
{
    using namespace simucpp;
    using namespace std;
    Simulator sim1, sim2;

    constexpr double v1 = 1;
    constexpr double v2 = 2;
    constexpr double L = 5;

    // 方法一
    MInput *s1shipdist = new MInput(&sim1);  //输出敌舰向y轴正方向移动的距离
    MFcn *s1fcnx = new MFcn(&sim1);  // 输出导弹的速度的x分量
    MFcn *s1fcny = new MFcn(&sim1);  // 输出导弹的速度的y分量
    MFcn *s1divider = new MFcn(&sim1);  // 输出导弹的x坐标的倒数
    MFcn *s1theta = new MFcn(&sim1);  // 输出导弹的角度
    MIntegrator *s1missilex = new MIntegrator(&sim1);  // 输出导弹的x坐标
    MIntegrator *s1missiley = new MIntegrator(&sim1);  // 输出导弹的y坐标
    MSum *sum1 = new MSum(&sim1);  // 输出敌舰y坐标减去导弹的y坐标
    MProduct *s1tan = new MProduct(&sim1);  // 输出导弹的角度的正切值

    s1fcnx->Set_Function([](double u){return v2*cos(u);});
    s1fcny->Set_Function([](double u){return v2*sin(u);});
    s1divider->Set_Function([](double u){return -1/u;});
    s1theta->Set_Function([](double u){return atan(u);});
    s1shipdist->Set_InputFunction([](double t){return v1*t;});
    s1missilex->Set_InitialValue(-L);
    s1missiley->Set_InitialValue(0);

    sim1.connect(s1theta, s1fcnx);
    sim1.connect(s1theta, s1fcny);
    sim1.connect(s1fcnx, s1missilex);
    sim1.connect(s1fcny, s1missiley);
    sim1.connect(s1shipdist, sum1);
    sim1.connect(s1missiley, sum1);
    sum1->Set_InputGain(-1);
    sim1.connect(s1missilex, s1divider);
    sim1.connect(sum1, s1tan);
    sim1.connect(s1divider, s1tan);
    sim1.connect(s1tan, s1theta);

    sim1.Initialize();
    sim1.Set_ShowWave(false);
    double t, missilex, missiley;
    cout.precision(12);
    do {
        sim1.Simulate_OneStep();
        t = sim1.Get_t();
        missilex = s1missilex->Get_OutValue();
        missiley = s1missiley->Get_OutValue();
        // cout << missilex << ", "<< missiley << endl;
    } while (missilex < 0);
    // cout << "t: " << t << endl;

    // 方法二
    MIntegrator *s2missiley = new MIntegrator(&sim2);  // 输出导弹的y坐标
    MIntegrator *s2z = new MIntegrator(&sim2);  // 输出导弹的y分量速度
    MFcn *s2fcn1 = new MFcn(&sim2);  // 输出√(1+z^2)
    MInput *s2in = new MInput(&sim2);  // 输出(v1)(v2x)
    MProduct *s2zd = new MProduct(&sim2);  // 输出导弹的y分量加速度

    sim2.Set_t(-L);  // 方法一的自变量是t，方法二的自变量是x
    s2fcn1->Set_Function([](double u){return sqrt(1+u*u);});
    s2in->Set_InputFunction([](double t){return -v1/v2/t;});

    sim2.connect(s2in, s2zd);
    sim2.connect(s2z, s2fcn1);
    sim2.connect(s2fcn1, s2zd);
    sim2.connect(s2zd, s2z);
    sim2.connect(s2z, s2missiley);

    sim2.Initialize();
    sim2.Set_ShowWave(false);
    double missileyd;
    do {
        sim2.Simulate_OneStep();
        t = sim2.Get_t();
        missileyd = s2z->Get_OutValue();
        missiley = s2missiley->Get_OutValue();
        // cout << t << ", "<< missiley << endl;
    } while (t < 0);

    return 0;
}