// 一機無限大母線系統 + 仮想発電機 のシミュレーション

#include "trs5-inc.mm"

void dat_output();
void rem_dat_output_16();
List bifudiag_output();
Real outvf(), secf();
Integer is_identical();

Func void main() {
     Matrix x;
     Real t;
     Real t_fin_bifu;
     Matrix tmpx;

     Integer USE_GRAPHICS;
     Integer conv, p, breakflg;
     Integer fd;

     settimer();
     USE_GRAPHICS= no;
     breakflg= 0;

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

     x = Z(5,1);		// 状態変数は5次元

//   各ベクトル、行列の領域確保
     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);

     printf("# timestep= %.15f\n",h);		// 数値積分の時間刻みを出力
     //// set_ini_to_equlibrium(f, jm, x, eps, xbase);
     rem_dat_output_16(x, "#initial_state");	// 初期状態を出力 (１６桁)

     fd= 1; // 既に開いている標準出力を用いる。

//   数値積分のループ -----------------------
     TV1= 0.9;
     t_fin_bifu= 1500.0;
     printf("# free parameter is GV1.\n");
     ftrs5(t,x);  rem_prm_output();		// globalパラメータの値を出す
     for (GV1 = 10; GV1 <= 300.01; GV1 = GV1+0.5){
	  x(Nw1) = x(Nw1)+1.0e-7; // 分岐を促進
	  for(t=0.0; ;){
	       x = rngkut4(t,x,ftrs5,h);  t = t + h;
	       // {t, x} = rk4(ftrs5,t,x,h);		//  １ステップ進める。
	       tmpx = x; // xが間違って変更されないように xtmp に複製して使う
	       // 定常状態に収束 または 既定の時間に達した場合は 次へ。
	       {conv, p}= bifudiag_output(GV1, secf(x), outvf(x), t, t_fin_bifu, is_identical, tmpx, fd);
	       if (conv==1){ break; }
	       if ( x(Ndlt1)>1000.0 ){ breakflg= 1;  break; } // 脱調
	  }
	  if (breakflg){ break; }
     }
//   ----------------------------------------
     printf("# CPU time= %g(%gh)\n",gettimer(),gettimer()/3600.0);
     rem_dat_output_16(x, "#final_state" );	// 最終値を出力 (１６桁)
}


// secf(x)=0 がポアンカレ断面となる
Func Real secf(x)
     Matrix x;
{
     return x(Nw1);
}

// 分岐図の縦軸になる値
Func Real outvf(x)
     Matrix x;
{
     return x(Ndlt1)*180.0/PI;
}



// 同じ点とみなしていいかどうかの判定
Func Integer is_identical(x1, x2, mul)
     Matrix x1, x2;
     Real mul; // 普通は1.0。甘くしたいときに大きくする。
{
     Real distance, threshold;

     // UlV によるアトラクタの伸縮を考慮
     threshold= mul * 1.0e-4*(UlV1-1.90929);

     distance= norm( Array(x2-x1)/Array(xbase) );

     if ( distance < threshold ){
//	  printf("- %g %g\n",distance, threshold);
	  return 1;
     }else{
	  return 0;
     }
}



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");
}