#define SMALLRATE 20  /* 細かく計算するときの細かさ */

#define CHAOS 5

#define TFIN_MUL 2.0 /* 細かく計算するときに計算時間を増やす */

// このファイルで定義
Integer find_bifurcation_point();

Func List record_2p_bifurcation_point(p, prm1, prm2, step1, step2, x, tfin)
     Integer p;
     Real prm1, prm2, step1, step2;
     Matrix x;
     Real tfin;
{
     static Integer notfirsttime;
     static Real lastprm2;
     static Real prm1bak;
     static Matrix xbak;
     static Integer pbak, smallco;
     static Real step1base, step2base, tfinbase;
     static Integer lastp;
     Integer bp;

     if (notfirsttime==0){ // はじめて呼び出される時
	  // 初期化
	  prm1bak= 0.0;  xbak= x;  pbak=-1;
	  step1base= step1;  step2base= step2;  tfinbase= tfin;
     }
     // prm2（外側のループ) が更新された場合
     if (notfirsttime==0 || lastprm2 != prm2){
	  lastp= -1;
	  smallco= 0;
	  tfin= tfinbase;
     }
     notfirsttime=1;
     lastprm2= prm2;

     printf("%g\t%g\t%d\tlastp=%d\tsmallco=%d\n",prm1,prm2, p, lastp, smallco);

     bp= find_bifurcation_point(p, lastp, smallco, prm1, prm2);
     if ( smallco>0 ){ // 細かく正向きに再計算中の時
	  if (smallco==SMALLRATE+1){ // 正向き再計算の最初
	       step1= step1base/Real(SMALLRATE);
	  }
	  smallco--;
	  if (smallco==0){ // 再計算終了
	       if (bp==CHAOS){ // カオスの場合は ///
		    // もう１区間細かく計算。//
		    smallco= SMALLRATE+1;
		    lastp= p;
		    return { prm1, x, p, step1, step2, tfin };
	       }
	       step1= step1base;
	       tfin= tfinbase;
	       // 後で戻す為に覚える
	       prm1bak= prm1;  xbak= x; pbak= p; 
	       lastp= p;
	       return { prm1, x, p, step1, step2, tfin };
	  }
     }
//     if ( smallco<0 ){ // 細かく逆向きに再計算中の時
//	  smallco++;
//	  if (smallco==0){
//	       // 正向きに方向転換
//	       step1= 0.0; 
//	       smallco= SMALLRATE+1;
//	       lastp= -1;
//	       //前の状態に戻す
//	       prm1= prm1bak; x= xbak; p= pbak;
//	  }
//     }
     if ( smallco==0 && bp>0 ){
	  //分岐点があれば step を細かくして再計算
	  printf("##%g\t%g\t%d\t : bifu. point found. (%d)\n",prm1, prm2, p, bp);
//	  if (bp<CHAOS){ // 逆向きへ
//	       step1= -step1base/Real(SMALLRATE);
//	       tfin= tfinbase*TFIN_MUL;
//	       smallco= -SMALLRATE;
//	       lastp= -1;
//	  }else{ // 正向きのみ
	       step1= 0.0;
	       tfin= tfinbase*TFIN_MUL;
	       smallco= SMALLRATE+1;
	       lastp= -1;
	       //前の状態へ戻す
	       prm1= prm1bak; x= xbak; p= pbak;
//	  }
     }
     if ( smallco==0 ){
	  // 後で戻す為に覚える
	  prm1bak= prm1;  xbak= x; pbak= p; 
     }
     lastp= p;

     return { prm1, x, p, step1, step2, tfin };
}


// 分岐点の検出と記憶
Func Integer find_bifurcation_point(p, lastp, flag, prm1, prm2)
     Integer p;
     Integer lastp;
     Integer flag;
     Real prm1, prm2;
{
     static Integer notfirsttime;
     static Matrix Hopf0to1, Hopf1to0, flip1to2, flip2to1, stepout, chaos, unknown;

     if (notfirsttime==0){ // はじめて呼び出される時
	  // 初期化
	  Hopf0to1=Z(0); Hopf1to0=Z(0); flip1to2=Z(0); flip2to1=Z(0);
	  stepout=Z(0); chaos=Z(0); unknown=Z(0);
     }
     notfirsttime=1;
     if (p==-100){ //output_2p_bifurcation_point_to_matfiles() から呼ばれた場合
	  print Hopf0to1 >> "Hopf0to1.mat";
	  print Hopf1to0 >> "Hopf1to0.mat";
	  print flip1to2 >> "flip1to2.mat";
	  print flip2to1 >> "flip2to1.mat";
	  print stepout >> "stepout.mat";
	  print chaos >> "chaos.mat";
	  print unknown >> "unknown.mat";
	  return -100;
     }

     if       (lastp==0 && p==1){ // Hopf
	  if (flag){
	       Hopf0to1= [ Hopf0to1, [[prm1][prm2][p][lastp]] ];
	       printf("##%g\t%g\t\t%d  %d Hopf0to1\n",prm1,prm2,p,lastp);
	  }
	  return 1;
     }else if (lastp==1 && p==0){ // Hopf
	  if (flag){
	       Hopf1to0= [ Hopf1to0, [[prm1][prm2][p][lastp]] ];
	       printf("##%g\t%g\t\t%d  %d Hopf1to0\n",prm1,prm2,p,lastp);
	  }
	  return 2;
     }else if (lastp==1 && p==2){ // flip
	  if (flag){
	       flip1to2= [ flip1to2, [[prm1][prm2][p][lastp]] ];
	       printf("##%g\t%g\t\t%d  %d flip1to2\n",prm1,prm2,p,lastp);
	  }
	  return 3;
     }else if (lastp==2 && p==1){ // flip
	  if (flag){ 
	       flip2to1= [ flip2to1, [[prm1][prm2][p][lastp]] ];
	       printf("##%g\t%g\t\t%d  %d flip2to1\n",prm1,prm2,p,lastp);
	  }
	  return 4;
     }else if (lastp!=-1 && p>12){ // カオス
	  if (flag){
	       chaos= [ chaos, [[prm1][prm2][p][lastp]] ];
	       printf("##%g\t%g\t\t%d  %d chaos\n",prm1,prm2,p,lastp);
	  }
	  return CHAOS;
     }else if (lastp!=-1 && p==-2){ // 脱調
	  if (flag){
	       stepout= [ stepout, [[prm1][prm2][p][lastp]] ];
	       printf("##%g\t%g\t\t%d  %d stepout\n",prm1,prm2,p,lastp);
	  }
	  return 6;
     }else if (lastp!=-1 && lastp!=p){ // 何かの分岐かな?
	  if (flag){
	       unknown= [ unknown, [[prm1][prm2][p][lastp]] ];
	       printf("##%g\t%g\t\t%d  %d (unknown)\n",prm1,prm2,p,lastp);
	  }
	  return 100;
     }else{
	  return 0;
     }
}


Func void output_2p_bifurcation_point_to_matfiles(){
     Matrix darr;
     Integer dummy;
	  
     darr=Z(0);

     dummy= find_bifurcation_point(-100, darr, 0, 0.0, 0.0);
}