// This may look like C code, but it is really -*- C++ -*-
// RTSS --- Railway Total System Simulator
// (c) TAKAGI Ryo (at Kogakuin University)
// simulate.cc --- simulation routine
// -----
// ChangeLog:
// 2008. 1. 18
//  Cleared up the CSV file write routine in the feeder::simccalresult()
//  member function.
// 2007. 11. 21
//  Redirected log messages to l_file.
// -----


#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <iostream>
#include <fstream>

#include "kidenrt.hh"
#include "StrPrintf.hh"

using std :: cerr ;
using std :: endl ;
using std :: ofstream ;


FILE * estimss ;

//#define DEBUG

#define TPANVPSEP 50    // パンタ点電圧 histogram は TPANVPSEP Ｖ おきに
#define TPANVPLOW 900   // TPANVPLOW Ｖ から
#define TPANVPHIGH 1800 // TPANVPHIGH Ｖ まで計算
const int TPANVPSTEP = ( TPANVPHIGH - TPANVPLOW ) / TPANVPSEP + 2 ;
                        // step: 10以上あることが望ましい

double tpanvp [ TPANVPSTEP ] ;    // 力行車パンタ点電圧ヒストグラム
double tpanvpeak ;           // 力行車パンタ点電圧ピーク
double tpanvbottom ;         // 力行車パンタ点電圧ボトム
double tpanvmean = 0.0 ;     // 力行車パンタ点電圧平均値
double tpanvmv ;
double tpanvmvt ;
bool sw_tpanvpbset ;         // 力行車パンタ点電圧ピーク・ボトムセット済み

double tregendead = 0.0 ;    // 回生失効時間

#define TSSCSEP 500     // 変電所電流 histogram は TSSCSEP Ａ おきに
#define TSSCLOW -5000   // TSSCLOW - TSSCSEPＡ から
#define TSSCHIGH 16000  // TSSCHIGH Ａ まで計算
const int TSSCSTEP = ( TSSCHIGH - TSSCLOW ) / TSSCSEP + 2 ;
                        // step: 10以上あることが望ましい
double tssc [ TSSCSTEP ] ;                // 変電所電流のヒストグラム

#define GRTSSCSEP 500   // 合計回生電流 histogram は GRTSSCSEP Ａ おきに
#define GRTSSCLOW -15000    // GRTSSCLOW - GRTSSCSEPＡ から
#define GRTSSCHIGH 0    // GRTSSCHIGH Ａ まで計算
const int GRTSSCSTEP = ( GRTSSCHIGH - GRTSSCLOW ) / GRTSSCSEP + 2 ;
                        // step: 10以上あることが望ましい
double grtssc [ GRTSSCSTEP ] ;            // 変電所合計回生電流 histogram

double rpower [ 2 ] , rpsupp [ 2 ] ;     // 主回路回生電力, 回生可能電力
double rtpower [ 2 ] , ptpower [ 2 ] ;   // 回生電力，力行電力
double ppower [ 2 ] , spower [ 2 ] ;     // 主回路力行電力, 変電所入力電力
double srpower [ 2 ] ;                   // 変電所回生電力
double mpower [ 2 ] ;                    // 車両のパワー（動輪周）
double mbpower [ 2 ] ;                   // 車両のブレーキパワー（動輪周）
double mepower [ 2 ] ;                   // mpower のうち電気分
double mbepower [ 2 ] ;                  // mbpower のうち電気分
double renrgy = 0.0 ;                    // 主回路回生エネルギー
double resupp = 0.0 ;                    // 回生可能エネルギー
double penrgy = 0.0 ;                    // 主回路力行エネルギー
double rtenrgy = 0.0 ;                   // 回生エネルギー
double ptenrgy = 0.0 ;                   // 力行エネルギー
double senrgy = 0.0 ;                    // 変電所入力エネルギー
double srenrgy = 0.0 ;                   // 変電所回生エネルギー
double menrgy = 0.0 ;                    // 車両のエネルギー（動輪周）
double meenrgy = 0.0 ;                   // menrgy の電気分
double mbenrgy = 0.0 ;                   // 車両のブレーキエネルギー（動輪周）
double mbeenrgy = 0.0 ;                  // mbenrgy の電気分
double * sqc1 , * nsqc1 , * psqc1 ;      // 電流の２乗値(+/-, -, +)
double * sqc2 , * nsqc2 , * psqc2 ;      // 電流の２乗値(+/-, -, +)
double * rmsc , * nrmsc , * prmsc ;      // ＲＭＳ電流値(+/-, -, +)
double ptime = 0.0 ;                     // 力行時間
double xptime = 0.0 ;                    // 「力行」状態である時間
double ddhh = 0.0 ;                      // 積分時間、最初は 0.0

void prhistog ( double , int , int , int , double * , ofstream & ) ;

// Initialization routine (non-member function)
void
shokika
( int nss )
{

  // Initialising substation current histogram
  for ( int i = 0 ; i < TSSCSTEP ; ++ i )
    tssc [ i ] = 0.0 ;

  // Initialising substation total regen current histogram
  for ( int i = 0 ; i < GRTSSCSTEP ; ++ i )
    grtssc [ i ] = 0.0 ;

  // Initialising pantograph voltage histogram of accelerating trains
  for ( int i = 0 ; i < TPANVPSTEP ; ++ i )
    tpanvp [ i ] = 0.0 ;

  // Regenerated main circuit power, "regeneratable" power
  rpower [ 0 ] = rpsupp [ 0 ] = 0.0 ;

  // Main circuit accelerating power consumption, substation input power
  ppower [ 0 ] = spower [ 0 ] = 0.0 ;

  // Regenerating power, accelerating power
  rtpower [ 0 ] = ptpower [ 0 ] = 0.0 ;

  // Substation regen power, train pantograph voltage
  srpower [ 0 ] = tpanvmv = 0.0 ;

  // Power at wheel rim
  mpower [ 0 ] = 0.0 ;
  mepower [ 0 ] = 0.0 ;

  // Brake power at wheel rim
  mbpower [ 0 ] = 0.0 ;
  mbepower [ 0 ] = 0.0 ;

  // Square of currents
  sqc1 = new double [ nss ] ;
  sqc2 = new double [ nss ] ;

  // RMS currents
  rmsc = new double [ nss ] ;

  // Square of negative currents
  nsqc1 = new double [ nss ] ;
  nsqc2 = new double [ nss ] ;

  // RMS negative currents
  nrmsc = new double [ nss ] ;

  // Square of positive currents
  psqc1 = new double [ nss ] ;
  psqc2 = new double [ nss ] ;

  // RMS positive currents
  prmsc = new double [ nss ] ;

  for ( int i = 0 ; i < nss ; ++ i )
  {
    rmsc [ i ] = sqc1 [ i ] = sqc2 [ i ] = 0.0 ;
    nrmsc [ i ] = nsqc1 [ i ] = nsqc2 [ i ] = 0.0 ;
    prmsc [ i ] = psqc1 [ i ] = psqc2 [ i ] = 0.0 ;
  }

  sw_tpanvpbset = false ;
}


// Substitution function (non-member)
void
irekae
( double & a , double & b )
{
  a = b ;
  b = 0.0 ;
}


void sepstep
 ( double clval ,
   int lw ,
   int sep ,
   int stp ,
   double * x )
{
  double fld = clval - double ( lw ) + double ( sep ) ;
  int n = ( int ) floor ( fld / sep ) ;
  if ( n < 0 )
  {
    n = 0 ;
  }
  if ( n > stp - 1 )
  {
    n = stp - 1 ;
  }
#ifdef DEBUG_X
  l_ofs << StrPrintf ( "sepstep, n, step = %i %i" , n , stp ) << endl ;
#endif
  x [ n ] += ddhh ;
}


// Subroutine around the circuit analysis
void
feeder :: simccsub
( int & cctt )
{
#ifdef DEBUG_X
  l_ofs << "before ccal..." << endl ;
#endif

  // Calculate Y-matrix
  ycal () ;

#ifdef DEBUG_X
  l_ofs << "before ycal..." << endl ;
#endif

  if ( g_sw_substn_ctl_a )
  {
    // Substation voltage control, type A
    simcc_ssctl_a ( cctt ) ;
  }
  else
  {
    simcccore ( cctt ) ;
  }

#ifdef DEBUG_X
  l_ofs << "end of ccal..." << endl ;
#endif

  if ( cctt )
  {
    // Calculation ended successfully

    // ----
    // Calculate voltage and current of the ESS
    // ----
    simCalculateESDVoltagesAndCurrents () ;

    if ( g_sw_ss_direction )
    {
      // Only when this switch is true branch current vector is calculated
      simssdirec () ;
    }
    if ( g_sw_ccalresult )
    {
      // Output circuit analysis result at each iteration
      simccalresult () ;
    }
  }
  else
  {
    // Calculation abnormally terminated
    l_ofs << "ccal failed!" << endl ;
  }
}



// Subroutine around circuit analysis, "the core"
void
feeder :: simcccore
( int & cctt )
{
  // emcc ( false ) ;

#ifdef DEBUG_X
  l_ofs << "before ccal..." << endl ;
#endif

  // 回路計算：正常終了時 ccal() の返り値 true
  cctt = ccal () ;
//   if ( ! ccal () )
//   {
//     emcc ( true ) ;
//     l_ofs << "ccal first try failed!" << endl ;
//     cctt = ccal () ;   // Re-try
//   }
//   else
//   {
//     cctt = 1 ;    // cctt = 1: 計算成功！
//   }
}



// Output circuit analysis result at each iteration
void
feeder :: simccalresult
()
{
  l_ofs << "CCALRESULT" << endl ;
  for ( int i = 0 ; i < ncss ; ++ i )
  {
    if ( i < nss )
      l_ofs << "SUB " ;
    else
      l_ofs << "TRN " ;
    l_ofs << StrPrintf ( "%3d %8.1f %8.1f" , i , getCSD ( i ) -> volt () ,
                         getCSD ( i ) -> ampere () ) ;
    if ( i < nss )
    {
      if ( g_sw_substn_ctl_a )
      {
	// 変電所リアルタイム電圧制御・タイプA
	// sub [ i ] . getsubcntl () を印刷
	l_ofs << StrPrintf ( " %8.6f" , getSUB ( i ) . getsubcntl () ) ;
      }
      if ( g_sw_ss_direction )
      {
	l_ofs << "  2" << endl ;
	for ( int j = 0 ; j < getSUB ( i ) . feedcon () . number () ; ++ j )
	{
          l_ofs << "  FLN  " ;
          l_ofs << StrPrintf ( "%3d  %8.1f  %8.1f" ,
                               getSUB ( i ) . feedcon () [ j ] . nfline ,
                               getSUB ( i ) . feedcon () [ j ] . sttcurr ,
                               getSUB ( i ) . feedcon () [ j ] . endcurr )
                << endl ;
	}
      }
      else
      {
	l_ofs << "  -1" << endl ;
      }
    }
    else
    {
      l_ofs << StrPrintf ( "  %3i %2i %5.2f %8.4f %d" ,
                           getCAR ( i - nss ) . direc () ,
                           getCAR ( i - nss ) . carstat () ,
                           getCAR ( i - nss ) . velocity () ,
                           getCAR ( i - nss ) . position () ,
                           getCAR ( i - nss ) . gradnumber () )
            << endl ;
    }
  }

  // 神の声: 野木F君、この関数はこうやって改造するんだよ -- 高木
  // CSV: 時刻, SUB, TRN ごとに電圧と電流と位置(列車のみ)を出力
  if ( use_c_file )
  {
    string __quot = g_sw_csv_without_quotation ? "" : "\"" ;
    if ( nogi_header_required )
    {
      nogi_header_required = false ;
      csvout_nogi << __quot << "Time" << __quot ;
      for ( int i = 0 ; i < ncss ; ++ i )
      {
        getCSD ( i ) -> writeHeaderToCSV ( csvout_nogi ) ;
      }
      csvout_nogi << endl ;
    }
    csvout_nogi << " " << t_sim  ;
    for ( int i = 0 ; i < ncss ; ++ i )
    {
      getCSD ( i ) -> writeToCSV ( csvout_nogi ) ;
    }
    csvout_nogi << endl ;
  }

  l_ofs << "ENDCCALRESULT" << endl ;
}



// 変電所方向別電流を得る
void
feeder :: simssdirec
()
{
  //*************************************************************
  // 接続行列 hh および ブランチアドミタンス行列 dy は ccal () 関 
  // 数の実行中に与えられるものとし，ここでは改めて計算しない    
  //*************************************************************
  // ブランチ電流ベクトル = dy * ht * ノード電圧行列 を計算

  static matrix tmp_va ;
  tmp_va . renewmn ( ncss ) ;               // 記憶領域確保

  for ( int i = 0 ; i < ncss ; ++ i )
    tmp_va ( i ) = getCSD ( i ) -> volt () ;     // 電圧を代入
  static matrix tmp_htvn ;

  ht = feedy . h () ;                     // *** ここで ht 行列の抽出 ***
  dy = feedy . branchadmittance () ;      // *** ここで dy 行列の抽出 ***

  tmp_htvn = ht * tmp_va ;
  tmp_va = dy * tmp_htvn ;     // ブランチ電流 ( この順序で計算:速い )
  bran_node ck ;

  // 電流ゼロリセット
  // Zero-resetting the directional current at substations.
  for ( int i = 0 ; i < nss ; ++ i )
  {
    for ( int j = 0 ; j < getSUB ( i ) . feedcon () . number () ; ++ j )
    {
      int i_fln = getSUB ( i ) . feedcon () [ j ] . nfline ;
      double x_pos = getSUB ( i ) . feedcon () [ j ] . pos ;
      getSUB ( i ) . directional_curr_start ( i_fln , x_pos , 0 ) ;
      getSUB ( i ) . directional_curr_end ( i_fln , x_pos , 0 ) ;
    }
  }
  
  for ( int i = 0 ; i < dy . row () ; ++ i )
  {
    // ----
    // This gives the start and end nodes of Branch i
    // ----
    ck = feedy . getElement ( i ) ;

    if ( ck . startnode < nss )
    {
      getSUB ( ck . startnode ) . directional_curr_end  // branch の起点は
	( ck . nfline , ck . startpos , tmp_va ( i ) ) ;     // 変電所の終点側！
    }
    if ( ck . endnode < nss )
    {
      getSUB ( ck . endnode ) . directional_curr_start  // branch の終点は
	( ck . nfline , ck . endpos , tmp_va ( i ) ) ;     // 変電所の起点側！
    }
  }
}



// 積分計算まわりのサブルーチン
void
feeder :: siminteg
()
{
  irekae ( rpower [ 1 ] , rpower [ 0 ] ) ;               // 回生電力
  irekae ( ppower [ 1 ] , ppower [ 0 ] ) ;               // 力行電力
  irekae ( rtpower [ 1 ] , rtpower [ 0 ] ) ;             // 回生電力
  irekae ( ptpower [ 1 ] , ptpower [ 0 ] ) ;             // 力行電力
  irekae ( rpsupp [ 1 ] , rpsupp [ 0 ] ) ;               // 回生可能電力
  irekae ( mpower [ 1 ] , mpower [ 0 ] ) ;	// 動輪周パワー
  irekae ( mepower [ 1 ] , mepower [ 0 ] ) ;	// mpower のうち電気分
  irekae ( mbpower [ 1 ] , mbpower [ 0 ] ) ;	// 動輪周ブレーキパワー
  irekae ( mbepower [ 1 ] , mbepower [ 0 ] ) ;	// mbpower のうち電気分
  for ( int i = 0 ; i < ncar ; ++ i )
  {
#ifdef MORIYAIEN
    if ( getCAR ( i ) . position () >= MORIYA )
    {
      mpower [ 0 ] += getCAR ( i ) . getMechanicalPower () ;
      mepower [ 0 ] += getCAR ( i ) . getMechanicalPower ()
        * getCAR ( i ) . getRatioOfElectricEffortInMechPower () ;
      if ( getCAR ( i ) . getMechanicalPower () < 0 )
      {
        mbpower [ 0 ] -= getCAR ( i ) . getMechanicalPower () ;
        mbepower [ 0 ] -= getCAR ( i ) . getMechanicalPower ()
          * getCAR ( i ) . getRatioOfElectricEffortInMechPower () ;
      }
    }
    if ( getCAR ( i ) . rnotchvl () < 0.0
	 && getCAR ( i ) . position () >= MORIYA )
    {
      rpower [ 0 ]
        += getCAR ( i ) . volt () * getCAR ( i ) . sibori () / 1000.0 ;
      rpsupp [ 0 ]
        += getCAR ( i ) . volt () * getCAR ( i ) . asibori () / 1000.0 ;
    }
    else if ( getCAR ( i ) . rnotchvl () > 0.0
	      && getCAR ( i ) . position () >= MORIYA )
    {
      ppower [ 0 ] -= getCAR ( i ) . volt () * getCAR ( i ) . ampere () / 1000.0 ;
    }
#else
    mpower [ 0 ] += getCAR ( i ) . getMechanicalPower () ;
    mepower [ 0 ] += getCAR ( i ) . getMechanicalPower ()
      * getCAR ( i ) . getRatioOfElectricEffortInMechPower () ;
    if ( getCAR ( i ) . getMechanicalPower () < 0 )
    {
      mbpower [ 0 ] -= getCAR ( i ) . getMechanicalPower () ;
      mbepower [ 0 ] -= getCAR ( i ) . getMechanicalPower ()
        * getCAR ( i ) . getRatioOfElectricEffortInMechPower () ;
    }
    if ( getCAR ( i ) . rnotchvl () < 0.0 )
    {
      rpower [ 0 ]
        += getCAR ( i ) . volt () * getCAR ( i ) . sibori () / 1000.0 ;
      rpsupp [ 0 ]
        += getCAR ( i ) . volt () * getCAR ( i ) . asibori () / 1000.0 ;
    }
    else if ( getCAR ( i ) . rnotchvl () > 0.0 )
    {
      ppower [ 0 ] -= getCAR ( i ) . volt () * getCAR ( i ) . sibori () / 1000.0 ;
    }
#endif
    if ( getCAR ( i ) . ampere () > 0.0 )
    {
      rtpower [ 0 ] += getCAR ( i ) . volt () * getCAR ( i ) . ampere () / 1000.0 ;
    }
    else
    {
      ptpower [ 0 ] -= getCAR ( i ) . volt () * getCAR ( i ) . ampere () / 1000.0 ;
    }
  }

  irekae ( spower [ 1 ] , spower [ 0 ] ) ;               // 変電所入力電力
  irekae ( srpower [ 1 ] , srpower [ 0 ] ) ;             // 変電所回生電力
  for ( int i = 0 ; i < nss ; ++ i )
  {
    irekae ( getSUB ( i ) . power () [ 1 ] , getSUB ( i ) . power () [ 0 ] ) ;
    irekae ( getSUB ( i ) . rpower () [ 1 ] , getSUB ( i ) . rpower () [ 0 ] ) ;
    double suba = getSUB ( i ) . ampere () ;
    spower [ 0 ] += getSUB ( i ) . volt () * suba / 1000.0 ;
    if ( suba < 0.0 )
    {
      srpower [ 0 ] -= getSUB ( i ) . volt () * suba / 1000.0 ;
    }
    sqc2 [ i ] = sqc1 [ i ] ;                      // 変電所ＲＭＳ電流
    sqc1 [ i ] = suba * suba ;
    psqc2 [ i ] = psqc1 [ i ] ;                    // 変電所ＲＭＳ電流
    if ( suba >= 0.0 )                        // ( positive )
      psqc1 [ i ] = suba * suba ;
    else psqc1 [ i ] = 0.0 ;
    nsqc2 [ i ] = nsqc1 [ i ] ;                    // 変電所ＲＭＳ電流
    if ( suba < 0.0 )                         // ( negative )
      nsqc1 [ i ] = suba * suba ;
    else nsqc1 [ i ] = 0.0 ;
    getSUB ( i ) . power () [ 0 ]                       // 変電所毎エネルギー
      = getSUB ( i ) . volt () * suba / 1000.0 ;    // の入出力
    if ( suba < 0.0 )
      getSUB ( i ) . rpower () [ 0 ] = getSUB ( i ) . volt () * suba / 1000.0 ;
    getSUB ( i ) . energy () += ddhh
      * ( getSUB ( i ) . power () [ 0 ] + getSUB ( i ) . power () [ 1 ] ) / 2.0 ;
    getSUB ( i ) . renergy () += ddhh
      * ( getSUB ( i ) . rpower () [ 0 ] + getSUB ( i ) . rpower () [ 1 ] ) / 2.0 ;
    if ( g_sw_ss_direction )
    {
      // 変電所方向別エネルギーの入出力
      for ( int j = 0 ; j < getSUB ( i ) . feedcon () . number () ; ++ j )
      {
	irekae ( getSUB ( i ) . feedcon () [ j ] . sttpower [ 1 ] ,
		 getSUB ( i ) . feedcon () [ j ] . sttpower [ 0 ] ) ;
	irekae ( getSUB ( i ) . feedcon () [ j ] . endpower [ 1 ] ,
		 getSUB ( i ) . feedcon () [ j ] . endpower [ 0 ] ) ;
	// towards the beginning of line
        getSUB ( i ) . feedcon () [ j ] . sttpower [ 0 ] =
          getSUB ( i ) . feedcon () [ j ] . sttcurr * getSUB ( i ) . volt () ;
        getSUB ( i ) . feedcon () [ j ] . sttcurr = 0.0 ;
        getSUB ( i ) . feedcon () [ j ] . sttenrgy
          += ddhh * 
          ( getSUB ( i ) . feedcon () [ j ] . sttpower [ 0 ]
            + getSUB ( i ) . feedcon () [ j ] . sttpower [ 1 ] ) / 2000.0 ;
	// towards the end of line
        getSUB ( i ) . feedcon () [ j ] . endpower [ 0 ] =
          getSUB ( i ) . feedcon () [ j ] . endcurr * getSUB ( i ) . volt () ;
        getSUB ( i ) . feedcon () [ j ] . endcurr = 0.0 ;
        getSUB ( i ) . feedcon () [ j ] . endenrgy
          += ddhh * 
          ( getSUB ( i ) . feedcon () [ j ] . endpower [ 0 ]
            + getSUB ( i ) . feedcon () [ j ] . endpower [ 1 ] ) / 2000.0 ;
      }
    }
  }

  // 台形法による積分
  renrgy += ddhh * ( rpower [ 0 ] + rpower [ 1 ] ) / 2.0 ;
  resupp += ddhh * ( rpsupp [ 0 ] + rpsupp [ 1 ] ) / 2.0 ;
  penrgy += ddhh * ( ppower [ 0 ] + ppower [ 1 ] ) / 2.0 ;
  senrgy += ddhh * ( spower [ 0 ] + spower [ 1 ] ) / 2.0 ;
  rtenrgy += ddhh * ( rtpower [ 0 ] + rtpower [ 1 ] ) / 2.0 ;
  ptenrgy += ddhh * ( ptpower [ 0 ] + ptpower [ 1 ] ) / 2.0 ;
  srenrgy += ddhh * ( srpower [ 0 ] + srpower [ 1 ] ) / 2.0 ;
  menrgy += ddhh * ( mpower [ 0 ] + mpower [ 1 ] ) / 2.0 ;
  mbenrgy += ddhh * ( mbpower [ 0 ] + mbpower [ 1 ] ) / 2.0 ;
  meenrgy += ddhh * ( mepower [ 0 ] + mepower [ 1 ] ) / 2.0 ;
  mbeenrgy += ddhh * ( mbepower [ 0 ] + mbepower [ 1 ] ) / 2.0 ;
  for ( int i = 0 ; i < nss ; ++ i ) {
    rmsc [ i ] += ddhh * ( sqc1 [ i ] + sqc2 [ i ] ) / 2.0 ;
    prmsc [ i ] += ddhh * ( psqc1 [ i ] + psqc2 [ i ] ) / 2.0 ;
    nrmsc [ i ] += ddhh * ( nsqc1 [ i ] + nsqc2 [ i ] ) / 2.0 ;
  }

  // 力行車パンタ点電圧
  tpanvmv = 0.0 ;
  tpanvmvt = 0.0 ;
  for ( int i = 0 ; i < ncar ; ++ i )
  {
    // 「力行」状態である時間
    if ( getCAR ( i ) . carstat () == Power )
      xptime += ddhh ;

    // rnotch 正のときのみ力行時間
    if ( getCAR ( i ) . rnotchvl () > 0.0 )
    {
      ptime += ddhh ;
      tpanvmvt += ddhh ;
      tpanvmv += getCAR ( i ) . volt () ;
      sepstep ( getCAR ( i ) . volt () , TPANVPLOW , TPANVPSEP ,
		TPANVPSTEP , tpanvp ) ;
      if ( sw_tpanvpbset )
      {
	if ( getCAR ( i ) . volt () < tpanvbottom )
	{
	  tpanvbottom = getCAR ( i ) . volt () ;
	}
	if ( getCAR ( i ) . volt () > tpanvpeak )
	{
	  tpanvpeak = getCAR ( i ) . volt () ;
	}
      }
      else
      {
	sw_tpanvpbset = true ;
	tpanvbottom = tpanvpeak = getCAR ( i ) . volt () ;
      }
    }
  }
  tpanvmean += ddhh * tpanvmv ;
  setpeak () ;                                  // 変電所電流のピーク
  double grtsx = 0.0 ;

  // 変電所電流のヒストグラム
  for ( int i = 0 ; i < nss ; ++ i )
  {
    //      if ( getSUB ( i ) . ampere () > ispeak )
    //          ispeak = getSUB ( i ) . ampere () ;
    sepstep ( getSUB ( i ) . ampere () , TSSCLOW , TSSCSEP , TSSCSTEP , tssc ) ;
    if ( getSUB ( i ) . ampere () < 0 )
      grtsx += getSUB ( i ) . ampere () ;
  }
  sepstep ( grtsx , GRTSSCLOW , GRTSSCSEP , GRTSSCSTEP , grtssc ) ;
}



// シミュレーションプログラムの本体
void
feeder :: simulation
()
{

  deltah = 0.25 ;                                      // 時間刻み幅
  cal_deltah = deltah ;        // 計算用時間刻み幅
  readfile () ;						// 変数初期化
  disttrain () ;					// 列車配置初期化

  t_sim = 0.0 ;                                        // 時刻初期化
  double xx , xx_stt ;
  for ( int k = 0 ; k < ndp ; ++ k )
  {
    getDPT ( k ) . prtinit () ;
  }

  // simulation 開始条件まで空回し
  l_ofs << "free run till simulation starts..." << endl ;
  g_sw_freerun_cycle = true ;
  if ( nfree == 0 )
  {
    g_sw_freerun_cycle = false ;
    goto FREEEND ;
  }
  if ( ! g_sw_nosetviold )
  {
    for ( int k = 0 ; k < ncss ; ++ k )
    {
      getCSD ( k ) -> teta () = 0 ;
      getCSD ( k ) -> teta_vi () ;
    }
    setviold () ;
  }
  while ( 1 )
  {
    l_ofs << StrPrintf ( "t_sim = %.2f" , t_sim ) << endl ;
    int cctt ;
    simccsub ( cctt ) ;                           // 回路計算
    if ( ! cctt )
    {
      // 計算異常終了時
      cal_deltah += deltah ;
      restorevi () ;
    }
    else
    {
      // -----
      // Renew the electrical state variables that all objects have.
      // すべてのオブジェクトが保有する電気的状態変数を更新。
      // -----
      simRenewElectricalStates () ;
      cal_deltah = deltah ;                              // 積分刻み幅復元
#ifdef DEBUG_X
      l_ofs << "before tellimp..." << endl ;
#endif
      tellimp () ;                  // 列車の前後SS . の Cc ,Cl 値教える
#ifdef DEBUG_X
      l_ofs << "after tellimp , before autodreg" << endl ;
#endif
      for ( int j = 0 ; j < ncar ; ++ j )
      {	// 列車 V1 , V2 値を Cc ,Cl 値で変える
	getCAR ( j ) . autodreg () ;
        getCAR ( j ) . regendead ( ddhh , tregendead ) ;
      }
    }
#ifdef DEBUG_X
    l_ofs << "after autodreg" << endl ;
#endif
    if ( ! g_sw_nosetviold )
    {
      setviold () ;
    }
#ifdef DEBUG_X
    l_ofs << "before runtrain ... " << endl ;
#endif
    runtrain () ;                                     // 列車を動かす
#ifdef DEBUG_X
    l_ofs << "runtrain end" << endl ;
#endif
    for ( int j = 0 ; j < ndp ; ++ j )
      getDPT ( j ) . setbuf () ;
    t_sim += deltah ;
#ifdef DEBUG_X
    l_ofs << "time incremented ..." << endl ;
#endif

    // シミュレーションサイクル数
    if ( t_sim >= t_cyc * n_simcycle )
      ++ n_simcycle ;

    // 空回しサイクル数
    if ( t_sim >= t_cyc * nfree )
    {
      g_sw_freerun_cycle = false ;
      break ;
    }
  }

 FREEEND :	// いきなりここに来ることもある
  // ----
  // Setting origin of distance.
  // 距離原点セット
  // ----
  xx_stt = 0.0 ;
  for ( int i = 0 ; i < ncar ; ++ i )
  {
    xx_stt += getCAR ( i ) . runlen () ;
  }

  // ----
  // Setting ESS Initial SOC.
  // ----
  for ( int i = 0 ; i < ncss ; ++ i )
  {
    getCSD ( i ) -> setStartSOC () ;
  }

  // simulation 開始
  shokika ( nss ) ;
  bool last = false ;                  // 回路計算だけで exit なら true
  double tsim_stt = t_sim ;            // シミュレーション開始時刻
  while ( 1 )
  {
    l_ofs << StrPrintf ( "t_sim = %.2f" , t_sim ) << endl ;
    int cctt ;
    simccsub ( cctt ) ;                               // 回路計算
    if ( ! cctt )
    {
      // 計算異常終了時は積分しない
      ddhh += deltah ;
      cal_deltah = ddhh ;
      restorevi () ;
    }
    else
    {
      // -----
      // Renew the electrical state variables that all objects have.
      // すべてのオブジェクトが保有する電気的状態変数を更新。
      // -----
      simRenewElectricalStates () ;
      tellimp () ;                  // 列車前後SS . の Cc ,Cl 値教える
      for ( int j = 0 ; j < ncar ; ++ j )
      {
	// 列車 V1 , V2 値を Cc ,Cl 値で変更 ( Ea_cvauto0 用 )
	getCAR ( j ) . autodreg () ;

        getCAR ( j ) . regendead ( ddhh , tregendead ) ;
      }
      siminteg () ;                 // シミュレーション結果積分 ( 台形法 )
      ddhh = deltah ;                              // 積分刻み幅復元
      cal_deltah = ddhh ;
    }

    // 最後 && 回路演算正常終了なら終了
    if ( last && cctt )
      break ;

    for ( int j = 0 ; j < ndp ; ++ j )
      getDPT ( j ) . setbuf () ;
    if ( ! g_sw_nosetviold )
    {
      setviold () ;                 // 追加 92.10.13
    }

    runtrain () ;                                     // 列車を動かす
    t_sim += deltah ;

    // シミュレーションサイクル数
    if ( t_sim >= t_cyc * n_simcycle )
      ++ n_simcycle ;

    // 終了条件の判定 ( サイクル数で )
    if ( t_sim >= t_cyc * ( nfree + nsim )
	 + ( nfree == 0 ? deltah : 0.0 ) )
      last = true ;
  }

  // シミュレーション距離の計算
  xx = 0.0 ;
  for ( int i = 0 ; i < ncar ; ++ i )
    xx += getCAR ( i ) . runlen () ;

  // 計算結果の出力
  double ts_l = t_sim - tsim_stt ;
  ofstream fresult ;
  fresult. open ( r_file ) ;
  if ( fresult . fail () )
  {
    cerr  << "Error: Result file \"" << r_file << "\" open failed." << endl ;
    l_ofs << "Error: Result file \"" << r_file << "\" open failed." << endl ;
    exit ( 1 ) ;
  }
  fresult
    << "     ***** SIMULATION RESULT by RTSS ver.2.6 *****" << endl << endl
    << StrPrintf ( "Total Powering time          =%9.2f" , ptime ) << endl
    << StrPrintf ( "Total Power-Status time      =%9.2f" , xptime ) << endl ;
  if ( g_sw_per_car_out )
  {
    fresult
      << StrPrintf ( "Regeneratable energy per car =%11.3f MJ /%12.3f kW" ,
                     resupp / ( 1000.0 * nsim * ncar ) ,
                     resupp / ( ts_l * nsim * ncar ) ) << endl
      << StrPrintf ( "Regenerated energy per car   =%11.3f MJ /%12.3f kW" ,
                     rtenrgy / ( 1000.0 * nsim * ncar ) ,
                     rtenrgy / ( ts_l * nsim * ncar ) ) << endl
      << StrPrintf ( " (Main circult only)         =%11.3f MJ /%12.3f kW" ,
                     renrgy / ( 1000.0 * nsim * ncar ) ,
                     renrgy / ( ts_l * nsim * ncar ) ) << endl
      << StrPrintf ( "Powering energy per car      =%11.3f MJ /%12.3f kW" ,
                     ptenrgy / ( 1000.0 * nsim * ncar ) ,
                     ptenrgy / ( ts_l * nsim * ncar ) ) << endl
      << StrPrintf ( " (Main circuit only)         =%11.3f MJ /%12.3f kW" ,
                     penrgy / ( 1000.0 * nsim * ncar ) ,
                     penrgy / ( ts_l * nsim * ncar ) ) << endl
      << StrPrintf ( "Energy at wheel rim          =%11.3f MJ /%12.3f kW" ,
                     menrgy / ( 1000.0 * nsim * ncar ) ,
                     menrgy / ( ts_l * nsim * ncar ) ) << endl
      << StrPrintf ( " (electrical only)           =%11.3f MJ /%12.3f kW" ,
                     meenrgy / ( 1000.0 * nsim * ncar ) ,
                     meenrgy / ( ts_l * nsim * ncar ) ) << endl
      << StrPrintf ( "Braking energy at wheel rim  =%11.3f MJ /%12.3f kW" ,
                     mbenrgy / ( 1000.0 * nsim * ncar ) ,
                     mbenrgy / ( ts_l * nsim * ncar ) ) << endl
      << StrPrintf ( " (electrical only)           =%11.3f MJ /%12.3f kW" ,
                     mbeenrgy / ( 1000.0 * nsim * ncar ) ,
                     mbeenrgy / ( ts_l * nsim * ncar ) ) << endl ;
  }
  else
  {
    fresult
      << StrPrintf ( "Regeneratable energy         =%11.3f MJ /%12.3f kW" ,
                     resupp / 1000.0 , resupp / ts_l ) << endl
      << StrPrintf ( "Regenerated energy           =%11.3f MJ /%12.3f kW" ,
                     rtenrgy / 1000.0 , rtenrgy / ts_l ) << endl
      << StrPrintf ( " (Main circult only)         =%11.3f MJ /%12.3f kW" ,
                     renrgy / 1000.0 , renrgy / ts_l ) << endl
      << StrPrintf ( "Powering energy              =%11.3f MJ /%12.3f kW" ,
                     ptenrgy / 1000.0 , ptenrgy / ts_l ) << endl
      << StrPrintf ( " (Main circuit only)         =%11.3f MJ /%12.3f kW" ,
                     penrgy / 1000.0 , penrgy / ts_l ) << endl
      << StrPrintf ( "Energy at wheel rim          =%11.3f MJ /%12.3f kW" ,
                     menrgy / 1000.0 , menrgy / ts_l ) << endl
      << StrPrintf ( " (electrical only)           =%11.3f MJ /%12.3f kW" ,
                     meenrgy / 1000.0 , meenrgy / ts_l ) << endl
      << StrPrintf ( "Braking energy at wheel rim  =%11.3f MJ /%12.3f kW" ,
                     mbenrgy / 1000.0 , mbenrgy / ts_l ) << endl
      << StrPrintf ( " (electrical only)           =%11.3f MJ /%12.3f kW" ,
                     mbeenrgy / 1000.0 , mbeenrgy / ts_l ) << endl ;
  }
  fresult
    << StrPrintf ( "Energy consumed by trains    =%11.3f MJ /%12.3f kW" ,
                   ( ptenrgy - rtenrgy ) / 1000.0 ,
                   ( ptenrgy - rtenrgy ) / ts_l ) << endl
    << StrPrintf ( "Feeder loss                  =%11.3f MJ /%12.3f kW" ,
                   ( senrgy - ptenrgy + rtenrgy ) / 1000.0 ,
                   ( senrgy - ptenrgy + rtenrgy ) / ts_l ) << endl
    << StrPrintf ( "Substation inflow energy     =%11.3f MJ /%12.3f kW" ,
                   ( senrgy + srenrgy ) / 1000.0 ,
                   ( senrgy + srenrgy ) / ts_l ) << endl
    << StrPrintf ( "Substation outflow energy    =%11.3f MJ /%12.3f kW" ,
                   srenrgy / 1000.0 , srenrgy / ts_l ) << endl
    << StrPrintf ( "Substation input energy      =%11.3f MJ /%12.3f kW" ,
                   senrgy / 1000.0 , senrgy / ts_l ) << endl
    << StrPrintf ( "Regeneration Invalidness Rate =%9.6f" ,
                   1.0 - renrgy / resupp ) << endl ;
  fresult << StrPrintf ( "Regeneration Disabled time    =%9.2f" ,
                         tregendead ) << endl ;
  fresult << StrPrintf ( "Regeneration Rate             =%9.6f" ,
                         rtenrgy / ptenrgy ) << endl
          << StrPrintf ( " (Main circuit only)          =%9.6f" ,
                         renrgy / penrgy ) << endl ;
  double meanvel ;
  meanvel = xx - xx_stt ;
  fresult << StrPrintf ( "Run Length = %10.3f" , meanvel ) << endl ;
  meanvel *= 3600.0 / ( ncar * t_cyc * nsim ) ;
  fresult << StrPrintf ( "Velocity = %10.3f" , meanvel ) << endl ;

  // ----
  // Output SOC-related results
  // ----
  ostringstream oss_soc ;
  ostringstream oss_head ;
  oss_soc << endl << endl
          << "Energy[SOC] of ESS's: init, max, min, final [MJ, %]" << endl ;
  for ( int i = 0 ; i < ncss ; ++ i )
  {
    oss_head . str ( "" ) ;
    if ( i < nss )
    {
      oss_head << StrPrintf ( " SUB %3d:" , i ) ;
    }
    else
    {
      oss_head << StrPrintf ( " TRN %3d:" , i ) ;
    }
    getCSD ( i ) -> printSOCStoredResults ( fresult , oss_soc , oss_head ) ;
  }

  fresult
    << endl << endl
    << "Pantograph voltage peak/bottom/mean value of powering cars:" << endl
    << StrPrintf
    ( "  %11.3f [V peak] / %11.3f [V bottom] / %11.3f [V mean]" ,
      tpanvpeak , tpanvbottom , tpanvmean / ptime ) << endl << endl
    << "Pantograph voltage of powering cars" << endl ;
  prhistog ( ptime , TPANVPLOW , TPANVPSEP , TPANVPSTEP , tpanvp , fresult ) ;
  fresult << "Substation Current Histogram" << endl ;
  double tsx = ts_l * nss ;
  prhistog ( tsx , TSSCLOW , TSSCSEP , TSSCSTEP , tssc , fresult ) ;
  fresult << "Total Substation Regenerative Current Histogram" << endl ;
  prhistog ( tsx , GRTSSCLOW , GRTSSCSEP , GRTSSCSTEP , grtssc , fresult ) ;
  fresult << "Substation Peak Current" << endl ;
  short int hyoji = 0 ;
  for ( int i = 0 ; i < nss ; ++ i )
  {
    if ( hyoji == 2 )
    {
      fresult << StrPrintf
        ( " No.%2i: %8.1f" , i , getSUB ( i ) . peakcur () ) << endl ;
      hyoji = 0 ;
    }
    else
    {
      fresult << StrPrintf
        ( " No.%2i: %8.1f" , i , getSUB ( i ) . peakcur () ) ;
      ++ hyoji ;
    }
  }
  if ( hyoji != 0 ) fresult << endl ;
  fresult << "Substation Negative Peak Current" << endl ;
  hyoji = 0 ;
  for ( int i = 0 ; i < nss ; ++ i )
  {
    if ( hyoji == 2 )
    {
      fresult << StrPrintf
        ( " No.%2i: %8.1f" , i , getSUB ( i ) . npeakcur () ) << endl ;
      hyoji = 0 ;
    }
    else
    {
      fresult << StrPrintf
        ( " No.%2i: %8.1f" , i , getSUB ( i ) . npeakcur () ) ;
      ++ hyoji ;
    }
  }
  if ( hyoji != 0 ) fresult << endl ;
  fresult << "Substation Maximum Voltage" << endl ;
  hyoji = 0 ;
  for ( int i = 0 ; i < nss ; ++ i )
  {
    if ( hyoji == 2 )
    {
      fresult << StrPrintf
        ( " No.%2i: %8.1f" , i , getSUB ( i ) . peakvolt () ) << endl ;
      hyoji = 0 ;
    }
    else
    {
      fresult << StrPrintf
        ( " No.%2i: %8.1f" , i , getSUB ( i ) . peakvolt () ) ;
      ++ hyoji ;
    }
  }
  if ( hyoji != 0 ) fresult << endl ;
  fresult << "Substation Minimum Voltage" << endl ;
  hyoji = 0 ;
  for ( int i = 0 ; i < nss ; ++ i )
  {
    if ( hyoji == 2 )
    {
      fresult << StrPrintf
        ( " No.%2i: %8.1f" , i , getSUB ( i ) . bottomvolt () ) << endl ;
      hyoji = 0 ;
    }
    else
    {
      fresult << StrPrintf
        ( " No.%2i: %8.1f" , i , getSUB ( i ) . bottomvolt () ) ;
      ++ hyoji ;
    }
  }
  if ( hyoji != 0 ) fresult << endl ;
  if ( g_sw_kilowatt_subout )
  {
    fresult << "Substation Energy in [kW] (total, inflow, outflow)" << endl ;
    for ( int i = 0 ; i < nss ; ++ i )
    {
      fresult << StrPrintf
        ( " No.%2i: %10.3f %10.3f %10.3f" , i ,
          getSUB ( i ) . energy () / ts_l , 
          ( getSUB ( i ) . energy () - getSUB ( i ) . renergy () ) / ts_l ,
          - getSUB ( i ) . renergy () / ts_l ) << endl ;
    }
  }
  else
  {
    fresult << "Substation Energy in [MJ] (total, inflow, outflow)" << endl ;
    for ( int i = 0 ; i < nss ; ++ i )
    {
      fresult << StrPrintf
        ( " No.%2i: %10.3f %10.3f %10.3f" , i ,
          getSUB ( i ) . energy () / 1000 , 
          ( getSUB ( i ) . energy () - getSUB ( i ) . renergy () ) / 1000 ,
          - getSUB ( i ) . renergy () / 1000 ) << endl ;
    }
  }
  if ( g_sw_ss_direction )
  {
    // 方向別出力をしたときのみ
    fresult << endl << "Substation Directional Energy and  Peak Current in "
            << "[MJ][A]" << endl ;
    for ( int i = 0 ; i < nss ; ++ i )
    {
      fresult << StrPrintf ( "No.%2i: 2" , i ) << endl ;
      for ( int j = 0 ; j < getSUB ( i ) . feedcon () . number () ; ++ j )
      {
        fresult
          << StrPrintf ( " Feedline No.%2i STT: %11.3f [MJ] , " , j ,
                         getSUB ( i ) . feedcon () [ j ] . sttenrgy / 1000 )
          << StrPrintf ( "%11.3f [A] , %12.3f [A]" ,
                         getSUB ( i ) . feedcon () [ j ] . endpeak ,
                         - getSUB ( i ) . feedcon () [ j ] . endnpeak )
          << endl
          << StrPrintf ( " Feedline No.%2i END: %11.3f [MJ] , " , j ,
                         getSUB ( i ) . feedcon () [ j ] . endenrgy / 1000 )
          << StrPrintf ( "%11.3f [A] , %12.3f [A]" ,
                         getSUB ( i ) . feedcon () [ j ] . sttpeak ,
                         - getSUB ( i ) . feedcon () [ j ] . sttnpeak )
          << endl ;
      }
    }
  }
  fresult << endl << "Substation RMS Current (total, positive, negative)"
          << endl ;
  for ( int i = 0 ; i < nss ; ++ i )
  {
    fresult << StrPrintf
      ( " No.%2i: %8.1f %8.1f %8.1f" , i , sqrt ( rmsc [ i ] / ts_l ) ,
        sqrt ( prmsc [ i ] / ts_l ) , sqrt ( nrmsc [ i ] / ts_l ) ) << endl ;
  }
  fresult << endl << endl
          << StrPrintf ( "simulated time = %.3f" , ts_l ) << endl ;
  fresult . close () ;
}


// histogram 出力ルーチン
void
prhistog
( double tsx ,
  int lw ,
  int sep ,
  int stp ,
  double * tblhyoji ,
  ofstream & fresult )
{
  fresult << StrPrintf ( "      -%6i:%.6f  %6i-%6i:%.6f  %6i-%6i:%.6f" ,
                         lw , tblhyoji [ 0 ] / tsx , lw , lw + sep ,
                         tblhyoji [ 1 ] / tsx , lw + sep , lw + sep * 2 ,
                         tblhyoji [ 2 ] / tsx ) << endl ;
  int npp = int ( floor ( ( stp - 1 ) / 3 ) ) ;
  for ( int i = 1 ; i <= npp ; ++ i )
  {
    int k1 = i * 3 ;
    int k2 = k1 + 1 ;
    int k3 = k2 + 1 ;
    int v1 = lw + sep * ( k1 - 1 ) ;
    int v2 = lw + sep * k1 ;
    int v3 = lw + sep * k2 ;
    int v4 = lw + sep * k3 ;
    if ( i != npp )
      fresult << StrPrintf ( "%6i-%6i:%.6f  %6i-%6i:%.6f  %6i-%6i:%.6f" ,
                             v1 , v2 , tblhyoji [ k1 ] / tsx , v2 , v3 , 
                             tblhyoji [ k2 ] / tsx , v3 , v4 ,
                             tblhyoji [ k3 ] / tsx ) << endl ;
    else if ( k2 == stp )
      fresult << StrPrintf ( "%6i-      :%.6f" , v1 , tblhyoji [ k1 ] / tsx )
              << endl << endl ;
    else if ( k3 == stp )
      fresult << StrPrintf ( "%6i-%6i:%.6f  %6i-      :%.6f" , v1 , v2 ,
                             tblhyoji [ k1 ] / tsx , v2 ,
                             tblhyoji [ k2 ] / tsx ) << endl << endl ;
    else
      fresult << StrPrintf ( "%6i-%6i:%.6f  %6i-%6i:%.6f  %6i-      :%.6f" ,
                             v1 , v2 , tblhyoji [ k1 ] / tsx , v2 , v3 ,
                             tblhyoji [ k2 ] / tsx , v3 , tssc [ k3 ] / tsx )
              << endl << endl ;
  }
}



// -----
// Renew the electrical state variables that all objects have.
// すべてのオブジェクトが保有する電気的状態変数を更新。
// -----
void
feeder :: simRenewElectricalStates
()
{
  for ( int i = 0 ; i < ncss ; ++ i )
  {
    getCSD ( i ) -> renewElectricalStates () ;
  }
}



// -----
// Calculate voltages and currents of all ESD's
// -----
void
feeder :: simCalculateESDVoltagesAndCurrents
()
{
  for ( int i = 0 ; i < ncss ; ++ i )
  {
    getCSD ( i ) -> calculateESDVoltageAndCurrent () ;
  }
}
