#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#define GNUPLOT_PATH "C:/gnuplot/bin/pgnuplot.exe"

/*定数*********************************************************************/

#define PI 3.14159265358979
#define c 2.99792458e+8	//[m/s]
#define G 6.67259e-11	//[m^3/(s^2 kg) ]

#define M 1.9891e+37	//太陽質量の10^7倍[kg]
#define q0 1.0			//二つの質量は同じ
#define kappa 0.05		//重力波エネルギーは質量エネルギーの5％

static double r_S = 2.0*G*M/(c*c);								//2.953518e+10
static double L_peak = c*c*c*c*c/G/1000.0*q0*q0;				//3.629185e+49
static double t_GW = kappa*M*c*c/(c*c*c*c*c/G/1000.0*q0*q0)/3;	//2.462969e+03 kappa*M*c*c/L_peak
static double t1 = 7.0*2.0*G*M/(c*c)/c;//kappa*M*c*c/(c*c*c*c*c/G/1000.0*q0*q0);		//2.462969e+03 t_GW
static double r_max = 100000.0*2.0*G*M/(c*c);					//100000*r_S
static double r_min = 100.0*2.0*G*M/(c*c);						//100*r_S
static double theta = 0.0;

/*関数***********************************************************************/

//GWのルミノシティ
double L_GW(double t){
	if(t<0){
		return L_peak / pow(fabs(q0*t-t1)/t1, 5/4);
	}else if(0<=t && t<t1){
		return L_peak;
	}else if(t>=t1){
		return L_peak * exp(-c*(t-t1)/2.5/r_S);
	} }

//遅延時間
double t_ret(double t, double r, double phi){
	return t - r/c*(1-sin(theta)*cos(phi)); }

//冷却時間
double t_c(double r){
	return 100; }

//GW加熱で発生するエネルギー
double He_heat(double t, double r, double phi){
	return 8.0/3.0 *G/(c*c*c) *5.0/2.0 *L_GW(t_ret(t,r,phi))/(4.0*PI*r*r); }

//GW加熱で発生したエネルギーの放射フラックス（すぐに放射する場合）
double F_GW1(double t, double r, double phi ){
	return He_heat(t, r, phi); }

//上記の関数にrをかけたもの
double F_GW1_I(double t, double r, double phi ){
	return F_GW1(t, r, phi)*r; }
/*
double Int1(double t){
	return	(4.0*t1)/(q0*t_c(r))
			*(
				1.0/pow(1.0-q0*t/t1, 0.25)
				+ 4.0/3.0*t1/q0/t_c(r)*pow(1.0-q0*t/t1, 0.75)
				- 4.0/3.0*t1/q0/t_c(r)*pow(1.0-q0*(t-t_c(r))/t1, 0.75)
			); }

//GW加熱で発生したエネルギーの放射フラックス（冷却時間で放射する場合）
double F_GW2(double t, double r, double phi ){
	if(t<0){
		return	8.0/3.0 *G/(c*c*c) *5.0/2.0 *L_peak/(4.0*PI*r*r) *Int1(t) *exp(-t/t_c(r));
	}else if(0<=t && t<t1){
		return	8.0/3.0 *G/(c*c*c) *5.0/2.0 *L_peak/(4.0*PI*r*r) *(Int1(0) + exp(t/t_c(r)) - 1.0) *exp(-t/t_c(r));
	}else if(t>=t1){
		return	8.0/3.0 *G/(c*c*c) *5.0/2.0 *L_peak/(4.0*PI*r*r) *(Int1(0) + exp(t1/t_c(r)) - 1.0 + 1.0/(1.0/t_c(r) - c/2.5/r_S)*(exp()) ) *exp(-t/t_c(r));
	} }

//上記の関数にrをかけたもの
double F_GW2_I(double t, double r, double phi ){
	return F_GW2(t, r, phi)*r; }
*/
//降着円盤から発生するエネルギー
double F_disk(double r){
	return 3.0/8.0/PI *G*M/(r*r*r); }

//上記の関数にrをかけたもの
double F_disk_I(double r){
	return F_disk(r)*r; }

/*メイン*********************************************************************/

int main(){
	printf("t_GW=%le r_S=%le L_peak=%le \n", t_GW, r_S, L_peak);
	int i, j, k, N, D;
	double  h1, h2, S_GW1, S_GW2, S_disk, t;

	printf("How about N ? （積分の分割数 1000とかで）\n");
	scanf("%d",&N);
	printf("N=%d \n" ,N);

	printf("How about D ? （時間の刻み数 10とかで）\n");
	scanf("%d",&D);
	printf("D=%d \n" ,D);

	h1 = (r_max-r_min)/N;
	h2 = 2*PI/N;

//S_diskの計算（一重積分）
	S_disk = 0.0;
	for(i=0; i<N; i=i+1){
		S_disk = S_disk + h1*F_disk_I(r_min+i*h1);
	}
	S_disk = 2 * PI * S_disk;	

//date.datファイル初期化
	FILE *fp;
	if( ( fp=fopen("data.dat", "wt") ) == NULL ) {
		printf( "ファイルがオープンできません\n" );
		exit( EXIT_FAILURE );
	}
	fclose(fp);

	for(k=-D; k<D+1; k=k+1){
		S_GW1 = 0.0;
		S_GW2 = 0.0;
	//時間をlogスケールで区切る
		t = 0.0;
		if(k<0){		t = -pow(10.0, -4.0*(double)k/(double)D);
		}else if(k==0){	t = 0;
		}else if(k>0){	t = pow(10.0, 4.0*(double)k/(double)D);
		}
	//S_GWの計算（二重積分）
		for(i=0; i<N; i=i+1){
			for(j=0; j<N; j=j+1){
				S_GW1 = S_GW1 + h1*h2*F_GW1_I(r_S/c*t, r_min+i*h1, j*h2);
				//S_GW2 = S_GW2 + h1*h2*F_GW2_I(r_S/c*t, r_min+i*h1, j*h2);
			}
		}
		//printf("t=%8.1lf S_GW1=%le S_GW2=%le S_disk=%le rate1=%lf rate2=%lf \n" , t, S_GW1, S_GW2, S_disk, S_GW1/S_disk, S_GW2/S_disk);
		printf("t=%8.1lf S_GW1=%le S_disk=%le rate1=%lf \n" , t, S_GW1, S_disk, S_GW1/S_disk);
	//data.datファイル書き込み
		FILE *fp;
		if( ( fp=fopen("data.dat", "a+t") ) == NULL ) {
			printf( "ファイルがオープンできません\n" );
			exit( EXIT_FAILURE );
		}
		fprintf(fp, "%5.2lf  %lf \n" , 4.0*(double)k/(double)D, S_GW1/S_disk);
		fclose(fp);
	}
//gnuplotでグラフ描画
	FILE *gp;
	gp = popen(GNUPLOT_PATH" -persist","w");
	fprintf(gp, "set logscale y \n");
	fprintf(gp, "plot \"data.dat\" using 1:2 \n");
	fprintf(gp, "replot \"data.dat\" using 1:2 with lines \n");
	pclose(gp);
}