#include #include #include #include #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