//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
//  Program example.cc
//  Fractal分解を利用した Symplectic integrator による
//  物理計算のプログラム例
//  Copyright (C) by 渡辺尚貴 1998年 8月14日
//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=

#include <stdio.h>
#include <math.h>
#include "vector2.h"
#include "nxgraph.h"

// fracdeco.cc の出力結果をこのファイルにでも入れておいて下さい。
#include "symplectic.h"

//---- 諸定数の設定(好み合わせて変更して下さい)
#define dT             (1.0/1024)  // 基本 時間刻幅 
#define FINAL_T        (2048.0)    // 最終時刻

#define POINTS_PER_DOT (16)        // 縦1dot当たりのプロット数

// プロットする際のエネルギー軸の拡大率(かなり大きな値が必要です)
#define MAG            (128*256.0) // 2次の分割向け
//#define MAG            (128*128*256.0) // 4次の分割向け
//#define MAG            (128*128*128*128*128.0) // 6次の分割向け
//#define MAG            (128*128*128*128*256.0) // 8次の分割向け


//---- 諸定数の定義(手を触れないで下さい)
#define SCALE_WIDTH    (512)       // 横軸の長さ(dot)
#define SCALE_HEIGHT   (256)       // 縦軸の長さ(dot)
#define XMARGIN        (40)        // 画面の枠と軸との隙間の横幅(dot)
#define YMARGIN        (16)        // 画面の枠と軸との隙間の縦幅(dot)

#define WIN_WIDTH    (SCALE_WIDTH +2*XMARGIN)
#define WIN_HEIGHT   (SCALE_HEIGHT+2*YMARGIN)
#define XO           (XMARGIN)
#define YO           (YMARGIN+SCALE_HEIGHT/2)
#define STEP_PER_POINT (int(FINAL_T/SCALE_WIDTH/POINTS_PER_DOT/dT))


//---- 座標軸や目盛などを描く関数
void MakeSheet( void )
{
  double ymax, yrange, y_value;
  int y_height;

    // 縦軸と横軸を描く
  NXDrawLine( XO, YO-SCALE_HEIGHT/2, XO, YO+SCALE_HEIGHT/2 );
  NXDrawLine( XO, YO+SCALE_HEIGHT/2, XO+SCALE_WIDTH, YO+SCALE_HEIGHT/2 );

    // 拡大率(MAG)に合わせて縦軸の目盛を選んで描く
  ymax   = (SCALE_HEIGHT/2)/MAG;
  yrange = pow( 10, (int)floor( log10( ymax ) ) );

  if( 5*yrange < ymax ){
    y_value  = 5*yrange;
  }else if( 2*yrange < ymax ){
    y_value  = 2*yrange;
  }else{
    y_value  = yrange;
  }
  y_height = int(y_value*MAG);

    // 縦軸の＋側に目盛を入れる
  NXDrawMoveto( XO-3, YO-y_height );
  NXDrawLinerel( 6, 0 );
  NXDrawPrintf( 0, YO-y_height+4,"%+5.0e", +y_value );

    // 縦軸の原点に目盛を入れる
  NXDrawPrintf( 0, YO+4, "%0" );

    // 縦軸のー側に目盛を入れる
  NXDrawMoveto( XO-3, YO+y_height );
  NXDrawLinerel( 6, 0 );
  NXDrawPrintf( 0, YO+y_height+4,"%+5.0e", -y_value );

    // 縦軸の変数名を描く
  NXDrawPrintf( XO-16, YO-SCALE_HEIGHT/2, "Energy" );

    // 横軸の原点に目盛を入れる
  NXDrawPrintf( XO-6, YO+SCALE_HEIGHT/2+16, "0" );

    // 横軸の中間に目盛を入れる
  NXDrawMoveto( XO+SCALE_WIDTH/2, YO+SCALE_HEIGHT/2-3 );
  NXDrawLinerel( 0, 6 );
  NXDrawPrintf( XO+SCALE_WIDTH/2-12, YO+SCALE_HEIGHT/2+16, "1000" );

    // 横軸の右端に目盛を入れる
  NXDrawMoveto( XO+SCALE_WIDTH, YO+SCALE_HEIGHT/2-3 );
  NXDrawLinerel( 0, 6 );
  NXDrawPrintf( XO+SCALE_WIDTH-12, YO+SCALE_HEIGHT/2+16, "2000" );

    // 横軸の変数名を描く
  NXDrawPrintf( XO+SCALE_WIDTH-16, YO+SCALE_HEIGHT/2-8, "Time" );

    // Symplectic法の種類を表す数値を描く
  NXDrawPrintf( WIN_WIDTH/2, YMARGIN, "m=%d  r=%d", symp_m, symp_r );

  NXFlush();
}

//---- 時間を横軸、エネルギーを縦軸にプロットする関数
void Plot( double t, double E )
{
  NXDrawPoint( XO+int(SCALE_WIDTH*t/FINAL_T), YO - int(MAG*E) );
}

//---- 質点に働く力を計算する関数
Vector2 Force( const Vector2& q )
{
  return( Vector2(-q.x*q.y*q.y, -q.x*q.x*q.y ) );
}

//---- 質点のエネルギーを計算する関数
double Energy( const Vector2& q, const Vector2& p )
{
  return( 0.5*(p.x*p.x + p.y*p.y + q.x*q.x*q.y*q.y) - 2.0 );
}


int main( void )
{
    // 質点の位置 q と運動量 p の変数の定義と初期化
  Vector2 q(2,1), p(0,0);

  NXOpenWindow( "example of symplectic integrator", WIN_WIDTH, WIN_HEIGHT );

    // 座標軸などの描画
  MakeSheet();

    // 時間発展のループ
  for( double t=0.0 ; t<FINAL_T ; ){
      // いくらかstepを踏んでからプロットするためのループ
    for( int step=0 ; step<STEP_PER_POINT ; step++ ){
        // Symplectic integratorの心臓部
      for( int r=0 ; r<symp_r ; r++ ){
	Vector2 f =  Force( q );   // 力の計算
	p += dT * symp_u[r] * f;   // Duの効果による運動量の更新
	q += dT * symp_k[r] * p;   // Dkの効果による位置の更新
      }
      t+=dT;
    }

      // エネルギーの時間変化を描く
    Plot( t, Energy( q, p ) );
  }

  printf("プログラムは終了です。クリックしてください。\n");
  XEvent ev;
  NXCheckEvent( NX_WAIT, ev );
  NXCloseWindow();
}