#define LASTSEC_NUM 40
#define LASTXSTACK_NUM 200
#define TRANSIENT_TIME (tfin*0.5)
// 周期を数えるときは点が同一かどうかの判定を甘くする
#define PERIOD_COUNT_MUL 50.0

Integer all_identical();
Integer count_sec_lastx();

Func List bifudiag_output(fprm, ref, outv, t, tfin, is_identical, xx, fd)
     Real fprm, ref, outv, t, tfin;
     Integer is_identical();
     Matrix xx;
     Integer fd;
{
     static Real lastref, lastoutv, lastt, lastfprm;
     Real  sec_t, sec_outv;
     static Matrix sec_lastx;
     static Real sec_lastt, sec_lastfprm;
     static Integer sec_co, sec_lastx_filled;
     Matrix tmpx, sec_x;
     static Matrix lastx;
     static Integer lastx_co, lastx_filled;
     static Matrix lastx_stack;
     Integer k, p;
     tmpx= xx; // xxを間違って変更しないように tmpx に複製して使う
     // はじめての場合は初期化
     if (length(sec_lastx)==0){  sec_lastx=Z(length(tmpx),LASTSEC_NUM);  }
     if (length(lastx_stack)==0){ lastx_stack=Z(length(tmpx),LASTXSTACK_NUM); }

     // 過渡状態の時はなにもしない
     if ( t<=TRANSIENT_TIME ){ return {0, -1}; }

     if ( fprm!=lastfprm ){
	  // 今のパラメータ値が初めてなら初期化
	  lastx_co= 0;
	  lastx_filled= 0;
     }else{
	  // lastx_stack に覚えた点全部が同一と判断される場合は平衡点。
	  if ( lastx_filled && all_identical(tmpx,lastx_stack,is_identical) ){
	       fprintf(fd,"%g %g %g %g\n",fprm,outv,t-lastt,t);
	       fprintf(fd,"##%g %g %g %g p= 0\n",fprm,outv,t-lastt,t);
//	       printf("##%g %g %g %g p= 0\n",fprm,outv,t-lastt,t); ///////////////////
	       return {1, 0};
	  }
//	  printf("-- %g %g\n",lastref,ref); ///////////////////////////////
	  // 断面上での値
	  if ( lastref>0.0 && ref<=0.0 ){
	       // 線形補間
	       sec_outv= (lastoutv*ref-outv*lastref)/(ref-lastref);
	       sec_x= (lastx*ref-tmpx*lastref)/(ref-lastref);
	       sec_t= (lastt*ref-t*lastref)/(ref-lastref);
	       if ( fprm != sec_lastfprm ){
		    //断面上では，今のパラメータ値が初めてならば
		    sec_co= 0;
		    sec_lastx_filled= 0;
	       }else{
		    fprintf(fd,"%g %g %g   %g\n",fprm,sec_outv,sec_t-sec_lastt,sec_t);
//		    printf("%g %g %g   %g\n",fprm,sec_outv,sec_t-sec_lastt,sec_t); //////////////////////////////////////////////////
	       }
	       if (sec_co>LASTSEC_NUM){ sec_lastx_filled= 1; }
	       sec_lastx= shiftRight(sec_lastx);  sec_lastx(:,1)= sec_x;
	       sec_co++; 
	       //
	       sec_lastt= sec_t;
	       sec_lastfprm= fprm;
	  } // if ( lastref>0.0 && ref<=0.0 )
     } // if ( fprm!=lastfprm)

     // 既定の時間(tfin)が経った場合
     if (t>tfin+1.0e-10){ 
	  if (sec_lastx_filled){ // 断面上の点が十分多く取れている場合
	       p= count_sec_lastx(sec_x, sec_lastx, is_identical, PERIOD_COUNT_MUL);
	       fprintf(fd,"## prm=%g t=%g p=%d\n",fprm,t,sec_t, p);
//	       printf("## prm=%g t=%g p=%d\n",fprm,t,sec_t, p); ///////////////
//	       print sec_lastx >> String(fprm)+".mat"; //////////////
	       return {1, p};
	  }else{
	       return {1, -99}; // どういう状態?
	  }
	       
     }
     lastref= ref;
     lastoutv= outv;
     lastt= t;
     lastfprm = fprm;
     lastx = tmpx;
     //
     if (lastx_co>Cols(lastx_stack)){ lastx_filled= 1; }
     lastx_stack= shiftRight(lastx_stack);  lastx_stack(:,1)= lastx;
     lastx_co++;

     return {0, -1};
}

Func Integer all_identical(x0, xstack, is_identical)
     // 全ての点が identical と判断されたときのみ真を返す。
     // すなわち、周期解については２周期以上は収束判定をあきらめる。
     Matrix x0, xstack;
     Integer is_identical();
{
     Integer k;

     for(k=1; k<=Cols(xstack); k++){
	  if ( !is_identical(x0, xstack(:,k), 1.0) ){
	       return (1!=1);
	  }
     }
     return (1==1);
}



Func Integer count_sec_lastx(sec_x, sec_lastx, is_identical, mul)
     Matrix sec_x, sec_lastx;
     Integer is_identical();
     Real mul;
{
     Matrix sectmp;
     Integer k, m;

     sectmp= [sec_x, sec_lastx ];
     
     for(k=1; k<=LASTSEC_NUM; k++){
	  if (k>= Cols(sectmp)){ break; }
	  for(m=k+1; m<=(LASTSEC_NUM+1); m++){
	       if (m > Cols(sectmp)){ break; }
	       if ( is_identical(sectmp(:,k),sectmp(:,m),mul) ){
		    sectmp= [ sectmp(:,1:(m-1)), sectmp(:,(m+1):) ];
		    m--;
	       }
	  }
     }
     return Cols(sectmp);
}