// 一機無限大母線系統 + AVR のシミュレーション

#include "trs5-inc.mm"
#include "tincas2_inc.mm"

void dat_output();
void rem_dat_output_16();


Func void main() {
     Matrix x;
     Real t;
     Integer co;
     Integer USE_GRAPHICS;

     USE_GRAPHICS= yes;
//     USE_GRAPHICS= no;

     h=0.0;  t=0.0;  t_fin=0.0;

     co= 0;  	// auto clear の為のカウンタ
     update_interval= 100; // tincas canvas の更新時間間隔
     auto_clear= 1; // 定期的に canvas をクリアする。

     x = Z(5,1);		// 状態変数は最大５次元
//   各ベクトル、行列の領域確保
     eps=   Z(x);
     xbase= Z(x);
     x_ini= Z(x);

     set_parameters();		// パラメータのデフォルト値を設定

     x= x_ini;			// 初期状態を x_ini に設定

//   計算する次元 comp_dimにしたがって resize
     x =         x(1:comp_dim, 1);
     eps =     eps(1:comp_dim, 1);
     xbase = xbase(1:comp_dim, 1);

     /* *************** パラメータ値の仮変更はここで ****************** */
     TdoD1= 12.0;  TV1= 0.69;
//     UlV1= 1.0e5;  LlV1= -1.0e5;  GV1= 10.0;  Pms1= 1.6; TV1= 0.8;

     // 初期値をファイルから(存在すれば)読み込む
     if (access("./xini.mat")==0){ read x << "./xini.mat"; }
	   
     t_fin= 100.0;

     if (USE_GRAPHICS==yes){
	  pid = open_canvas([ 300, 300 ]);
	  set_canvas_xyrange([0.0,120.0, -0.01, 0.01]);//xmin,xmax, ymin,ymax
	  set_tincas_title(pid, "trs5");
	  set_tincas_update_interval(pid, update_interval);
     }else{
//	  printf("# %d %d\n",length(x)+1,Integer(t_fin/h)+2); //MAT data header
	  //// set_ini_to_equlibrium(f, jm, x, eps, xbase);
	  printf("# timestep= %.15f\n",h);	// 数値積分の時間刻みを出力
	  rem_dat_output_16(x, "#initial_state");  // 初期状態を出力 (１６桁)
	  ftrs5(t,x);  dat_output(x);		// 初期状態を出力
	  rem_prm_output();			// globalパラメータの値を出す
     }
//   数値積分のループ -----------------------
     while ( ( t < t_fin )||( USE_GRAPHICS==yes ) ){
	  x = rngkut4(t,x,ftrs5,h);  t = t + h;
//	  {t, x} = rk4(ftrs5,t,x,h);		//  １ステップ進める。
	  if (USE_GRAPHICS==yes){
	       pset( pid, gxy([x(1)/PI*180  x(2)]') );
	       {t, x}= kbd_event(t, x);
	       if (auto_clear && (co++)>auto_clear_interval){
		    clear_canvas(pid); co= 0;
	       }
	  }else{
	       dat_output(x);			// 現在の状態を出力
	  }
     }
//   ----------------------------------------
     if (USE_GRAPHICS==yes){
	  close_canvas(pid);
     }else{
	  rem_dat_output_16(x, "#final_state" );  // 最終値を出力 (１６桁)
     }
}


Func void dat_output(x)
     Matrix x;
{
     switch(systype){
     case FirstOrderAVR:
	  printf("%g %g %g %g\n",x(Ndlt1)/PI*180, x(Nw1), x(NeqD1), x(NgV1));
	  break;
     case WithVMdelay:
     default:
	  printf("%g %g\n",Vt1, Pe1);
	  break;
     }
}


Func void rem_dat_output_16(x, s)
     Matrix x;
     String s;
{
     switch(systype){
     case FirstOrderAVR:
	  printf("#%.16g %.16g %.16g %.16g %s",x(1)/PI*180, x(2), x(3), x(4), s);
	  break;
     case WithVMdelay:
	  printf("#%.16g %.16g %.16g %.16g %.16g %s",x(1), x(2), x(3), x(4), x(5), s);
     default:
	  break;
     }
     printf("\n");
}