#include <stdio.h>

#include "nrutil.h"

#define CNT 1000
#define VarNum	2 // 変数の個数

char *filename1 = "yout.dat";

FILE *fpdata1;

void func(double time, double y[], double dydt[]) // 状態空間表現式を入れる
{
	static double	damper = 1.00;
	static double	spring = 4.00;
	
	// 好きな任意の状態空間表現式に入れ替える
	dydt[0] = y[1];
	dydt[1] = - damper * y[1] - spring * y[0];
}

int main()
{
	int i;
	double time = 0.0;
	double kizami = 0.01;
	double y[VarNum],dydt[VarNum];

	if ((fpdata1 = fopen(filename1,"w")) == NULL){
		fprintf(stderr,"can not open %s\n",filename1);
		return 0;
	}

	// 出力の初期値設定
	y[0] = 1.0;
	y[1] = 0.0;

	for (i=0;i<CNT;i++){
		rk4(y,dydt,VarNum,time,kizami,y,func); // ルンゲクッタ実行
		fprintf(fpdata1,"%lf %lf %lf\n",time,y[0],y[1]);
	
		time += kizami; // 時刻のインクリメント
	}

	fclose(fpdata1);
	return 0;
}