プログラム

#include <iostream>
#include<math.h>
#include<stdio.h>
#define     PI      3.141592653589793238463
using namespace std;
//主翼揚力係数　NACA2312のデータを直線で近似
double CL1(double theta){
	double t=theta/2.0/PI*360;
	if (t<22) return(0.3+0.0708*t);
	else return(1.65-0.0588*(t-22));
}
//主翼抗力係数　　NACA2312のデータを放物線で近似
double CD1(double theta){
	double t=theta/2.0/PI*360;
	return(0.05+1.736e-3*(t+4)*(t+4));
}
//水平尾翼揚力係数
double CL2(double theta){
	double t=theta/2.0/PI*360;
	if (t<22) return(0.3+0.0708*t-0.3);
	else return(1.65-0.0588*(t-22)-0.3);
}
//水平尾翼抗力係数
double CD2(double theta){
	double t=theta/2.0/PI*360;
	return(0.05+1.736e-3*(t+4)*(t+4));
}

void main(void){
	const double mass=6.5e-3+0.6e-3+3.0e-3; //kg　機体重量　主翼＋尾翼＋胴体
	const double moment=6.5e-3*0.025*0.025+0.6e-3*0.18*0.18+1.5e-3/3*(0.2*0.2)*2; //kg・m^2　重心周りの回転モーメント
	const double X1=2.5e-2;　//重心から主翼揚力中心までの距離
	const double X2=18.0e-2;　//重心から尾翼揚力中心までの距離
	const double G=9.8; //重力加速度
	const double S1=0.04*0.30; //m^2　主翼面積
	const double S2=0.04*0.08; //m^2　水平尾翼面積
	const double dT=10e-6; //sec　計算の時間刻み
	const double rho_=1.293/2 ; //kg/m^3　空気の密度

	double vx=3.0,vz=0.0; // m/s 初速　x:水平 z:垂直　
	double x=0.0,z=2.0; // m 初期位置　
	double theta_a=0.0/360*2*PI; // rad　機体の水平面に対する角度
	double omega_a=0.0; //機体の回転の角速度
	double time; //時刻
	for(time=0; time <= 30; time+=dT){
		double v2=vx*vx+vz*vz;
		double theta_v=atan2(vz,vx); //機体速度ベクトルから機体が受ける相対風の角度を算出
		double theta_in=theta_a-theta_v; //迎え角の計算　（主翼角度に対する相対風の角度）
		double theta_1=atan2(CL1(theta_in),CD1(theta_in)); //進行方向に対する主翼の受ける揚力と効力による合力の方向
		double theta_1_=theta_1-theta_v;　//水平面に対する主翼の受ける揚力と効力による合力の方向
		double theta_2=atan2(CL2(theta_in),CD2(theta_in));　//進行方向に対する水平尾翼の受ける揚力と効力による合力の方向
		double theta_2_=theta_2-theta_v;　//水平面に対する水平尾翼の受ける揚力と効力による合力の方向
	//cout <<theta_v<<" "<<theta_1<<" "<<theta_2<<"\n";
		double ax=-rho_*v2*S1*(CL1(theta_in)+CD1(theta_in))*cos(theta_1_)/mass
				  -rho_*v2*S2*(CL2(theta_in)+CD1(theta_in))*cos(theta_2_)/mass;　//水平方向加速度
		double az=+rho_*v2*S1*(CL1(theta_in)+CD1(theta_in))*sin(theta_1_)/mass
				  +rho_*v2*S2*(CL2(theta_in)+CD1(theta_in))*sin(theta_2_)/mass 
				  -G; //垂直方向加速度
		double temp;
		if(omega_a>=0) temp=1;
		else temp=-1;

		double alpha_a=+X1*rho_*v2*S1*(CL1(theta_in)+CD1(theta_in))*sin(theta_1_+theta_a)/moment
					 -X2*rho_*v2*S2*(CL2(theta_in)+CD2(theta_in))*sin(theta_2_+theta_a)/moment
					 -temp* rho_*((omega_a*X1)*(omega_a*X1)*S1+(omega_a*X2)*(omega_a*X2)*S2)/moment;//重心周りの加角速度
		x+=vx*dT+0.5*ax*dT*dT; //x位置の変化
		vx+=ax*dT;　//x速度の変化
		z+=vz*dT+0.5*az*dT*dT;//z位置の変化
		vz+=az*dT;//z速度の変化
		theta_a+=omega_a*dT+0.5*alpha_a*dT*dT;//機体角度の変化
		omega_a+=alpha_a*dT;//機体角速度の変化
		if(z<0) break;//地面(z=0)に着いたら終了
		if( int(time/dT/1000)*1000==int(time/dT))		cout << x<<" "<<z<< " "<<cos(theta_a)<<" "<<sin(theta_a)<<"\n";//10msごとに機体の位置と角度を書き出す
	}
}