//////////////////////////////////
// 荷電粒子の磁場中での運動
// ①一様磁場中(-5.0 < z < 5.0)
// ②
// ③
// 2013/08/22-2013/09/16   T. Kojima
//////////////////////////////////

#include <GLUT/glut.h>
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <string.h> //文字列の長さを調べるため

#define ix   0.333  //荷電粒子の初期位置のx座標(固定)    |vy/Bz|とすればいい
#define iy   0.0    //荷電粒子の初期位置のy座標(固定)
#define iz   5.0    //荷電粒子の初期位置のz座標(固定)
#define ivx  0.0    //荷電粒子の初期速度のx成分
#define ivy  1.0    //荷電粒子の初期速度のy成分
#define ivz  (-0.5) //荷電粒子の初期速度のz成分
#define dt 0.05   //サンプリング時間
#define MAXPOINT 10000 //微小線の数


////Global variables////  回転半径r=v/Bで表せる。

int w, h;
int state0;
int i=0;
double *x_output, *y_output, *z_output, *vx_output, *vy_output, *vz_output;
int pointnum;
double point[MAXPOINT][3];
double distance, twist, elevation, azimuth;

//関数(モジュール)の準備/////////////////////////////////////////////
double Bx(double x, double y, double z){
  return 0.0;
}

double By(double x, double y, double z){
  return 0.0;
}

double Bz(double x, double y, double z){
  return -3.0;
}

double func_dx_dt(double x, double y, double z,
		  double vx, double vy, double vz){
  return vx;
}

double func_dy_dt(double x, double y, double z, 
		  double vx, double vy, double vz){
  return vy;
}

double func_dz_dt(double x, double y, double z, 
		  double vx, double vy, double vz){
  return vz;
}

double func_dvx_dt(double x, double y, double z, 
		   double vx, double vy, double vz){
  return vy*Bz(x, y, z) - vz*By(x, y, z);
}

double func_dvy_dt(double x, double y, double z,
		   double vx, double vy, double vz){
  return vz*Bx(x, y, z) - vx*Bz(x, y, z);
}

double func_dvz_dt(double x, double y, double z,
		   double vx, double vy, double vz){
  return vx*By(x, y, z) - vy*Bx(x, y, z);
}

void RK4(double x, double y, double z, double vx, double vy, double vz){
  double kx[4], ky[4], kz[4], kvx[4], kvy[4], kvz[4];
  kx[0]  = dt * func_dx_dt(x, y, z, vx, vy, vz);
  ky[0]  = dt * func_dy_dt(x, y, z, vx, vy, vz);
  kz[0]  = dt * func_dz_dt(x, y, z, vx, vy, vz);
  kvx[0] = dt * func_dvx_dt(x, y, z, vx, vy, vz);
  kvy[0] = dt * func_dvy_dt(x, y, z, vx, vy, vz);
  kvz[0] = dt * func_dvz_dt(x, y, z, vx, vy, vz);

  kx[1]  = dt * func_dx_dt(x  + kx[0]/2.0,  y  + ky[0]/2.0,  z  + kz[0]/2.0, 
			   vx + kvx[0]/2.0, vy + kvy[0]/2.0, vz + kvz[0]/2.0);
  ky[1]  = dt * func_dy_dt(x  + kx[0]/2.0,  y  + ky[0]/2.0,  z  + kz[0]/2.0, 
			   vx + kvx[0]/2.0, vy + kvy[0]/2.0, vz + kvz[0]/2.0);
  kz[1]  = dt * func_dz_dt(x  + kx[0]/2.0,  y  + ky[0]/2.0,  z  + kz[0]/2.0, 
			   vx + kvx[0]/2.0, vy + kvy[0]/2.0, vz + kvz[0]/2.0);
  kvx[1] = dt * func_dvx_dt(x + kx[0]/2.0,  y  + ky[0]/2.0,  z  + kz[0]/2.0, 
			   vx + kvx[0]/2.0, vy + kvy[0]/2.0, vz + kvz[0]/2.0);
  kvy[1] = dt * func_dvy_dt(x + kx[0]/2.0,  y  + ky[0]/2.0,  z  + kz[0]/2.0, 
			   vx + kvx[0]/2.0, vy + kvy[0]/2.0, vz + kvz[0]/2.0);
  kvz[1] = dt * func_dvz_dt(x + kx[0]/2.0,  y  + ky[0]/2.0,  z  + kz[0]/2.0, 
			   vx + kvx[0]/2.0, vy + kvy[0]/2.0, vz + kvz[0]/2.0);

  kx[2]  = dt * func_dx_dt(x  + kx[1]/2.0,  y  + ky[1]/2.0,  z  + kz[1]/2.0, 
			   vx + kvx[1]/2.0, vy + kvy[1]/2.0, vz + kvz[1]/2.0);
  ky[2]  = dt * func_dy_dt(x  + kx[1]/2.0,  y  + ky[1]/2.0,  z  + kz[1]/2.0, 
			   vx + kvx[1]/2.0, vy + kvy[1]/2.0, vz + kvz[1]/2.0);
  kz[2]  = dt * func_dz_dt(x  + kx[1]/2.0,  y  + ky[1]/2.0,  z  + kz[1]/2.0, 
			   vx + kvx[1]/2.0, vy + kvy[1]/2.0, vz + kvz[1]/2.0);
  kvx[2] = dt * func_dvx_dt(x + kx[1]/2.0,  y  + ky[1]/2.0,  z  +  kz[1]/2.0, 
			   vx + kvx[1]/2.0, vy + kvy[1]/2.0, vz + kvz[1]/2.0);
  kvy[2] = dt * func_dvy_dt(x + kx[1]/2.0,  y  + ky[1]/2.0,  z  + kz[1]/2.0, 
			   vx + kvx[1]/2.0, vy + kvy[1]/2.0, vz + kvz[1]/2.0);
  kvz[2] = dt * func_dvz_dt(x + kx[1]/2.0,  y  + ky[1]/2.0,  z  + kz[1]/2.0, 
			   vx + kvx[1]/2.0, vy + kvy[1]/2.0, vz + kvz[1]/2.0);

  kx[3]  = dt * func_dx_dt(x  + kx[2],  y  + ky[2],  z  + kz[2], 
			   vx + kvx[2], vy + kvy[2], vz + kvz[2]);
  ky[3]  = dt * func_dy_dt(x  + kx[2],  y  + ky[2],  z  + kz[2], 
			   vx + kvx[2], vy + kvy[2], vz + kvz[2]);
  kz[3]  = dt * func_dz_dt(x  + kx[2],  y  + ky[2],  z  + kz[2], 
			   vx + kvx[2], vy + kvy[2], vz + kvz[2]);
  kvx[3] = dt * func_dvx_dt(x + kx[2],  y  + ky[2],  z  +  kz[2], 
			   vx + kvx[2], vy + kvy[2], vz + kvz[2]);
  kvy[3] = dt * func_dvy_dt(x + kx[2],  y  + ky[2],  z  + kz[2], 
			   vx + kvx[2], vy + kvy[2], vz + kvz[2]);
  kvz[3] = dt * func_dvz_dt(x + kx[2],  y  + ky[2],  z  + kz[2], 
			   vx + kvx[2], vy + kvy[2], vz + kvz[2]);

  x = x + kx[0]/6.0 + kx[1]/3.0 + kx[2]/3.0 + kx[3]/6.0;
  y = y + ky[0]/6.0 + ky[1]/3.0 + ky[2]/3.0 + ky[3]/6.0;
  z = z + kz[0]/6.0 + kz[1]/3.0 + kz[2]/3.0 + kz[3]/6.0;
  vx = vx + kvx[0]/6.0 + kvx[1]/3.0 + kvx[2]/3.0 + kvx[3]/6.0;
  vy = vy + kvy[0]/6.0 + kvy[1]/3.0 + kvy[2]/3.0 + kvy[3]/6.0;
  vz = vz + kvz[0]/6.0 + kvz[1]/3.0 + kvz[2]/3.0 + kvz[3]/6.0;

  x_output = &x;
  y_output = &y;
  z_output = &z;
  vx_output = &vx;
  vy_output = &vy;
  vz_output = &vz;
}


//文字列描画モジュール///////////////////////////////////////////////////////
static void DrawString(char *str, int w, int h, int x0, int y0){
  int k;
  glDisable(GL_LIGHTING);
  // 平行投影にする
  glMatrixMode(GL_PROJECTION);
  glPushMatrix();
  glLoadIdentity();
  gluOrtho2D(0, w, h, 0);
  glMatrixMode(GL_MODELVIEW);
  glPushMatrix();
  glLoadIdentity();
  
  // 画面上にテキスト描画
  glRasterPos2f(x0, y0);
  int size = (int)strlen(str);
  for(k = 0; k < size; ++k){
    char ic = str[k];
    glutBitmapCharacter(GLUT_BITMAP_9_BY_15, ic);
  }
  
  glPopMatrix();
  glMatrixMode(GL_PROJECTION);
  glPopMatrix();
  glMatrixMode(GL_MODELVIEW);
}


//極ビューイング変換///////////////////////////////////////////////////////////
void polarview(double disatance, double twist, 
	       double elevation, double azimuth){
  glTranslated(0.0, 0.0, -distance);
  glRotated(twist, 0.0, 0.0, 1.0);
  glRotated(elevation, 1.0, 0.0, 0.0);
  glRotated(azimuth, 0.0, 1.0, 0.0);
}


//画面のリサイズ////////////////////////////////////////////////////////////////
void resize(int width, int height){
  glViewport(0,0,width,height); //ウィンドウ全体をビューポートにする
  glMatrixMode(GL_PROJECTION);  //透視変換行列の設定
  glLoadIdentity();             //透視変換行列の初期化
  gluPerspective(17.0, (double)width / (double)height, 1.0, 100.0);//透視変換(画角、縦横比、手前、奥)
  glMatrixMode(GL_MODELVIEW);   //モデルビュー変換行列の設定
  w = width;
  h = height;
}

//マウス////////////////////////////////////////////////////////////////////

void mouse(int button, int state, int x, int y){

    if(button == GLUT_LEFT_BUTTON && state == GLUT_DOWN){
      state0 = 2;
    }else if(button == GLUT_LEFT_BUTTON && state == GLUT_UP){ 
      state0 = 1;
    }else{

    }
}

//地面(磁力線)を書く/////////////////////////////////////////////////////////////

void mag_line(void){

  glColor3d(0.063, 0.306, 0.545) ;  //磁力線の色を設定
  glLineWidth(1.0);                 //磁力線の太さを設定
  glBegin(GL_LINES);                //以下glEndまでに指定された点を線で結ぶ
  glVertex3d(0.0, 0.0,  iz);
  glVertex3d(0.0, 0.0, -iz);
  glEnd();

}


//軌跡をなぞる関数///////////////////////////////////////////////////////////
void kiseki(double point[][3], int pointnum){
  int j;
  glEnable(GL_LINE_STIPPLE);   //glLineStippleを有効にする
  glLineStipple(1 , 0xF0F0);   //0xF0F0(破線)という線の繰り返しパターンを設定
  if(pointnum > 1){
    glColor3d(0.0, 0.0, 0.0);  //線の色を黒に設定した
    glLineWidth(1.0);          //線の幅(太さ)を1.0に設定した
    glBegin(GL_LINE_STRIP);    //以下で指定する(pointnum)個の点を線で結ぶ
    for( j=0 ; j<pointnum ; ++j){
      glVertex3dv(point[j]);   //点をdouble型の3次元ベクトル(配列)で指定する
    }
    glEnd();
  }
  glDisable(GL_LINE_STIPPLE);  //glLineStipple(破線)を使うのはここまで(こまめに電源を切るのと同じ)
}


//アイドリング時の関数/////////////////////////////////////////////////////////
void idle(void){
  glutPostRedisplay();   //glutDisplayFuncに登録したコールバック関数を呼び出す
}


//画面表示///////////////////////////////////////////////////////////////////
void display(void){
  static GLfloat lightpos[4] = {0.0, 5.0, 0.0, 0.0};  //光源の位置を決定
  static GLfloat lightcol[4] = {0.0, 0.6, 1.0, 0.0};  //光源の色を決定
  static GLfloat white[4] = {0.8, 0.8, 0.8, 1.0};     //球の色を決定
  static double vx, vy, vz;
  static double x,  y,  z;
  
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //画面初期化

  /*シーンの描画*/
  glMatrixMode(GL_MODELVIEW); //オブジェクト座標系をワールド座標系の原点に配置(Obj-1)
  glPushMatrix();             //Obj-1を保存してそこから派生したObj-1_1を作成
  polarview(distance, twist, elevation, azimuth);  //Obj-1_1がワールド座標系の中で向きを変える
  glPushMatrix();             //Obj-1_1を保存してそこから派生したObj-1_1_1を作成
  glTranslated(x, y, z);      //Obj-1_1_1がワールド座標系の中で動く
  glMaterialfv(GL_FRONT, GL_DIFFUSE, white);   //これからつくる物体の質感を設定
                                               //GL_FRONTは表側
                                               //GL_DIFFUSEは拡散反射係数
                                               //whiteは球の色(4次元color)
  glLightfv(GL_LIGHT0, GL_POSITION, lightpos); //ライトの位置を決定
  glLightfv(GL_LIGHT0, GL_DIFFUSE, lightcol);  //ライト(拡散光)の色を決定
  glEnable(GL_LIGHTING);     //全ライト有効
  glEnable(GL_LIGHT0);       //ライト0電源ON
  glutSolidSphere(0.07, 16, 8); //球描画(半径, Z軸まわり分割数, Z軸方向分割数)、Obj-1_1_1の原点に描画する。
  glDisable(GL_LIGHT0);      //ライト0電源OFF
  glDisable(GL_LIGHTING);    //全ライト無効
  glPopMatrix();             //Obj-1_1_1からObj-1_1に戻る

  /*位置をRK4で動かす*/
  if(state0==0){
    x = ix, y = iy, z = iz;  //位置初期に戻す
      vx = 0.0, vy = 0.0, vz = 0.0;
    i = 0;
  }else if(state0==1){
    if(x==ix && y==iy && z==iz && vx==0.0 && vy==0.0 && vz==0.0){
      vx = ivx, vy = ivy, vz = ivz; //初期速度を与えてやる
    }else if(z < -(iz)){
      x = ix, y = iy, z = iz;
      vx = ivx, vy = ivy, vz = ivz;
      i = 0;
      pointnum = 0;
    }else{
      RK4(x, y, z, vx, vy, vz); //入れた○が計算されて○_outputの形で出てくる
      vx = *vx_output;
      vy = *vy_output;
      vz = *vz_output;
      x  = *x_output;
      y  = *y_output;
      z  = *z_output;
      point[i][0] = x;
      point[i][1] = y;
      point[i][2] = z;
      i = i+1;
      pointnum = pointnum + 1;
      kiseki(point, pointnum);
    }
  }else if(state0==2){
    kiseki(point, pointnum);
  }else{

  }
  mag_line();     //Obj-1_1に磁力線を書いた
  glPopMatrix();  //Obj-1_1からObj-1(初期設定)へ戻った
  
  char str1[64];
  sprintf(str1, "distance:%.1lf", distance);  //変数である数字を文字列に変換しておく
  glColor3d(0.0, 0.0, 0.0);                   //文字色を黒に設定
  DrawString("enter_key : start", w, h, 10, 10);  
  DrawString("up_key    : move up", w, h, 10, 25);
  DrawString("down_key  : move down", w, h, 10, 40);
  DrawString("left_key  : move left", w, h, 10, 55);
  DrawString("right_key : move right", w, h, 10, 70);
  DrawString("z : zoom in", w, h, 10, 85);
  DrawString("x : zoom out", w, h, 10, 100);
  DrawString("r : restart", w, h, 10, 115);
  DrawString("q : exit", w, h, 10, 130);
  DrawString(str1, w, h, 10, 145);


  //glFlush();  //これまでは実は計画段階だった。ここで初めて上の計画通りに一気に描画。
  glutSwapBuffers();  //glFlushに加えて、2枚のバッファを入れ替える機能を追加したもの

}



//初期設定関数////////////////////////////////////////////////////////////////
void init(void){
  glClearColor(1.0,1.0,1.0,1.0); //ウィンドウカラー(赤,緑,青,透明度α値)0〜1で指定
  glEnable(GL_DEPTH_TEST);       //視点方向に物体が重なったときの距離判定を有効にする
  distance = 15.0;               //パラメータ設定
  twist = 0.0;                   //パラメータ設定
}


//キーボード//////////////////////////////////////////////////////////////////
void keyboard(unsigned char key, int x_screen, int y_screen){
  switch(key){
  case 'q':
  case 'Q':
  case '\033':
    exit(0);
  case 'r':
    state0 = 0;
    pointnum = 0;
    //x = ix, y = iy, z = iz; //球の位置初期設定
    break;
  case 'z':
    if(distance>1.0){
      distance = distance - 1.0;
    }else{ 

    }
    break;
  case 'x':
    distance = distance + 1.0;
    break;
  case '\015':
    state0 = 1;
    break;
  default:
    break;
  }
}


//キーボード(特殊)/////////////////////////////////////////////////////////////
void special(int key, int x, int y){
    if(key == GLUT_KEY_UP){
      elevation = elevation + 2.0;
    }else if(key == GLUT_KEY_DOWN){
      elevation = elevation - 2.0;
    }else if(key == GLUT_KEY_LEFT){
      azimuth = azimuth + 2.0;
    }else if(key == GLUT_KEY_RIGHT){
      azimuth = azimuth - 2.0;
    }else{

    }
}



//メイン関数/////////////////////////////////////////////////////////////////
int main(int argc, char *argv[]){
  glutInit(&argc, argv);         //glutライブラリの初期化
                                 //プラットフォームのウィンドウシステムが初期化される
                                 //argcはコマンドラインで付加したオプション
  glutInitWindowPosition(0,0);   //ウィンドウの初期位置を決定する
  glutInitWindowSize(1500,800);  //ウィンドウの初期サイズを決定する
  glutInitDisplayMode(GLUT_RGBA | GLUT_DEPTH | GLUT_DOUBLE); 
                                 //displayの表示モードを設定する
                                 //GLUT_RGBAは4次元色(Red,Green,Blue,α)を使えるバッファを用意
                                 //GLUT_DEPTHは視線方向に重なった物体の深さ情報を記憶できるバッファを用意
                                 //GLUT_DOUBLEはバッファを2枚用意して交互に表示するダブルバッファ(ちらつきを防ぐ)
  glutCreateWindow(argv[0]);     //ウィンドウを開く(name)

  glutDisplayFunc(display);      //ウィンドウ内に絵を描く関数を決める
                                 //glutCreateWindow〜glutMainLoopの間に置く。
                                 //必要なことはdisplay関数に書くことにしてます。
  glutReshapeFunc(resize);       //ウィンドウに変形があったときに使うリサイズ関数を指定
  glutKeyboardFunc(keyboard);    //キーが押されたときに使う関数を指定
  glutSpecialFunc(special);      //特殊キーが押されたときに使う関数を指定
  glutMouseFunc(mouse);          //マウス操作があったときに使う関数を指定
  glutIdleFunc(idle);            //何も操作がないときに実行される関数を指定
  init();                        //上で定めた初期化関数
  /*↑ココまでの関数は全て一回ずつ実行される*/

  glutMainLoop();                //イベント発生待機状態(ある時間間隔ごとに実行)
                                 //その都度、行われたイベントを感知し、そのイベントに対応した関数を呼び出す
                                 //この関数は終了まで延々と繰り返し実行される
  return 0;
}

//関数displayが機能するのはglutMainLoop()が終了するとき。
//なお、glutReshapeFuncで指定したresize関数はウィンドウを開いた瞬間に1回呼び出される