モデルロケット高度計算プログラム

モデルロケット高度計算プログラムを書いてみました。
高度の計算式、

をルンゲクッタ法で計算。 ここで、抗力の係数の式

と、燃焼薬が燃焼して減りながら変化する重量の式

を使う。

c6-5エンジンの推力変化

計算結果

代表面積、抗力係数、気体密度
は結構適当に与えているので検討の余地があります。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int dim = 2;
void get_ev(int dim,double A[][dim],double _Complex EV[],double _Complex V[][dim]);
void func(double A[][dim],double [], double []);
void runge(double A[][dim],double []);
void initialize(double []);
void determin_matrix(double A[][dim]);

const double Tmax = 8.;
const double dt = 0.001;
int imax;
int youso_max = 2;
double t;


double g = 9.8;
double C_D= 1;
double S  = 0.0019;
double rho= 1.225;
double m0 = 100;


//エンジン固有
double mp = 12.48;
double Imp = 10.;
double ave_F = 6.;
double F;


double k;
double w;
double T_p;


int main(void){
	imax = Tmax/dt;
	int i,j,n;
	
	FILE *fp;	
	fp = fopen("rocket.csv","w");

	double x[youso_max];
	double A[dim][dim];
	
	initialize(x);
	for(i=0;i<=imax;i++){
		
		fprintf(fp,"%e,",t);
		for(j=0;j<youso_max;j++){
			fprintf(fp,"%e,",x[j]);
		}
		fprintf(fp,"%e,",F);
		fprintf(fp,"\n");
		
		determin_matrix(A);
		runge(A,x);
		t = t+dt;		
	}
	
}

void determin_matrix(double A[][dim]){
	int i,j;
	double t1 = 0.2,t2=0.4,t3=1.8,t4=1.9;
	
	T_p = Imp / ave_F;
		
	if(t<=t1){
		F = 14.5*t/t1;
	}
	else if((t1<t)&&(t<=t2)){
		F = 14.5-9.*(t-t1)/(t2-t1);
	}
	else if((t2<t)&&(t<=t3)){
		F = 4.;
	}
	else if((t3<t)&&(t<=t4)){
		F = 4.*(1-(t-t3)/(t4-t3));
	}
	else{
		F = 0;
	}
	
	
	if(t <= T_p){
		w = (m0/1000. - (mp/1000. * t / T_p))*g;
	}
	else{
		w = (m0/1000. - mp/1000.)*g;	
	}
	k = C_D * S * 0.5 * rho;
	
	A[0][0]=0.; A[0][1]=1.; 
	A[1][0]=0.; A[1][1]=-g*k/w; 


}

void func(double A[][dim],double x[],double f[]){
	double C[2];
	double eta=0;	
	
	if(x[1] >= 0){
		eta = 1.0;
	}
	else{
		eta = -1.0;
	}

	determin_matrix(A);
	C[0] = 0.;
	C[1] = g*(F/w - 1.);
	

	f[0] = A[0][0]*x[0] +     A[0][1]     *x[1] + C[0];
	f[1] = A[1][0]*x[0] + eta*A[1][1]*x[1]*x[1] + C[1];
	return;
}

void initialize(double x[]){
	int i;
	t=0.; 	
	x[0]=0.;
	x[1]=0.;
}

void runge(double A[][dim],double x[]){

	double k1[youso_max],k2[youso_max],k3[youso_max],k4[youso_max];
	double f1[youso_max],f2[youso_max],f3[youso_max],f4[youso_max];
	double temp1[youso_max],temp2[youso_max],temp3[youso_max],temp4[youso_max];
	int i;
	//一つ目
	for(i=0;i<youso_max;i++){
		temp1[i] = x[i];
	}
	func(A,temp1,f1);
	for(i=0;i<youso_max;i++){
		k1[i] = f1[i]*dt;
	}

	//二つ目
	for(i=0;i<youso_max;i++){
		temp2[i] = x[i] + 0.5*k1[i];
	}
	func(A,temp2,f2);
	for(i=0;i<youso_max;i++){
		k2[i] = f2[i]*dt;
	}

	//三つ目
	for(i=0;i<youso_max;i++){
		temp3[i] = x[i] + 0.5*k2[i];
	}
	func(A,temp3,f3);
	for(i=0;i<youso_max;i++){
		k3[i] = f3[i]*dt;
	}

	//四つ目
	for(i=0;i<youso_max;i++){
		temp4[i] = x[i] + k3[i];
	}
	func(A,temp4,f4);
	for(i=0;i<youso_max;i++){
		k4[i] = f4[i]*dt;
	}
	
	for(i=0;i<youso_max;i++){
		x[i] = x[i] + (k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.;
	}
	 	

	return;
}