/**
 * 想定した推力-重量プロフィル描画データを得る
 *  Two-stage booster
 *  Payload = 1000 [Kg]
 *  Specific impuluse = 350 [s] (both stage)
 *  fuel mass fraction = 0.73  (both stage)
 *  desired velocity = 6200 [m/s] (#1stage 2067 :(1/3),#2stage 4133 :(2/3)
 *  Ballistic Missile total wight = 64816 [Kg]
 **/


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

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

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

	double xisp1,xisp2;			/* specific impulse [s] (stage 1 & stage 2) */
	double xmf1,xmf2;			/* fuel mass fraction (stage 1 & stage 2) */

	double wpay;				/* payload [Kg] */

	double delv,delv1,delv2;	/* ΔV,ΔV1,ΔV2 [m/s] */
	double amax1,amax2;			/* maximum axial acceleration [m/s^2] */

	double numer2,denom2,wp2;	/* 第２段推進剤質量(wp2)計算式の分子,分母 */
	double ws2;					/* structural weight of 2nd-stage [Kg] */
	double wtot2;				/* weight of 2nd-stage + payload [Kg] */
	double trst2;				/* thrust level of 2nd-stage {N] */
	double tb2;					/* burn time of 2nd-stage [s] */

	double numer1,denom1,wp1;	/* 第１段推進剤質量(wp1)計算式の分子,分母 */
	double ws1;					/* structural weight of 1st-stage [Kg] */
	double wtot;				/* total ballistic missile weight [Kg] */
	double trst1;				/* thrust level of 1st-stage [N] */
	double tb1;					/* burn time of 1st-stage [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] */
	double trst;				/* 推力 [N] */

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

	xisp1 = 350;
	xisp2 = 350;

	xmf1 = 0.73;
	xmf2 = 0.73;

	wpay = 1000;

	delv = 6200;
	delv1 = 0.3333*delv;
	delv2 = 0.6667*delv;

	amax1 = 10.0;
	amax2 = 10.0;

	numer2 = wpay * (exp(delv2 / (xisp2 * G)) - 1.0);
	denom2 = 1 / xmf2 - ((1.0 - xmf2) / xmf2) * exp(delv2 / (xisp2 * G));
	wp2 = numer2 / denom2;

	ws2 = wp2 * (1 - xmf2) / xmf2;
	wtot2 = wp2 + ws2 + wpay;

	trst2 = amax2 * (wpay + ws2);
	tb2 = xisp2 * wp2 / trst2;

	numer1 = wtot2 * (exp( delv1 / (xisp1 * G)) - 1.0);
	denom1 = 1 / xmf1 - ((1.0 - xmf1) / xmf1) * exp(delv1 / (xisp1 * G));
	wp1 = numer1 / denom1;

	ws1 = wp1 * (1 - xmf1) / xmf1;
	wtot = wp1 + ws1 + wtot2;

	trst1 = amax1 * (wtot2 + ws1);
	tb1 = xisp1 * wp1 / trst1;

	delvk = delv / 1000.;

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

	while( t <= ( tb1 + tb2 ) )
	{
 		vold = v;

		step = 1;

		do
		{
			if(t < tb1)
			{
				wgt = -wp1 * t / tb1 + wtot;
				trst = trst1;
			}
			else if(t <(tb1 + tb2))
			{
				wgt = -wp2 * t / tb2 + wtot2 + wp2 *tb1 / tb2;
				trst = trst2;
			}
			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;
}