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

モデルロケット高度計算プログラムを書いてみました。
高度の計算式、
\frac{d^2x}{dt^2} = -\frac{gk}{w}\left ( \frac{dx}{dt}\right )^2 +g (\frac{F}{w} -1)
をルンゲクッタ法で計算。 ここで、抗力の係数の式
k = \frac{\rho C_D S}{2}
と、燃焼薬が燃焼して減りながら変化する重量の式
w = m(t)g
を使う。

c6-5エンジンの推力変化
https://estesrockets.com/wp-content/uploads/2013/11/C6_thrustcurve_1_1.jpg

計算結果
f:id:takao_science:20201008233349p:plain

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

#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;
}

はてなブログに実装されたjavascriptを埋め込む方法【人生初のjavascriptのプログラミング】

人生で初めてjavascriptを触りました。
このサイトを参考にしました。

techacademy.jp

↓サンプルコートをそのまま書いた電卓です。

 

 

 

便利ですね。

javascriptがはてなブログ上でうまく動作しなかったので、JSFiddleというサイトから埋め込みコードを使用して、埋め込みました。

↓こちらのサイトを参考にしました。

cpoint-lab.co.jp


これで計算がブラウザ上で出来るようになるぞー！

 

【C言語】LAPACKを使った多変数の最小二乗法の計算方法【Cygwin】

Cygwinを用いて、多変数の最小二乗法を行ったので、備忘録として、下記に留める。

　I = Σ_{ij} a_{ij} T^i \omega^j

となる a_{ij}を求めるプログラム。（今回は5次まで）
I\ T\ \omegaの順でスペースで区切られた近似したいデータファイル(test.dat)を読み込んで計算します。
gcc sample.c -llapack -lblas -lm のコマンドでコンパイルします。

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

int dim,dim1,dim2,dim3;
void matrix_transposition(int,int,double  **A,double  **t_A);//行列転置
void matrix_initialize(int,int,double  **A); //行列要素すべて０
void matrix_inverse(int,double **A,double **B);
void matrix_inverse2(int,double **A,double **B);
void matrix_calculate(int,int,int, double **A,double **B,double **C);
double g_calc(int,double T,double omega);
void test_printf(int,int,double **A,char *);//テスト用表示

#define SIZE 21
#define NSIZE 30

int main(void){
  int i,j;
  FILE *fp;
  int ret;
  int cnt = 0;
  
 
  char fname[100] = "test.dat";
	//printf("filename=");
   //scanf("%s",&fname);
	
	
  fp = fopen( fname , "r" );
  if( fp == NULL ){
   printf( "%sファイルが開けません\n", fname );
    return -1;
  }

  while((ret = getc(fp)) != EOF) {
		 if(ret == '\n') cnt++;
  }
 
	
  //int NSIZE = cnt;
	printf("NSIZE = %d\n",NSIZE);
 static double S_total,S[NSIZE],SI[NSIZE],sigma,rS[SIZE];
	
	double **G;
	G = malloc(sizeof(double *) * SIZE);
	for (i=0;i<SIZE;i++) {
		G[i] = malloc(sizeof(double) * NSIZE);
	}
    double **t_G;
	t_G = malloc(sizeof(double *) * NSIZE);
	for (i=0;i<NSIZE;i++) {
		t_G[i] = malloc(sizeof(double) * SIZE);
	}
	double **temp;
	temp = malloc(sizeof(double *) * SIZE);
	for (i=0;i<SIZE;i++) {
		temp[i] = malloc(sizeof(double) * SIZE);
	}
	double **i_temp;
	i_temp = malloc(sizeof(double *) * SIZE);
	for (i=0;i<SIZE;i++) {
		i_temp[i] = malloc(sizeof(double) * SIZE);
	}
	double **temp2;
	temp2 = malloc(sizeof(double *) * SIZE);
	for (i=0;i<SIZE;i++) {
		temp2[i] = malloc(sizeof(double) * 1);
	}
	double **I;
	I = malloc(sizeof(double *) * NSIZE);
	for (i=0;i<NSIZE;i++) {
		I[i] = malloc(sizeof(double) * 1);
	}
	double *T;
	T = malloc(sizeof(double ) * NSIZE);
	double *omega;
	omega = malloc(sizeof(double ) * NSIZE);
	
	double **a;
	a = malloc(sizeof(double *) * SIZE);
	for (i=0;i<SIZE;i++) {
		a[i] = malloc(sizeof(double) * 1);
	}
	double **unit;
	unit = malloc(sizeof(double *) * SIZE);
	for (i=0;i<SIZE;i++) {
		unit[i] = malloc(sizeof(double) * SIZE);
	}
  rewind( fp );
	
  for(i=0;i<NSIZE;i++){
    ret = fscanf( fp, "%lf %lf %lf", &I[i][0],&T[i],&omega[i]) ;
  	//T[i] = T[i];
  	//omega[i] = omega[i];
  	//printf("%e %e %e\n",I[i][0],T[i],omega[i]);
  }
  fclose(fp);
	
	for(i=0;i<SIZE;i++){
		for(j=0;j<NSIZE;j++){
			G[i][j] = g_calc(i,T[j],omega[j]);
			//printf("G[%d][%d] =%e\n",i,j,G[i][j]);
		}
	}
	//test_printf(SIZE,NSIZE,G,"G");
	matrix_transposition(SIZE,NSIZE,G,t_G);
	//test_printf(SIZE,NSIZE,t_G,"t_G");
	matrix_calculate(SIZE,NSIZE,SIZE,G,t_G,temp);
	test_printf(SIZE,SIZE,temp,"temp");
	matrix_inverse(SIZE,temp,i_temp);
	//matrix_calculate(SIZE,SIZE,SIZE,i_temp,temp,unit);
	//test_printf(SIZE,SIZE,unit,"unit");
	//puts("pass3");
	matrix_calculate(SIZE,NSIZE,1,G,I,temp2);
	//test_printf(SIZE,SIZE,i_temp,"i_temp");
	//test_printf(SIZE,SIZE,i_temp,"i_temp");
	matrix_calculate(SIZE,SIZE,1,i_temp,temp2,a);
	
	for(i=0;i<SIZE;i++){
		printf("a[%d] = %e\n",i,a[i][0]);
	}
	
	
	for(j=0;j<NSIZE;j++){
		SI[j]=0.;
		for(i=0;i<SIZE;i++){
			SI[j] += a[i][0]* g_calc(i,T[j],omega[j]);
		}
		S[j] = (I[j][0] - SI[j])*(I[j][0] - SI[j])/((double) NSIZE);
		rS[j] = sqrt(S[j]*((double) NSIZE));
		printf("%e\n",rS[j]);
	}
	
	S_total = 0.;
	for(j=0;j<NSIZE;j++){
		S_total += S[j]; 
	}
	
	
	sigma = sqrt(S_total);
	printf("standard deviation = %e",sigma);
	
	
	free( G );	
	free( t_G );	
	free( temp );	
	free( i_temp );	
	free( temp2 );	
	free( I );	
	free( T );	
	free( omega );
	free( a );	
}

double g_calc(int i,double T,double omega){
	double ans;
	switch (i) {
	case 0:
		ans =1.0;
		return ans;
		break;

	case 1:
		ans = T;
		return ans;
		break;
		
	case 2:
		ans = omega;
		return ans;
		break;
		
	case 3:
		ans =T*T;
		return ans;
		break;
		
	case 4:
		ans =T*omega;
		return ans;
		break;
		
	case 5:
		ans =omega*omega;
		return ans;
		break;
		
	case 6:
		ans = T*T*T;
		return ans;
		break;

	case 7:
		ans = T*T*omega;
		return ans;
		break;
		
	case 8:
		ans = T*omega*omega;
		return ans;
		break;
		
	case 9:
		ans =omega*omega*omega;
		return ans;
		break;
		
		
		
	case 10:
		ans =T*T*T*T;
		return ans;
		break;
		
	case 11:
		ans = T*T*T*omega;
		return ans;
		break;

	case 12:
		ans = T*T*omega*omega;
		return ans;
		break;
		
	case 13:
		ans = T*omega*omega*omega;
		return ans;
		break;
		
	case 14:
		ans =omega*omega*omega*omega;
		return ans;
		break;
		
		
	case 15:
		ans =T*T*T*T*T;
		return ans;
		break;
		
	case 16:
		ans = T*T*T*T*omega;
		return ans;
		break;

	case 17:
		ans = T*T*T*omega*omega;
		return ans;
		break;
		
	case 18:
		ans = T*T*omega*omega*omega;
		return ans;
		break;
		
	case 19:
		ans =T*omega*omega*omega*omega;
		return ans;
		break;
		
	case 20:
		ans =omega*omega*omega*omega*omega;
		return ans;
		break;	
		
	default:
		break;
	}
}


void matrix_calculate(int dim1,int dim2,int dim3, double **A,double **B,double **C){
	int i,j,k;
	matrix_initialize(dim1,dim3,C);
	for(i=0;i<dim1;i++){
		for(j=0;j<dim3;j++){
	  		for(k=0;k<dim2;k++){
				C[i][j]=C[i][j]+A[i][k]*B[k][j];
	  		}
	 	}
	}
}


void matrix_inverse2(int dim,double **A,double **B){
	
	double buf; //一時的なデータを蓄える
	int i,j,k; //カウンタ
	int n=dim;  //配列の次数
	double inv_a[dim][dim];
 	double C[dim][dim];
	
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
 			C[i][j]= A[i][j];
 		}
	}
	
	//単位行列を作る
	for(i=0;i<n;i++){
 		for(j=0;j<n;j++){
 			if(i==j){ B[i][j]=1.0;}
 			else{B[i][j] = 0.0;}
 		}
	}
	//掃き出し法
	for(i=0;i<n;i++){
		buf=1/C[i][i];
		for(j=0;j<n;j++){
			C[i][j]*=buf;
 			B[i][j]*=buf;
		}
		for(j=0;j<n;j++){
 			if(i!=j){
			buf=C[j][i];
			for(k=0;k<n;k++){
				C[j][k]-=C[i][k]*buf;
				B[j][k]-=B[i][k]*buf;
  			}
			}
		}
	}
}

void matrix_inverse(int dim,double **A,double **B){
  	int i,j;
  	long   m       = dim ; // 行のサイズ
  	long   n       = dim ; // 列のサイズ
  	long   lda     = dim ; // mと同じ値
	double *C;
	C = malloc(sizeof(double) * dim);
	
 	for(i=0;i<dim;i++){
	 for(j=0;j<dim;j++){
		C[dim*i+j] = A[i][j];
	 }
	}
	//test_printf(dim,dim,B,"B");
  	long   info           ; 
  	long   ipiv[dim]     ;
  	long   lwork   = dim ;
  	double work[dim]     ;

  	dgetrf_( &m, &n, C, &lda, ipiv, &info);
	if(info !=0){puts("error_LU");}
  	dgetri_( &n, C, &lda, ipiv, work, &lwork, &info);
	if(info !=0){puts("error_inverse");}

	for(i=0;i<dim;i++){
		for(j=0;j<dim;j++){
			B[i][j] = C[dim*i+j];
		}
	}
}

//行列要素すべて0
void matrix_initialize(int dim1,int dim2,double **A){
	int i,j;
	for(i=0;i<dim1;i++){
		for(j=0;j<dim2;j++){
			A[i][j]=0.0;
		}
	}
}

//行列転置
void matrix_transposition(int dim1,int dim2,double **A,double **t_A){
	int i,j;
	for(i=0;i<dim1;i++){
		for(j=0;j<dim2;j++){
			t_A[j][i]=A[i][j];
		}
	}
}

//テスト用表示
void test_printf(int dim1,int dim2,double **A,char *name){
	int i,j;
	for(i=0;i<dim1;i++){
		for(j=0;j<dim2;j++){
			printf("%s[%d,%d] = %e \n",name,i,j,A[i][j]);
		}	
	}
	printf("\n");
}

f:id:takao_science:20200320111613p:plain