数値計算法レポ2オイラー版(OK) - 阪大理学部科学計算・開発サークル -Circle for Science Caluculation and Development-

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

#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");
}
 

コメント募集中