「数値計算法レポ2ルンゲクッタ版(OK)」の編集履歴(バックアップ)一覧はこちら
追加された行は緑色になります。
削除された行は赤色になります。
-正直,[[NSMRのやつ>http://www47.atwiki.jp/cscd/pages/43.html]]のほうがいいプログラムだと思う。配列使ってないという意味で。
-[[ファイル名の自動化]]を使用している。
-後半はgnuplot用のマクロファイルを作っていて,それをgvで表示するところまで入れてある。
#highlight(linenumber,c){{
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
int i,j;
double x[9999],y[9999];
double kx11[9999],kx12[9999],kx13[9999],kx14[9999];
double ky11[9999],ky12[9999],ky13[9999],ky14[9999];
double kx21[9999],kx22[9999],kx23[9999],kx24[9999];
double ky21[9999],ky22[9999],ky23[9999],ky24[9999];
double vx[9999],vy[9999];
double r,dt,p;
FILE *fpOUT;
FILE *fpPLT;
double Z;//ターゲットの原子番号
double m,a,c,hbarc;//入射粒子の質量,微細構造定数,hバーかける光速
double A;
double ax,ay;
// http://questionbox.jp.msn.com/qa1138947.html
FILE *my_fopenw(int x)
{
// int x;
if (x >= 0 /* && x < 2 */) {
char file[256];
sprintf(file, "./RK%d.dat",x);
return fopen(file, "w");
}
return (FILE *)NULL;
}
FILE *my_fopenr(int x)
{
// int x;
if (x >= 0 /* && x < 2 */) {
char file[256];
sprintf(file, "./RK%d.dat",x);
return fopen(file, "r");
}
return (FILE *)NULL;
}
int main(){
//パラメーター,定数 はじまり
dt = 1.0;//zs
c = 3.0*100.0; // speed of light 3E08 fm/zs
p = 5.3;//MeV/c
a = 1.0/137.0;//fine constructure constant
hbarc = 197.0; //[MeV*fm]=197MeV*fm
m = 3.727*1000.0; //α粒子の質量=3.727379109(93) GeV/c2
Z = 79.0; //for Au nuclear
A = Z*a*hbarc/m;
x[0] = -2000.0;//fm
vx[0]=p*c/m;
vy[0]=0;
j=0;
//おわり
// 衝突係数の入力
printf("Enter the impact parameter (fm):\n");
scanf("%lf",&y[0]);
// ファイル名の設定↓
while((fpOUT = my_fopenr(j)) != NULL){
j++;
}
fpOUT = my_fopenw(j);
printf("Create new file RK%d.dat\n",j);
// ファイル名の設定↑
fprintf(fpOUT,"%lf %lf\n",x[0],y[0]);
for(i = 1;i <= 9999;i++){
r = sqrt(x[i-1]*x[i-1]+y[i-1]*y[i-1]);
kx11[i-1] = A*x[i-1]/(r*r*r);
kx21[i-1] = vx[i-1];
kx12[i-1] = A*(x[i-1]+kx21[i-1]*dt/2)/(r*r*r);
kx22[i-1] = vx[i-1]+kx11[i-1]*dt/2;
kx13[i-1] = A*(x[i-1]+kx22[i-1]*dt/2)/(r*r*r);
kx23[i-1] = vx[i-1]+kx12[i-1]*dt/2;
kx14[i-1] = A*(x[i-1]+kx23[i-1]*dt)/(r*r*r);
kx24[i-1] = vx[i-1]+kx12[i-1]*dt;
vx[i] = (kx11[i-1]+2*kx12[i-1]+2*kx13[i-1]+kx14[i-1])*dt/6 + vx[i-1];
x[i] = (kx21[i-1]+2*kx22[i-1]+2*kx23[i-1]+kx24[i-1])*dt/6 + x[i-1];
ky11[i-1] = A*y[i-1]/(r*r*r);
ky21[i-1] = vy[i-1];
ky12[i-1] = A*(y[i-1]+ky21[i-1]*dt/2)/(r*r*r);
ky22[i-1] = vy[i-1]+ky11[i-1]*dt/2;
ky13[i-1] = A*(y[i-1]+ky22[i-1]*dt/2)/(r*r*r);
ky23[i-1] = vy[i-1]+ky12[i-1]*dt/2;
ky14[i-1] = A*(y[i-1]+ky23[i-1]*dt)/(r*r*r);
ky24[i-1] = vy[i-1]+ky12[i-1]*dt;
vy[i] = (ky11[i-1]+2*ky12[i-1]+2*ky13[i-1]+ky14[i-1])*dt/6 + vy[i-1];
y[i] = (ky21[i-1]+2*ky22[i-1]+2*ky23[i-1]+ky24[i-1])*dt/6 + y[i-1];
fprintf(fpOUT,"%lf %lf\n",x[i],y[i]);
}
fclose(fpOUT);
// Plot with GNUPLOT
fpPLT = fopen("RK.plt", "w");
fprintf(fpPLT,"set term post eps colour\n");
fprintf(fpPLT,"set output \"RK.eps\"\n");
fprintf(fpPLT,"set xlabel \"x [fm]\"\n");
fprintf(fpPLT,"set ylabel \"y [fm]\"\n");
fprintf(fpPLT,"set key below\n");
fprintf(fpPLT,"plot [-2000:2000] [-2000:2000]\"RK0.dat\" w l");
for(i=1;i<=j;i++){
fprintf(fpPLT,", \"RK%d.dat\" w l", i);
}
fprintf(fpPLT,"\n");
fprintf(fpPLT,"set term wxt");
fclose(fpPLT);
system("gnuplot RK.plt");
system("gv RK.eps");
}
}}
#image(RK.gif,title=result)
コメント募集中
#comment_num2()
-正直,[[NSMRのやつ>http://www47.atwiki.jp/cscd/pages/43.html]]のほうがいいプログラムだと思う。配列使ってないという意味で。
-[[ファイル名の自動化]]を使用している。
-後半はgnuplot用のマクロファイルを作っていて,それをgvで表示するところまで入れてある。
#highlight(linenumber,c){{
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
int i,j;
double x[9999],y[9999];
double kx11[9999],kx12[9999],kx13[9999],kx14[9999];
double ky11[9999],ky12[9999],ky13[9999],ky14[9999];
double kx21[9999],kx22[9999],kx23[9999],kx24[9999];
double ky21[9999],ky22[9999],ky23[9999],ky24[9999];
double vx[9999],vy[9999];
double r,dt,p;
FILE *fpOUT;
FILE *fpPLT;
double Z;//ターゲットの原子番号
double m,a,c,hbarc;//入射粒子の質量,微細構造定数,hバーかける光速
double A;
double ax,ay;
// http://questionbox.jp.msn.com/qa1138947.html
FILE *my_fopenw(int x)
{
// int x;
if (x >= 0 /* && x < 2 */) {
char file[256];
sprintf(file, "./RK%d.dat",x);
return fopen(file, "w");
}
return (FILE *)NULL;
}
FILE *my_fopenr(int x)
{
// int x;
if (x >= 0 /* && x < 2 */) {
char file[256];
sprintf(file, "./RK%d.dat",x);
return fopen(file, "r");
}
return (FILE *)NULL;
}
int main(){
//パラメーター,定数 はじまり
dt = 1.0;//zs
c = 3.0*100.0; // speed of light 3E08 fm/zs
p = 5.3;//MeV/c
a = 1.0/137.0;//fine constructure constant
hbarc = 197.0; //[MeV*fm]=197MeV*fm
m = 3.727*1000.0; //α粒子の質量=3.727379109(93) GeV/c2
Z = 79.0; //for Au nuclear
A = Z*a*hbarc/m;
x[0] = -2000.0;//fm
vx[0]=p*c/m;
vy[0]=0;
j=0;
//おわり
// 衝突係数の入力
printf("Enter the impact parameter (fm):\n");
scanf("%lf",&y[0]);
// ファイル名の設定↓
while((fpOUT = my_fopenr(j)) != NULL){
j++;
}
fpOUT = my_fopenw(j);
printf("Create new file RK%d.dat\n",j);
// ファイル名の設定↑
fprintf(fpOUT,"%lf %lf\n",x[0],y[0]);
for(i = 1;i <= 9999;i++){
r = sqrt(x[i-1]*x[i-1]+y[i-1]*y[i-1]);
kx11[i-1] = A*x[i-1]/(r*r*r);
kx21[i-1] = vx[i-1];
kx12[i-1] = A*(x[i-1]+kx21[i-1]*dt/2)/(r*r*r);
kx22[i-1] = vx[i-1]+kx11[i-1]*dt/2;
kx13[i-1] = A*(x[i-1]+kx22[i-1]*dt/2)/(r*r*r);
kx23[i-1] = vx[i-1]+kx12[i-1]*dt/2;
kx14[i-1] = A*(x[i-1]+kx23[i-1]*dt)/(r*r*r);
kx24[i-1] = vx[i-1]+kx12[i-1]*dt;
vx[i] = (kx11[i-1]+2*kx12[i-1]+2*kx13[i-1]+kx14[i-1])*dt/6 + vx[i-1];
x[i] = (kx21[i-1]+2*kx22[i-1]+2*kx23[i-1]+kx24[i-1])*dt/6 + x[i-1];
ky11[i-1] = A*y[i-1]/(r*r*r);
ky21[i-1] = vy[i-1];
ky12[i-1] = A*(y[i-1]+ky21[i-1]*dt/2)/(r*r*r);
ky22[i-1] = vy[i-1]+ky11[i-1]*dt/2;
ky13[i-1] = A*(y[i-1]+ky22[i-1]*dt/2)/(r*r*r);
ky23[i-1] = vy[i-1]+ky12[i-1]*dt/2;
ky14[i-1] = A*(y[i-1]+ky23[i-1]*dt)/(r*r*r);
ky24[i-1] = vy[i-1]+ky12[i-1]*dt;
vy[i] = (ky11[i-1]+2*ky12[i-1]+2*ky13[i-1]+ky14[i-1])*dt/6 + vy[i-1];
y[i] = (ky21[i-1]+2*ky22[i-1]+2*ky23[i-1]+ky24[i-1])*dt/6 + y[i-1];
fprintf(fpOUT,"%lf %lf\n",x[i],y[i]);
}
fclose(fpOUT);
// Plot with GNUPLOT
fpPLT = fopen("RK.plt", "w");
fprintf(fpPLT,"set term post eps colour\n");
fprintf(fpPLT,"set output \"RK.eps\"\n");
fprintf(fpPLT,"set xlabel \"x [fm]\"\n");
fprintf(fpPLT,"set ylabel \"y [fm]\"\n");
fprintf(fpPLT,"set key below\n");
fprintf(fpPLT,"plot [-2000:2000] [-2000:2000]\"RK0.dat\" w l");
for(i=1;i<=j;i++){
fprintf(fpPLT,", \"RK%d.dat\" w l", i);
}
fprintf(fpPLT,"\n");
fprintf(fpPLT,"set term wxt");
fclose(fpPLT);
system("gnuplot RK.plt");
system("gv RK.eps");
}
}}
実行結果の例
#image(RK.gif,title=result)
コメント募集中
#comment_num2()