// dirty_main.cpp : Numerical Recipes の rk4.c を使って
// ２次振動系の微分方程式を解く（「杉戸」問題の正しいアプローチ）
// 2002/04/23 by SK

// 注意：Numerical Recipesのソースを１行も書き直さないで，
//       C のソースをC++プログラムとして使い倒し，
//       さらにディフォルトの浮動小数 float を，私の好きな double として使う。
//       このため，以下では汚い技を２つ使うからね！

// MS-Visual C++ 使用者へ：注意その１
// メニューの ツール＞オプション＞ディレクトリのところで
// 表示するディレクトリ＝「インクルードファイル」の欄に
// 次の２つのパスを書き込んでおくこと。
// ヘッダファイル（include）の在り処：
//     E:\SKON\C\NUMRECIPI\INCLUDE（近藤の場合は）
// C のソースと nrutil.c プログラムの在り処：
//     E:\SKON\C\NUMRECIPI\SRC（近藤の場合は）
// 上の２つ，どちらも「インクルードファイル」の欄だからね！
// これで，Visual C++ が勝手に見つけてくれる。

// float を double に読み替える。（あまりにも強引，きったねー！その１）
#define float double 

// 次の３行は，ほとんどいつも必要！！
#include "nr.h"			
#include "nrutil.h"
#include "nrutil.c"

// 使いたい Numerical Recipes の C ソースを取り込む！
#include "rk4.c"
// このファイル dirty_main.cpp の拡張子が .cpp だから，
// 取り込んでしまえばコンパイラには C++ プログラムに見える。
// だまし討ち。（あまりにも強引，きったねー！その２）

// MS-Visual C++ 使用者へ：注意その２
// ワークスペースの外部依存関係の中に
//     nrutil.c や rk4.c 
// が見つかるが，Source Files へ絶対に移動してはいけない！！！
// （もしそうすると，二重定義でコンパイラに叱られる。） 
// 一方，***.h ファイルは，Header Files へ移動してもいいけど，
// 外部依存関係の中に入れておいたほうがよろしい！（Visual C++ が面倒みてくれるから。）

// ******* 以上で Numerical Recipes を使う準備ができた！！ ********


//  ここから，杉戸問題を解こう！　つまり，２次の微分方程式を解く
//  cos, sin の振動をシミュレーションし，ルンゲ・クッタの精度を確認する。

#define PI 3.14159265358979323846

// cos, sin 振動する系の微分方程式の右辺を準備する。周期 = 1.0 sec
void func_sin(double t, double y[], double dydt[])
{
	dydt[1] = - 2.0 * PI * y[2];  // dy1/dt = d/dt(cos) = -omega*sin, frequency = 1 Hz 
	dydt[2] =   2.0 * PI * y[1];  // dy2/dt = d/dt(sin) = +omega*cos,
}

#include <stdio.h>
#include <stdlib.h>	// exit()
#include <math.h>	// sin(), cos() <== 計算誤差を評価するため！

int main(void)
{
	double t=0.0, h=0.01/1.0, tstart=0.0, towari=5.0;
	FILE *d_2ndA2;
	double *x, *dxdt;
	int Norder = 2;

	x = vector(1,Norder);       // rk4() の引数にあわせて，２次元ベクトルを確保する。
	dxdt = vector(1,Norder);
	x[1] = 1.0;		x[2] = 0.0; //  x[1] = cos(0), x[2]=sin(0) 状態変数の初期値

	if( (d_2ndA2 = fopen("d_2ndA2.dat","w")) == NULL) {
		printf("\nerror: ファイル d_2ndA2.dat が開けません\n");
		exit(110);
	}
	printf("\ntime = %f x1 = %f x2 = %f :err %g %g", t, x[1], x[2], x[1]-cos(2.0*PI*t), x[2]-sin(2.0*PI*t) );
	fprintf(d_2ndA2, "t x1 x2 x1err x2err \n");
	fprintf(d_2ndA2, "%f %f %f %g %g \n", t, x[1], x[2], x[1]-cos(2.0*PI*t), x[2]-sin(2.0*PI*t) );

	while (t < towari) {
		// Numerical Recipes の rk4() は次の２行をペアにして使うこと。
		func_sin(t, x, dxdt);		// rk4 の第１ステップの準備
		rk4(x, dxdt, 2, t, h, x, func_sin);
		// 上の２行はペア（可変刻みの時に備えてこうなっている。）

		t = t + h;
		printf("\ntime = %f x1 = %f x2 = %f :err %g %g", t, x[1], x[2], x[1]-cos(2.0*PI*t), x[2]-sin(2.0*PI*t) );
		fprintf(d_2ndA2, "%f %f %f %g %g \n", t, x[1], x[2], x[1]-cos(2.0*PI*t), x[2]-sin(2.0*PI*t) );
	}
	printf("\n");
	fclose(d_2ndA2);

	free_vector(dxdt, 1, 2);     // vectorを使ったら，後片付けすること！！！！
	free_vector(x, 1, 2);
	return 0;					// おわり！
}

// 刻み幅h と 計算誤差の関係！　違うプログラムでやったので異なるかもしれない。
// 刻み幅 h		t=10[s] での誤差（x1, x2 の大きいほう）		有効桁数
//	0.01*20		-0.742036		振幅=0.257964 使用不可		0
//	0.01*10		-0.0675592		振幅=0.957799 減衰大きい	1.5
//	0.01*8		-0.0300467		振幅=0.986073 減衰			2
//	0.01*4		-0.00204127		振幅=0.999564				3.5
//	0.01*2		-0.000129831	振幅=0.999986				4
//	0.01		-8.14902e-6		振幅=1.000000				5
//  0.01/2		-5.09854e-7		振幅=						6.5
//	0.01/4		-3.18752E-8		振幅=						7
//	0.01/8		-1.99501E-9		振幅=						8
//	0.01/16		-1.24307E-10	振幅=						9
//	0.01/32		2.08992e-11		振幅=						10
// end of file