编程高手进来帮我找下这程序问题出在哪……

文号t16766

8755

阿童木升级版2009/07/03软件综合 IP:重庆

很久很久以前，我发表了一篇文章《浅析模型火箭轨道》提出利用计算机解决“一个受与其速度成一非2k/π的角度的加速度的物体的运动分析”（https://www.kechuang.org/t/16723）

由于本人编程能力基本没有（顶多看得懂程序）故请朋友帮忙，现在大概是写出来了，但是有严重BUG，根本无法正常计算。希望各位大虾帮忙……谢谢了

BUG主要见下图：

源码如下：
#include<iostream>
#include<math.h>
#include<iomanip>
using namespace std;
double Ft,M0,Mt,aj,a,g,p,A,D,t,Tz,i,Fj;
double Tj,x,y,v0x,v0y,M,fjf,ax,ay,vmx,vmy,xz,yz,fx,fy;
double sin(double shuju);
double cos(double shuju);
double atan(double shuju);

void init()
{
     cout<<"输入火箭推力Ft，单位牛顿（N）"<<endl;
     cin>>Ft;
     if(Ft<=0){cout<<"输入数据为非正常数据，请重新输入各个已知量。";return init();}
     cout<<"输入火箭初始质量M0，单位千克（KG）"<<endl;
     cin>>M0;
     if(M0<=0){cout<<"输入数据为非正常数据，请重新输入各个已知量。";return init();}
     cout<<"输入火箭终了质量Mt，单位千克（KG）"<<endl;
     cin>>Mt;
     if(Mt<=0||Mt>=M0) {cout<<"输入数据为非正常数据，请重新输入各个已知量。";return init();}
     cout<<"输入火箭发射倾角a，单位度（°）"<<endl;
     cin>>aj;
     if(aj<=0||aj>=180){cout<<"输入数据为非正常数据，请重新输入各个已知量。";return init();}
     a=(aj*3.1415926535897932384626)/180;
     g=9.8056;
     p=1.29;
     cout<<"输入气阻有效面积A，单位平方米（m2）"<<endl;
     cin>>A;
     if(A<=0) {cout<<"输入数据为非正常数据，请重新输入各个已知量。";return init();}
     cout<<"输入气阻常数D"<<endl;
     cin>>D;
     if(D<=0){cout<<"输入数据为非正常数据，请重新输入各个已知量。";return init();}
     cout<<"输入火箭推力有效时间t,单位秒（S）"<<endl;
     cin>>t;
     if(t<=0){cout<<"输入数据为非正常数据，请重新输入各个已知量。";return init();}
     cout<<"输入时间间隔Tz,必须小于火箭推力有效时间，单位秒（s）"<<endl;
     cin>>Tz;
     if(Tz<=0||Tz>t){cout<<"输入数据为非正常数据，请重新输入各个已知量。";return init();}
     x=0;
     y=0;
     Tj=0;
     v0x=0;
     v0y=0;
}

void jisuan1()
{
     M=((Mt-M0)/t)*Tj+M0;
     Fj=(-(M*M*M*M)+sqrt((M*M*M*M)+2*p*A*D*t*t*Ft*M*M))/p*A*D*Tj*Tj;
     ax=(fjf*cos(a))/M;
     ay=(fjf*sin(a)-M*g)/M;
     vmx=ax*Tz;
     vmy=ay*Tz;
     xz=v0x*t+((1/2)*ax*t*t);
     yz=v0y*t+((1/2)*ay*t*t);
     x=x+xz;
     y=y+yz;
     if(y<0){cout<<"Y坐标小于0，计算结束。"<<endl;return;}
     Tj=Tj+Tz;
     cout<<"输出时刻："<<setprecision(10)<<Tj<<endl;
     cout<<"输出横坐标："<<setprecision(10)<<x<<endl;
     cout<<"输出纵坐标："<<setprecision(10)<<y<<endl;
     v0x=vmx;
     v0y=vmy;
     a=atan(v0y/v0x);
     fx=(1/2)*p*A*D*v0x*v0x*cos(a);
     fy=(1/2)*p*A*D*v0y*v0y*cos(a);
     ax=fx/Mt;
     ay=(fy+M*g)/Mt;
     xz=v0x*(Tz-t)+(1/2)*ax*(Tz-t)*(Tz-t);
     yz=v0y*(Tz-t)+(1/2)*ay*(Tz-t)*(Tz-t);
     x=x+xz;
     y=y+yz;
     Tj=Tj+Tz;
     if(y<0){cout<<"Y坐标小于0，计算结束。"<<endl;return;}
     cout<<"输出时刻："<<setprecision(10)<<Tj<<endl;
     cout<<"输出横坐标："<<setprecision(10)<<x<<endl;
     cout<<"输出纵坐标："<<setprecision(10)<<y<<endl;


}

int main()
{
    init();
    for(i=0;i<=t;i+=Tz)
    {
       jisuan1();
    }
    system("pause");
    return 0;
}
另附全部组件（程序，代码，论文）

鄙人QQ：563184214
邮箱：BeyoundATOM@XXXXXXXXX