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

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

#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％
#define theta 0.0

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);	//2.462969e+03 kappa*M*c*c/L_peak
static double t1 = 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

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

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 He_heat(double t, double r, double phi){
	return 8/3 *G/(c*c*c) *5/2 *L_GW(t_ret(t, r, phi))/(4*PI*r*r);
}

double F_GW(double t, double r, double phi ){
	return He_heat(t, r, phi);
}

double F_disk(double r){
	return 3.0/8.0/PI *G*M/(r*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_GW, 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 = 0.0;
	for(i=0; i<N; i=i+1){
		for(j=0; j<N; j=j+1){
			S_disk = S_disk + h1*h2*F_disk(r_min+i*h1)*(r_min+i*h1);
		}
	}

	for(k=-D; k<D+1; k=k+1){
		S_GW = 0.0;
		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);
		}

		for(i=0; i<N; i=i+1){
			for(j=0; j<N; j=j+1){
				S_GW = S_GW + h1*h2*F_GW(r_S/c*t, r_min+i*h1, j*h2)*(r_min+i*h1);
			}
		}
		printf("t=%10.3lf S_GW=%le S_disk=%le rate=%lf \n" ,t, S_GW, S_disk, fabs(S_GW/S_disk));

	}


}