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

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

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 fx(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;

void main()
{
  double xr[N];
  double x[NMAX], y[NMAX], z[NMAX];             /*  粒子の座標 */
  double vx[NMAX], vy[NMAX], vz[NMAX];          /* 粒子の速度 */
  double t, dt;                                 /* 時間, 時間ステップ */
  float xplot[NMAX], yplot[NMAX], zplot[NMAX];  /* for plot */
  double x1, x2, y1, y2;
  double xg1, yg1, zg1;
  double xg2, yg2, zg2;
  double vxg2, vyg2, vzg2;
  double mg1, mg2, eps, eps2;
  double dx1, dy1, dz1;
  double dx2, dy2, dz2;
  double r1, r2, r13, r23;
  double xp, yp, zp;
  double pi;
  double r, rg, vcirc;
  double ekin, epot, etot;
  int npar, nstep, np, i, istep, iseed;

  /* 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=-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;
      x[np]=xp;
      y[np]=yp;
      z[np]=zp;
      /*      vcirc=sqrt(mg1/r); */
      vcirc=sqrt(mg1/sqrt(r*r+eps2));
      vx[np]=-vcirc*y[np]/r;
      vy[np]=vcirc*x[np]/r;
      vz[np]=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++) {
      dx1=x[i]-xg1;
      dy1=y[i]-yg1;
      dz1=z[i]-zg1;
      r1=sqrt(dx1*dx1+dy1*dy1+dz1*dz1+eps2);
      dx2=x[i]-xg2;
      dy2=y[i]-yg2;
      dz2=z[i]-zg2;
      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;
      ekin=0.5*(vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i]);
      epot=-mg1/r1-mg2/r2;
      etot=etot+ekin+epot;

      xr[0]=x[i];
      xr[1]=y[i];
      xr[2]=z[i];
      xr[3]=vx[i];
      xr[4]=vy[i];
      xr[5]=vz[i];

      rungen(t, xr, dt);

      x[i]=xr[0];
      y[i]=xr[1];
      z[i]=xr[2];
      vx[i]=xr[3];
      vy[i]=xr[4];
      vz[i]=xr[5];
      xplot[i]=x[i];
      yplot[i]=y[i];
      zplot[i]=z[i];
    }

    printf("istep= %d,  etot= %lf\n", istep, etot);
      
    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;
}

/**********************************************************************/ 
/*   3次元運動方程式を解くためのルンゲ・クッタ法のサブルーチン */
void rungen(double t, double x[], double step)
{
  double f(double t, double x[], int i);
  double h=step/2.0,                         /* the midpoint */
         t1[N], t2[N], t3[N],                /* temporary storage */
         k1[N], k2[N], k3[N],k4[N];          /* for Runge-Kutta  */
  int i;
 
  for (i=0; i<N; i++) t1[i] = x[i]+0.5*(k1[i]=step*f(t, x, i));
  for (i=0; i<N; i++) t2[i] = x[i]+0.5*(k2[i]=step*f(t+h, t1, i));
  for (i=0; i<N; i++) t3[i] = x[i]+    (k3[i]=step*f(t+h, t2, i));
  for (i=0; i<N; i++) k4[i] =                 step*f(t + step, t3, i);
  
  for (i=0; i<N; i++) x[i] += (k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
}
/*--------------------------------------------------------------------*/

/* definition of equations - this is the harmonic oscillator */
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*/
}