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

#define pai 3.14159265359
#define hc 1.97315E-8   /* keV cm */
#define cc 2.99793E+8    /* m/sec */
#define eps 1.0
#define eps0 11.83	/* silicon */
#define NA 6.02E+23     /* アボガドロ数 */
#define Me 5.11006E+2 	/* keV */
#define Alpha 7.2972E-3  /* 微細構造定数 */
#define SigT 6.65E-25  /* cm2 トムソン散乱断面積 */
#define Sig0 5.79E-28  /* cm2 σ_0/Z^2 */
#define Denergy 0.1  /* energy step */
#define Num 1000  /* step 数 */


/* a, b 小さいほう */
double fmin(double a, double b)
{
  if (a>b) return b;
  else return a;
}

/* 個数密度[1/cm3] */
double nkosu(double A, double rou)
{
  double n=1;

  n=NA*rou/A; /* [1/cm3] */
  return n;
}

/* プラズマ振動数[keV] Z:原子番号 nrou:個数密度[1/cm3] */
double hwplasma(int Z, double nrou)
{
  double n=1, hwp=1, x=1;

  n=((double) Z)*nrou; /* [1/cm3] */
  x=pow(hc,3)*n*Alpha/Me;
  hwp=4*M_PI*sqrt(x);  /* [keV] */
  return hwp;
}

/* hw[keV]に対する吸収長[cm] Z: 原子番号 nrou:個数密度[1/cm3] */ 
double absplngth(int Z, double nrou, double hw)
{
  double la=0.1, er=1.0;

  er=hw/Me;
  la=0.1767767/SigT/pow(Alpha,4)/pow(Z,5)/nrou*pow(er,3.5);

  return la;
}

int main(void)
{
	FILE *fp;
	int Z=10, j=0;
	long loop=100;
	double rou=1.0, nrou=1.0, A=20, l=1, la=1, le=1;
	double e=0.0, ekv=1, gamma=1.0, d=3.14E-10, sg=0, maxe=1.0;
	double I=1.0, Itot=0.0, Ptot=0.0;
	double hw=1.0, dhw=0.1, hwp=1, myu=1, sigma=1;
	char name[64]="";
	time_t first, second;

	printf("Output file = ? "); scanf("%s",name);
	if ((fp = fopen(name, "wt"))== NULL){
		printf("Cannot open output file.\n");
		exit(1);
	}

	printf("electron energy [MeV] = ? "); scanf("%lf",&e);
	ekv=(1.0e+3)*e; /* E [keV] */
	gamma=ekv/Me;  /* γfactor */

	printf("物質のパラメータ\n");
	printf("原子番号 = ? "); scanf("%d",&Z);
	printf("質量数 = ? "); scanf("%lf",&A);
	printf("密度 [g/cm3] = ? "); scanf("%lf",&rou);
	printf("厚さ [cm] = ? "); scanf("%lf",&l);
	printf("計算範囲の最大 [MeV] = ? \n"); scanf("%lf",&maxe);
	if(maxe>e) maxe=e;
	printf("計算範囲は 1keV -> 1MeV \n");
	printf("Δhω [keV] = ? "); scanf("%lf",&dhw);
	loop=(long) ((1.0e+3)*maxe/dhw);

	fprintf(fp, "## E = %lf [MeV]\n", e); 
	fprintf(fp, "## Z=%d, A=%lf, ρ=%lf[g/cm3], thick=%lf[cm]\n", Z, A, rou, l); 
	first = time(NULL);  /* システム時間の取得 */

	nrou=nkosu(A, rou);  /* 原子個数密度 */
	hwp=hwplasma(Z,nrou);  /* [keV] プラズマ振動数 */
	hw=1.0; Itot=0; Ptot=0;
	for (j=0; j<loop; j++){
	  hw=1.0+dhw*((double) j);  /* [keV] */
	  myu=1+pow(hwp/hw,2)*pow(gamma,2);  /* screening 補正項 */
	  la=absplngth(Z, nrou, hw);   /* [cm] 吸収長 */
	  /*	  le=fmin(l, la); */
	  le=la*(1.0-exp(-1.0*l/la));
	  sigma=16.0/3.0*pow(Z,2)*Sig0*nrou;  /* [cm2 1/cm3] */
	  I=sigma*log(183*pow(Z,-0.3333)*sqrt(myu))/myu*dhw/hw*le;
	  fprintf(fp,"%lf, %lg\n", hw, I);
	  Itot=Itot+I;
	  Ptot=Ptot+I*hw;
	  hw=hw+dhw;
	}

	second = time(NULL); /* システム時間の取得 */

	printf("total output = %lg\n",Itot);
	printf("total emission energy = %lg [keV]\n",Ptot);
	printf("calculation time is: %lf seconds\n",difftime(second,first));

	fclose(fp);
	getchar();
	exit(0);
}