// 2-parameter 分岐図

#include "trs5-inc.mm"

#define CHOTTO 0.0000001

void dat_output();
void rem_dat_output_16();
List bifudiag_output();
Real outvf(), secf();
Integer is_identical();
Real distp();
void makepsfile();
List record_2p_bifurcation_point();
void output_2p_bifurcation_point_to_matfiles();

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

     Integer USE_GRAPHICS;
     Integer conv, p;
     Real step1, step1base;
     Real step2, step2base;
     Integer fd;
     String fname;

     settimer();
     USE_GRAPHICS= no;

     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();		// パラメータのデフォルト値を設定

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

     // ********** 一時的なパラメータ値変更はここで。*****************
//     UlV1 = 1.0e5;  LlV1 = -1.0e5;

     p= -1;
     t_fin_bifu= 1500.0;
     step1base= ( 1.0/distp(1.0,0.0) )/50.0;	// 内側のループ
     step2base= ( 1.0/distp(0.0,1.0) )/50.0;	// 外側のループ
     //step2base= ( 1.0/distp(0.0,1.0) )/10.0;
     step1= step1base;  step2= step2base;

// ****** 分岐図の可変パラメータ。...変更するときは,下の distp()も変更。***
// PRM1
#define PRM1      TdoD1
#define PRM1STR  "TdoD1"
// PRM2
#define PRM2      TV1
#define PRM2STR  "TV1"
     
     printf("# free parameters are %s and %s.\n",PRM1STR,PRM2STR);
     ftrs5(t,x);  rem_prm_output();		// globalパラメータの値を出す
     for (PRM2= 1.2;  PRM2 > 0.05-CHOTTO;  PRM2= PRM2-step2){ // prm2 : 外側
	  x= x_ini;  step1= step1base;
	  fname= "/dev/null";  // 分岐図はどこにも出力しない。
	  // fname= sprintf("bifudiag-%s=%g.dat",PRM2STR,PRM2)"; // 分岐図をファイルへ。
	  fd= fopen(fname,"w");
	  for (PRM1= 12; PRM1>1.0-CHOTTO; PRM1= PRM1-step1){// prm1: 内側
	       if (p==0){ x(Nw1) = x(Nw1)+1.0e-8; } // 平衡点からの分岐を促進
	       for(t=0.0; ;){ // -------------------------------------------
		    x = rngkut4(t,x,ftrs5,h);  t = t + h;
		    //{t, x} = rk4(ftrs5,t,x,h);	// １ステップ進める。
		    tmpx = x; // xが間違って変更されないように tmpx に複製
		    // 定常状態に収束 または 既定の時間に達した場合 conv==1.
		    {conv, p}= bifudiag_output(PRM1, secf(x), outvf(x), t, t_fin_bifu, is_identical, tmpx, fd);
		    if ( x(Ndlt1)>1000.0 ){ p= -2;  conv=1; } // 脱調
		    if ( conv==1 ){ break; }
	       } // --------------------------------------------------------
	       printf("t=%g\t",t); // **************************************
	       { PRM1, x, p, step1, step2, t_fin_bifu } =
		    record_2p_bifurcation_point(p, PRM1, PRM2, step1, step2, tmpx, t_fin_bifu);
	       if (p==-2){ break; }  // 脱調の時は縦軸のパラメータも次の値へ
	  } // --- prm1 loop ---
	  fclose(fd);
	  if (fname!= "/dev/null"){ makepsfile(fname); }
     } // --- prm2 loop ---

     output_2p_bifurcation_point_to_matfiles();

     printf("# CPU time= %g(%gh)\n",gettimer(),gettimer()/3600.0);
}


// パラメータ平面上での距離を返す
Func Real distp(PRM1, PRM2)
     Real PRM1, PRM2;
{
     Real x, y;

     x= PRM1/10.0;
     y= PRM2/1.0;

     return sqrt(x*x+y*y);
}


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



Func void makepsfile(fname)
     String fname;
{
     static Integer fd;
     static Integer notfirsttime;
     String gpfile;

     if (notfirsttime==0){
//	  gpfile= "/tmp/2pbifu-makeps.gp";
	  gpfile= "2pbifu-makeps.gp";
	  fd= fopen(gpfile,"w");
	  fprintf(fd,"set term post\n");
     }
//     fprintf(fd,"set term post\nset out \""+fname+".ps\"\n");
     fprintf(fd,"p \""+fname+"\"\n");
//     fclose(fd);
//     system("gnuplot "+gpfile+" > "+fname+".ps");
//     system("rm "+gpfile);
     notfirsttime=1;
}