数値計算法レポ2オイラー版(OK)


※上記の広告は60日以上更新のないWIKIに表示されています。更新することで広告が下部へ移動します。

正直,NSMRのやつのほうがいいプログラムだと思う。配列使ってないという意味で。詳しくは数値計算法レポ2ルンゲクッタ版(OK)で。

  1. #include<stdio.h>
  2. #include<stdlib.h>
  3. #include<math.h>
  4.  
  5.  
  6. int i;
  7.  
  8. double x[999],y[999],k1[999],k2[999],k3[999],k4[999];
  9. double vx[999],vy[999];
  10. double r,dt,p;
  11. FILE *fpOUT;
  12. FILE *fpPLT;
  13. double Z;//ターゲットの原子番号
  14. double m,a,c,hbarc;//入射粒子の質量,微細構造定数,hバーかける光速
  15. double A;
  16. double ax,ay;
  17. int j;
  18.  
  19. // http://questionbox.jp.msn.com/qa1138947.html
  20. FILE *my_fopenw(int x)
  21. {
  22. // int x;
  23. if (x >= 0 /* && x < 2 */) {
  24. char file[256];
  25. sprintf(file, "./euler%d.dat",x);
  26. return fopen(file, "w");
  27. }
  28. return (FILE *)NULL;
  29. }
  30.  
  31. FILE *my_fopenr(int x)
  32. {
  33. // int x;
  34. if (x >= 0 /* && x < 2 */) {
  35. char file[256];
  36. sprintf(file, "./euler%d.dat",x);
  37. return fopen(file, "r");
  38. }
  39. return (FILE *)NULL;
  40. }
  41.  
  42. int main(){
  43.  
  44. //パラメーター,定数
  45. dt = 10.0;//zs
  46. c = 3.0*100.0; // speed of light 3E08 fm/zs
  47. p = 5.3;//MeV/c
  48. a = 1.0/137.0;//fine constructure constant
  49. hbarc = 197.0; //[MeV*fm]=197MeV*fm
  50. m = 3.727*1000.0; //α粒子の質量=3.727379109(93) GeV/c2
  51. Z = 79.0; //for Au nuclear
  52. A = Z*a*hbarc/m;
  53. x[0] = -2000.0;//fm
  54. vx[0]=p*c/m;
  55. vy[0]=0;
  56.  
  57. j=0;
  58.  
  59. //おわり
  60.  
  61. // 衝突係数の入力
  62. printf("Enter the impact parameter (fm):\n");
  63. scanf("%lf",&y[0]);
  64.  
  65.  
  66. // ファイル名の設定↓
  67. while((fpOUT = my_fopenr(j)) != NULL){
  68. j++;
  69. }
  70. fpOUT = my_fopenw(j);
  71. printf("Create new file euler%d.dat\n",j);
  72. // ファイル名の設定↑
  73.  
  74. fprintf(fpOUT,"%lf %lf\n",x[0],y[0]);
  75. for(i = 1;i <= 999;i++){
  76. r = sqrt(x[i-1]*x[i-1]+y[i-1]*y[i-1]);
  77. ax = A*x[i-1]/(r*r*r);
  78. vx[i] = ax * dt + vx[i-1];
  79. x[i] = vx[i] * dt + x[i-1];
  80.  
  81. ay = A*y[i-1]/(r*r*r);
  82. vy[i] = ay *dt + vy[i-1];
  83. y[i] = vy[i] * dt + y[i-1];
  84. // printf("%lf %lf %lf %lf %lf %lf\n",x[i],y[i],vx[i],vy[i],ax,ay);
  85. // デバッグ用↑
  86. fprintf(fpOUT,"%lf %lf\n",x[i],y[i]);
  87. }
  88. fclose(fpOUT);
  89.  
  90.  
  91. // Plot with GNUPLOT
  92.  
  93. fpPLT = fopen("euler.plt", "w");
  94. fprintf(fpPLT,"set term post eps colour\n");
  95. fprintf(fpPLT,"set output \"euler.eps\"\n");
  96. fprintf(fpPLT,"set xlabel \"x [fm]\"\n");
  97. fprintf(fpPLT,"set ylabel \"[fm]\"\n");
  98. fprintf(fpPLT,"set nokey\n");
  99. //fprintf(fpPLT,"\n");
  100.  
  101. fprintf(fpPLT,"plot [-2000:2000] [-2000:2000]\"euler0.dat\" w l");
  102. for(i=1;i<=j;i++){
  103. fprintf(fpPLT,", \"euler%d.dat\" w l", i);
  104. }
  105. fprintf(fpPLT,"\n");
  106. fprintf(fpPLT,"set term wxt");
  107. fclose(fpPLT);
  108. system("gnuplot euler.plt");
  109. system("gv euler.eps");
  110. }
  111.  

コメント募集中
名前:
コメント:

すべてのコメントを見る
添付ファイル