「数値計算法レポ2オイラー版(OK)」の編集履歴(バックアップ)一覧はこちら
「数値計算法レポ2オイラー版(OK)」(2011/07/08 (金) 10:19:15) の最新版変更点
追加された行は緑色になります。
削除された行は赤色になります。
正直,[[NSMRのやつ>http://www47.atwiki.jp/cscd/pages/43.html]]のほうがいいプログラムだと思う。配列使ってないという意味で。詳しくは[[数値計算法レポ2ルンゲクッタ版(OK)]]で。
#highlight(linenumber,c){{
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
int i;
double x[999],y[999],k1[999],k2[999],k3[999],k4[999];
double vx[999],vy[999];
double r,dt,p;
FILE *fpOUT;
FILE *fpPLT;
double Z;//ターゲットの原子番号
double m,a,c,hbarc;//入射粒子の質量,微細構造定数,hバーかける光速
double A;
double ax,ay;
int j;
// http://questionbox.jp.msn.com/qa1138947.html
FILE *my_fopenw(int x)
{
// int x;
if (x >= 0 /* && x < 2 */) {
char file[256];
sprintf(file, "./euler%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, "./euler%d.dat",x);
return fopen(file, "r");
}
return (FILE *)NULL;
}
int main(){
//パラメーター,定数
dt = 10.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 euler%d.dat\n",j);
// ファイル名の設定↑
fprintf(fpOUT,"%lf %lf\n",x[0],y[0]);
for(i = 1;i <= 999;i++){
r = sqrt(x[i-1]*x[i-1]+y[i-1]*y[i-1]);
ax = A*x[i-1]/(r*r*r);
vx[i] = ax * dt + vx[i-1];
x[i] = vx[i] * dt + x[i-1];
ay = A*y[i-1]/(r*r*r);
vy[i] = ay *dt + vy[i-1];
y[i] = vy[i] * dt + y[i-1];
// printf("%lf %lf %lf %lf %lf %lf\n",x[i],y[i],vx[i],vy[i],ax,ay);
// デバッグ用↑
fprintf(fpOUT,"%lf %lf\n",x[i],y[i]);
}
fclose(fpOUT);
// Plot with GNUPLOT
fpPLT = fopen("euler.plt", "w");
fprintf(fpPLT,"set term post eps colour\n");
fprintf(fpPLT,"set output \"euler.eps\"\n");
fprintf(fpPLT,"set xlabel \"x [fm]\"\n");
fprintf(fpPLT,"set ylabel \"[fm]\"\n");
fprintf(fpPLT,"set nokey\n");
//fprintf(fpPLT,"\n");
fprintf(fpPLT,"plot [-2000:2000] [-2000:2000]\"euler0.dat\" w l");
for(i=1;i<=j;i++){
fprintf(fpPLT,", \"euler%d.dat\" w l", i);
}
fprintf(fpPLT,"\n");
fprintf(fpPLT,"set term wxt");
fclose(fpPLT);
system("gnuplot euler.plt");
system("gv euler.eps");
}}}
コメント募集中
#comment_num2()
正直,[[NSMRのやつ>http://www47.atwiki.jp/cscd/pages/43.html]]のほうがいいプログラムだと思う。配列使ってないという意味で。詳しくは[[数値計算法レポ2ルンゲクッタ版(OK)]]で。
#highlight(linenumber,c){{
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
int i;
double x[999],y[999],k1[999],k2[999],k3[999],k4[999];
double vx[999],vy[999];
double r,dt,p;
FILE *fpOUT;
FILE *fpPLT;
double Z;//ターゲットの原子番号
double m,a,c,hbarc;//入射粒子の質量,微細構造定数,hバーかける光速
double A;
double ax,ay;
int j;
// http://questionbox.jp.msn.com/qa1138947.html
FILE *my_fopenw(int x)
{
// int x;
if (x >= 0 /* && x < 2 */) {
char file[256];
sprintf(file, "./euler%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, "./euler%d.dat",x);
return fopen(file, "r");
}
return (FILE *)NULL;
}
int main(){
//パラメーター,定数
dt = 10.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 euler%d.dat\n",j);
// ファイル名の設定↑
fprintf(fpOUT,"%lf %lf\n",x[0],y[0]);
for(i = 1;i <= 999;i++){
r = sqrt(x[i-1]*x[i-1]+y[i-1]*y[i-1]);
ax = A*x[i-1]/(r*r*r);
vx[i] = ax * dt + vx[i-1];
x[i] = vx[i] * dt + x[i-1];
ay = A*y[i-1]/(r*r*r);
vy[i] = ay *dt + vy[i-1];
y[i] = vy[i] * dt + y[i-1];
// printf("%lf %lf %lf %lf %lf %lf\n",x[i],y[i],vx[i],vy[i],ax,ay);
// デバッグ用↑
fprintf(fpOUT,"%lf %lf\n",x[i],y[i]);
}
fclose(fpOUT);
// Plot with GNUPLOT
fpPLT = fopen("euler.plt", "w");
fprintf(fpPLT,"set term post eps colour\n");
fprintf(fpPLT,"set output \"euler.eps\"\n");
fprintf(fpPLT,"set xlabel \"x [fm]\"\n");
fprintf(fpPLT,"set ylabel \"[fm]\"\n");
fprintf(fpPLT,"set nokey\n");
//fprintf(fpPLT,"\n");
fprintf(fpPLT,"plot [-2000:2000] [-2000:2000]\"euler0.dat\" w l");
for(i=1;i<=j;i++){
fprintf(fpPLT,", \"euler%d.dat\" w l", i);
}
fprintf(fpPLT,"\n");
fprintf(fpPLT,"set term wxt");
fclose(fpPLT);
system("gnuplot euler.plt");
system("gv euler.eps");
}
}}
コメント募集中
#comment_num2()