/*  銀河の衝突のシミュレーション */

#define N 6
#define NMAX 10000

#include "cpgplot.h"

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* void rungen(double t, double x[], double dt); */

void  force(double dt, double x, double y, double z) ;
double  f(double t, double x[N], int i);

/*  double fx(double t,double x,double y,double z,double vx,double vy,double vz); */
/*  double fy(double t,double x,double y,double z,double vx,double vy,double vz); */
/*  double fz(double t,double x,double y,double z,double vx,double vy,double vz); */
/*  double fvx(double t,double x,double y,double z,double vx,double vy,double vz); */
/*  double fvy(double t,double x,double y,double z,double vx,double vy,double vz); */
/*  double fvz(double t,double x,double y,double z,double vx,double vy,double vz); */

double fx1, fy1, fz1;
double xg1, yg1, zg1;
double xg2, yg2, zg2;
double vxg2, vyg2, vzg2;
double mg1, mg2, eps, eps2;

void main()
{
  double star[NMAX][6];             /*  粒子の座標 */
  double t, dt;                                 /* 時間, 時間ステップ */
  float xplot[NMAX], yplot[NMAX], zplot[NMAX];  /* for plot */
  double x1, x2, y1, y2;
  double xp, yp, zp;
  double pi;
  double r, rg, vcirc;
  double ekin, epot, etot;
  int npar, nstep, np, i, istep, iseed, j;
  double t0[N], t1[N], t2[N], t3[N],                /* temporary storage */
         k1[N], k2[N], k3[N],k4[N];          /* for Runge-Kutta  */
  double starx,stary,starz;

  /* mg1, mg2: 銀河の質量
     fx1, fy1, fz1: 粒子に働く力 */

  char STR[8];
  
  printf("粒子数: npar ?\n");
  scanf("%d", &npar);
  printf("計算ステップ数: nstep ?\n");
  scanf("%d", &nstep);
  printf("時間刻み: dt ?\n");
  scanf("%lf", &dt);
  printf("衝突銀河の質量: mg2 ?\n");
  scanf("%lf", &mg2);
  printf("衝突パラメータ: yg2 ?\n");
  scanf("%lf", &yg2);
  printf("衝突銀河の速度: vxg2 ?\n");
  scanf("%lf", &vxg2);

  pi= 3.1415926;
  mg1=1.0;
  /*  mg2=1.0; */
  eps=0.1;
  eps2=eps*eps;

  /*  初期条件の設定 */
  x1=-5.0;
  x2=5.0;
  y1=-5.0;
  y2=5.0;
      
  /*  銀河1の位置 */
  xg1=0.0;
  yg1=0.0;
  zg1=0.0;

  /*  銀河2の位置 */
  xg2=-5.0;
  /*  yg2=2.0; */
  zg2=0.0;
  /*  銀河2の速度 */
  /* vxg2=1.0; */
  vyg2=0.0; 
  vzg2=0.0;
  
  rg=1.0;
  iseed=1;
  np=-1;
  
  srand(iseed);

  /*  粒子を半径rgの円内にランダムに分布させ、回転の速度を与える */
  for(i=0; i<NMAX; i++) {
    
    xp=-rg+2*rg*rand()/(RAND_MAX+1.0);
    yp=-rg+2*rg*rand()/(RAND_MAX+1.0);
    zp=0.0;
    r=sqrt(xp*xp+yp*yp+zp*zp);
    if(r < rg) {
      np=np+1;
      star[np][0]=xp;
      star[np][1]=yp;
      star[np][2]=zp;
      /*      vcirc=sqrt(mg1/r); */
      vcirc=sqrt(mg1/sqrt(r*r+eps2));
      star[np][3]=-vcirc*star[np][1]/r;
      star[np][4]=vcirc*star[np][0]/r;
      star[np][5]=0.0;
/*        printf("x=%f, y=%f, z=%f\n", x[np], y[np], z[np]); */
    }
    if(np == npar-1) break;
  }

  if(cpgbeg(0,"?",1,1) != 1) return;

  cpgpage();
  cpgvstd();
  cpgwnad(x1, x2, y1, y2);
 
  t=0.0;
  for(istep=0; istep<nstep; istep++) {
    /*    etot=0.0; */
      
    cpgbbuf();
    cpgeras();
    cpgsci(1);
    cpgbox("bcnst", 0.0, 0, "bcnst", 0.0, 0);

    for(i=0; i<npar; i++) {

      starx=star[i][0];
      stary=star[i][1];
      starz=star[i][2];
      t0[0]=star[i][0];
      t0[1]=star[i][1];
      t0[2]=star[i][2];
      t0[3]=star[i][3];
      t0[4]=star[i][4];
      t0[5]=star[i][5];

      force(0.0, starx, stary, starz);
      for (j=0; j<N; j++) {
	k1[j] = dt*f(t, t0, j);
	t1[j] = star[i][j]+0.5*k1[j]; 
      }

      starx=t1[0];
      stary=t1[1];
      starz=t1[2];
      
      force(dt/2.0, starx, stary, starz);
      for (j=0; j<N; j++) {
	k2[j] = dt*f(t, t1, j);
	t2[j] = t1[j]+0.5*k2[j]; 
      }

      starx=t2[0];
      stary=t2[1];
      starz=t2[2];
      
      force(dt/2.0, starx, stary, starz);
      for (j=0; j<N; j++) {
	k3[j] = dt*f(t, t2, j);
	t3[j] = t2[j]+0.5*k3[j]; 
      }

      starx=t3[0];
      stary=t3[1];
      starz=t3[2];

      force(dt,     starx, stary, starz);
      for (j=0; j<N; j++) {
	k4[j] = dt*f(t, t3, j);
	  /*	t4[j] = t3[j]+0.5*k3[j]; */
      }

      for (j=0; j<N; j++) {star[i][j] += (k1[j]+2*k2[j]+2*k3[j]+k4[j])/6.0; 
      } 

      xplot[i]=star[i][0];
      yplot[i]=star[i][1];
      zplot[i]=star[i][2];
    }
    
/*      printf("istep= %d,  etot= %lf\n", istep, etot); */
    printf("istep= %d\n", istep);
      
    cpgsci(3);
    cpgslw(15);
    cpgpt(npar, xplot, yplot, -1);
    cpgsci(2);
    cpgpt1((float)xg1,(float)yg1,22);
    cpgpt1((float)xg2,(float)yg2,22);
    cpgslw(1);
    sprintf(STR, "%d", istep); 
    cpgmtxt("istep", 2.0, 0.0, 0.0, STR);
    cpgebuf();
   
    /*  時間を進める */
    t=t+dt;

    /*  銀河2を動かす */
    xg2=xg2+vxg2*dt;
    yg2=yg2+vyg2*dt;
    zg2=zg2+vzg2*dt;
    
  }
  
  cpgend();
  
  return;
}

void  force(double dt, double x, double y, double z) 
{
/*    double xg1, yg1, zg1; */
  double xg2t, yg2t, zg2t;
/*    double vxg2, vyg2, vzg2; */
/*    double mg1, mg2, eps, eps2; */
  double dx1, dy1, dz1;
  double dx2, dy2, dz2;
  double r1, r2, r13, r23;

  dx1=x-xg1;
  dy1=y-yg1;
  dz1=z-zg1;
  r1=sqrt(dx1*dx1+dy1*dy1+dz1*dz1+eps2);
  xg2t=xg2+vxg2*dt;
  yg2t=yg2+vxg2*dt;
  zg2t=zg2+vxg2*dt;
  dx2=x-xg2t;
  dy2=y-yg2t;
  dz2=z-zg2t;
  r2=sqrt(dx2*dx2+dy2*dy2+dz2*dz2+eps2);
  r13=pow(r1,3);
  r23=pow(r2,3);
  fx1=-mg1*dx1/r13-mg2*dx2/r23;
  fy1=-mg1*dy1/r13-mg2*dy2/r23;
  fz1=-mg1*dz1/r13-mg2*dz2/r23;

}


/* definition of equations */
double  f(double t, double x[], int i)
{
  if (i == 0) return(x[3]);   /* dx/dt = vx */
  if (i == 1) return(x[4]);   /* dy/dt = vy */
  if (i == 2) return(x[5]);   /* dz/dt = vz */
  if (i == 3) return(fx1);    /* dvx/dt = fx1*/
  if (i == 4) return(fy1);    /* dvy/dt = fy1*/
  if (i == 5) return(fz1);    /* dvz/dt = fz1*/
}