/**
 * 想定した推力-重量プロフィル描画データを得る
 *  Single-stage booster
 *  Payload = 1000 [Kg]
 *  Specific impuluse = 350 [s]
 *  fuel mass fraction = 0.67
 *  desired velocity = 3400 [m/s]
 *  Ballistic Missile total wight = 16312 [Kg]
 **/


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

#define G (9.799)     /* 地球表面上における重力の加速度 [m/s^2]	*/

int main()
{
	int step;					/* 計算ステップ */

	double xisp;				/* specific impulse [s] */
	double xmf;					/* fuel mass fraction */

	double wpay;				/* payload [Kg] */

	double delv;				/* ΔV[m/s] */
	double amax;				/* maximum axial acceleration [m/s^2] */

	double numer,denom,wp;		/* 推進剤質量(wp2)計算式の分子,分母 */
	double ws;					/* structural weight [Kg] */
	double wtot;				/* Wp + Ws + payload [Kg] */
	double trst;				/* thrust level {N] */
	double tb;					/* burn time [s] */

	double delvk;				/* ΔV / 1000 ([km/s] 表現に変換) */

	double h;					/* 積分刻み幅 [s] */

	double t;					/* ブースト時間 [s] */

	double dout;				/* データ出力間隔 [s] */

	double v,vold;				/* ミサイル速度 [m/s], 新,旧 */

	double acc,accg;			/* ミサイル加速度 [m/s^2],重力(g)表現 */

	double vk;					/* v / 1000 ([km/s] 表現に変換) */
	double wgt;					/* ミサイル質量 [Kg] */

	FILE *fptr;
	fptr = fopen("DATFIL.TXT","w");

	xisp = 350;

	xmf = 0.67;

	wpay = 1000;

	delv = 3400;

	amax = 10.0;

	numer = wpay * (exp(delv / (xisp * G)) - 1.0);
	denom = 1 / xmf - ((1.0 - xmf) / xmf) * exp(delv / (xisp * G));
	wp = numer / denom;

	ws = wp * (1 - xmf) / xmf;
	wtot = wp + ws + wpay;

	trst = amax * (wpay + ws);
	tb = xisp * wp / trst;

	delvk = delv / 1000.0;

	h = 0.01;
	t = 0.0;
	dout = 0.0;
	v = 0.0;

	while( t <= (tb) )
	{
 		vold = v;

		step = 1;

		do
		{
			if(t < tb)
			{
				wgt = -wp * t / tb + wtot;
			}
			else
			{
				wgt = wpay;
				trst = 0.0;
			}

			acc = G * trst / wgt;

			switch(step)
			{
				case 1:
					step = 2;
					v = v + h * acc;
					t = t + h;
					break;

				case 2:
					step = 1;
					v = (vold + v) / 2 + 0.5 * h * acc;
 					dout = dout + h;

					if(dout >= 0.99999)
					{
 						dout = 0.0;
						accg = acc / G;
						vk = v / 1000.0;

						printf("%8.2f %8.2f %8.2f\n",t,v,accg);
						fprintf(fptr,"%8.2f %8.2f %8.2f\n",t,v,accg);
					}

					break;

				default:;

			} /* switch */

		}while( step >= 2 ); /* do - while */

	} /* while */

 	accg = acc / G;
	vk = v / 1000.;

 	printf("%8.2f %8.2f %8.2f\n",t,v,accg);

	fprintf(fptr,"%8.2f %8.2f %8.2f\n",t,v,accg);
	fclose (fptr);

	return 0;
}