мат. модель летательного аппарата

    public class PlaneSim
    {
        public double _cx0;
        public double _cy0;
        public double _cz0;
        public double _k1;
        public double _k2;
        public double _k3;
        public double _sq;   // площадь крыла
        public double _den;  // плотность воздуха
        public double _m;
        public double _g;
        public double _rz;
        public double _phi;

        public double _P;    // сила тяги
        public double _a;    // угол атаки
        public double _b;    // угол скольжения
        public double _y;    // угол крена

        public double _V;    // скорость
        public double _Q;
        public double _Cr;   // курс
        public double _H;    // высота
        public double _Lt;   // широта
        public double _Ln;   // долгота

        public double _px;
        public double _pz;

        public PlaneSim()
        {
            _cx0 = 0.034;
            _cy0 = 0.42;
            _cz0 = 0.0;
            _k1 = 0.038;
            _k2 = 0.09;
            _k3 = 0.01;
            _sq = 1.775;
            _den = 1.225;
            _m = 63.3;
            _g = 9.81;
            _rz = 6.378 * 1000000.0;
            _phi = 0.0;

            _P = 200.0;
            _a = 0.0;
            _b = 0.0;
            _y = 0.0;

            _V = 30.0;
            _Q = 0.2;
            _Cr = 0.0;
            _H = 200.0;
            _Lt = 0.0;
            _Ln = 0.0;
        }

        public void Integrate(double dt)
        {
            double sinQ, cosQ;
            double sinY, cosY;
            double tanLt, cosCr;
            double sinCr;
            double secLt;

            double f1, f2, f3, f4, f5, f6, f7, f8;
            double dV, dQ, dCr, dH, dLt, dLn;

            sinQ = Math.Sin(_Q);
            cosQ = Math.Cos(_Q);

            sinY = Math.Sin(_y);
            cosY = Math.Cos(_y);

            tanLt = Math.Tan(_Lt);
            cosCr = Math.Cos(_Cr);
            sinCr = Math.Sin(_Cr);

            secLt = 1.0 / Math.Cos(_Lt);

            // dV
            f1 = _cy0 + _k2 * _a;
            f2 = _cx0 + _k1 * f1 * f1 - (_cz0 + _k3 * _b) * _b;
            f3 = _P - f2 * (_den * _V * _V * _sq * 0.5);
            f4 = (f3 / _m) - _g * sinQ;
            dV = f4;

            // dQ
            f1 = _P * ((_a * Utils.DEG2RAD - _phi) * cosY + _b * sinY);
            f2 = _cy0 + _k2 * _a;
            f3 = (_cx0 + _k1 * f2 * f2) * _b * sinY;
            f4 = (_cy0 + _k2 * _a) * cosY;
            f5 = (_cz0 + _k3 * _b) * sinY;
            f6 = f1 - (f3 - f4 + f5) * (_den * _V * _V * _sq * 0.5);
            f7 = (f6 / (_m * _V)) - ((_g * cosQ) / _V) + ((_V * cosQ) / (_rz + _H));
            dQ = f7;

            // dCr
            f1 = _P * ((_a - _phi) * sinY - _b * cosY);
            f2 = _cx0 + _k1 * (_cy0 + _k2 * _a);
            f3 = f2 * f2 * _b * cosY;
            f4 = (_cy0 + _k2 * _a) * sinY;
            f5 = (_cz0 + _k3 * _b) * cosY;
            f6 = f1 + (f3 + f4 + f5) * (_den * _V * _V * _sq * 0.5);
            f7 = f6 / (_m * _V * cosQ);
            f8 = _V * (cosQ / (_rz + _H)) * tanLt * cosCr;
            dCr = -f7 - f8;

            // dH
            dH = _V * sinQ;

            // dLt
            dLt = (_V / (_rz + _H)) * cosQ * sinCr;

            // dLn
            dLn = (_V / (_rz + _H)) * secLt * cosQ * cosCr;

            // Euler
            _V += dV * dt;
            _Q += dQ * dt;
            _Cr += dCr * dt;
            _H += dH * dt;
            _Lt += dLt * dt;
            _Ln += dLn * dt;

            // Position
            _px = (_rz + _H) * _Ln;
            _pz = (_rz + _H) * _Lt;
        }
    }

            mPlaneSim.Integrate(mDeltaTime);

            float alpha = (float)(mPlaneSim._Q + mPlaneSim._a) * Utils.DEG2RAD;
            float cr = (float)mPlaneSim._Cr;
            float gamma = (float)mPlaneSim._y;

            mPlaneDir.X = -cr;
            mPlaneDir.Y = -alpha;
            mPlaneDir.Z = -gamma;

            float scale = 10.0f;
            mPlanePos.X = -(float)mPlaneSim._pz * scale;
            mPlanePos.Y = (float)mPlaneSim._H;
            mPlanePos.Z = (float)mPlaneSim._px * scale;

            PlaneMesh.Position = mPlanePos;
            PlaneMesh.Rotation = mPlaneDir;
            PlaneMesh.Draw();