数値計算法レポ2ルンゲクッタ版(OK)


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

  • 正直,NSMRのやつのほうがいいプログラムだと思う。配列使ってないという意味で。
  • ファイル名の自動化を使用している。
  • 後半はgnuplot用のマクロファイルを作っていて,それをgvで表示するところまで入れてある。

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

実行結果の例

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

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