GaAsの平衡分圧 (Ｃ言語)

　GaNの結晶成長の定量的な計算がしたくて、参考文献を下に、練習としてガリウム砒素(GaAs)で試してみました。水素をキャリアガスとした化学気相成長(CVD)法におけるGaAsの平衡分圧です。原料はトリメチルガリウム(TMG)とAs4。

【補足】
　練習目的：Newton法を用いた微分計算
　懸念事項：doubleの桁落ちにより、Gaの情報が損失
【参考文献】
　『エピタキシャル成長のメカニズム』（中嶋一雄）


#include stdio.h
#include math.h

struct formula {
double x;
};

int main(void)
{
double T, Temp, ek1, ek2, Pg, Pa, Ch;
double As4, Ga, As2, H2;
int i, j;

FILE *ga, *a4, *a2, *h2, *ch, *gp;
char *ga_file, *a4_file, *a2_file, *h2_file, *ch_file;
ga_file="out1.dat";
a4_file="out2.dat";
a2_file="out3.dat";
h2_file="out4.dat";
ch_file="out5.dat";
ga = fopen(ga_file,"w");
a4 = fopen(a4_file,"w");
a2 = fopen(a2_file,"w");
h2 = fopen(h2_file,"w");
ch = fopen(ch_file,"w");

//600-950度の温度範囲で計算
for(j=1;j<=7;j++){
Temp=600.0+(j-1)*50;
T=273.16+Temp;

//平衡定数K1, K2
//K1=a/(PG*Pa4^(1/4))
//K2=Pa2^2/Pa4
ek1=pow(10,(-6.49+1.91*pow(10,4)/T-33.5*pow(10,-2)*log10(T)));
ek2=pow(10,(-11.4-1.10*pow(10,4)/T+5.51*log10(T)));

//原料GaとAsの分圧[atm]
Pg=5.0*pow(10,-5);
Pa=5.0*pow(10,-4);
Ch=3.0*Pg;

//Pa4^(1/4)の５次方程式 A1*x^5+A2*x^4+A3*x^3+A4*x^2+A5*x+A6=0
struct formula A[6] = {4.0, 0.0, 2*sqrt(ek2), 0.0, Pg-Pa, -1/ek1};
double x0=1.0, x1=0.0, S=0.0, U=0.0;

//Pa4をNewton法で計算
form:
for(i=0;i<=4;i++){
S=S*x0+A[i].x;
U=U*x0+S;
}
S=S*x0+A[5].x;
x1=x0-S/U;

if((x0-x1)/x0>pow(10,-10)) {
x0=x1;
S=0.0;
U=0.0;
goto form;
}

//分圧の入力
else{
As4=pow(x1,4);
Ga=1.0/(ek1*pow(As4,0.25));
As2=sqrt(ek2*As4);
H2=1-(Ga+As4+As2+Ch);
fprintf(ga, "%lf \t%lf \n",Temp, Ga*760);
fprintf(a4, "%lf \t%lf \n",Temp, As4*760);
fprintf(a2, "%lf \t%lf \n",Temp, As2*760);
fprintf(h2, "%lf \t%lf \n",Temp, H2*760);
fprintf(ch, "%lf \t%lf \n",Temp, Ch*760);
}
}
fclose(ga);
fclose(a4);
fclose(a2);
fclose(h2);
fclose(ch);

// Gnuplotで表示
gp = popen("/Applications/gnuplot.app/gnuplot -persist","w");
fprintf(gp, "set logscale y\n");
fprintf(gp, "set xlabel \"Temperature [oC]\"\n");
fprintf(gp, "set ylabel \"Equilibrium partial pressure [Torr]\"\n");
fprintf(gp, "set xrange [590:950]\n");
fprintf(gp, "set yrange [0.000001:1500]\n");
fprintf(gp, "set style line 1 lt 1 lw 3\n");
fprintf(gp, "set style line 2 lt 2 lw 3\n");
fprintf(gp, "set style line 3 lt 3 lw 3\n");
fprintf(gp, "set style line 4 lt 4 lw 3\n");
fprintf(gp, "set style line 5 lt 5 lw 3\n");
fprintf(gp, "plot \"%s\" with lines linetype 1 title \"Ga\"\n", ga_file);
fprintf(gp, "replot \"%s\" with lines linetype 2 title \"As4\"\n", a4_file);
fprintf(gp, "replot \"%s\" with lines linetype 3 title \"As2\"\n", a2_file);
fprintf(gp, "replot \"%s\" with lines linetype 4 title \"H2\"\n", h2_file);
fprintf(gp, "replot \"%s\" with lines linetype 5 title \"CH3\"\n", ch_file);
pclose(gp);
return 0;
}