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

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

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