其中空气阻力系数k1应根据速度变化，还未找到解决方案

#include<iostream>
#include<iomanip>
using namespace std;
//#define PI 3.141592
//#define g 9.788

constexpr double PI = 3.141592;
constexpr double g = 9.7833;
constexpr double e = 2.718281;
/**
   * @brief 计算pitch
   * @param a							当前pitch角度
   * @param _bullet_speed				子弹速度
   * @param const cv::Mat& _pos_in_ptz	待击打目标位置信息
   * @return cv::Point3f				返回 Yaw Pitch 轴的偏移量和深度（m）
   * @author XX
   */
double solve(double v0, double a, double arr[2]) {
	a = a * PI / 180;
	//浮点数不使用科学计数法
	cout.setf(ios::fixed, ios::floatfield);
	//保留6位小数
	cout.precision(6);

	// k1 = k / m , k由空气密度rho（温度T，大气压强P）、迎风面积S（直径d、）、雷诺数（10m/s~30m/s综合速度v）、空气动态粘度系数mu所得
	// rho（20℃、101.3k） = 1.204 kg / m³ 
	// S（0.042m） = 0.2125 * 0.2125 * PI
	// mu = 1.825 e-5 = 0.00001825
	// 阻力系数C≈0.5
	//42.5mm弹丸	0.009342
	//16.8mm弹丸	0.019616	
	double k1 = 0.009342;
	//收敛速度系数（不宜太大，会爆炸）
	//double alpha = 0.1;
	double k = tan(a);
	double tempPoint = k * arr[0];
	double targetPoint = arr[1];
	if (arr[1] >= 0) {
		for (int i = 0; i < 20; i++) {
			double vx0 = v0 * cos(a);
			double vy0 = v0 * sin(a);
			double t = (pow(e, k1 * arr[0]) - 1) / (k1 * vx0);
			double vx = vx0 / (k1 * vx0 * t + 1);
			double h_max = vy0 * vy0 / (2 * g);
			double t_h = vy0 / g;
			double delta_t = t - t_h;
			double fall_h = 0.5 * g * delta_t * delta_t;
			double realPoint = h_max - fall_h;
			double deltaH = targetPoint - realPoint;
			tempPoint = tempPoint + deltaH;
			//k1 = k1 - alpha * deltaH;
			a = atan(tempPoint / arr[0]);
			std::cout << "第" << setw(2) << i + 1 << "次迭代: 仰角：" << a * 180 / PI << ", 临时目标点y值 : " << tempPoint
				<< ", 高度误差 : " << deltaH << std::endl;
		}
		return a;
	}
}

int main() {
	double arr[2] = { 20,1 };
	double v0 = 16;
	double a = 0;
	double b = solve(v0, a, arr);
}