//---------------------------------------------------------------------------

#include <vcl.h>
#pragma hdrstop

#include "main.h"
#include "results.h"
#include "help.h"
#include <stdio.h>
#include <math.h>
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
        : TForm(Owner)
{
}
//---------------------------------------------------------------------------
void __fastcall TForm1::BitBtn_execClick(TObject *Sender)
{
double kp, k1, m1, m2, h, dt, a, b, d, z_old, z_new,
  ainv, binv, dinv, dis1=0., dis2=0. , dt_org,
  vel1=0., vel2=0., accr1=0., accr2=0. ,
  f1, f2, z, dt2, d1_old=0., d2_old=0.,
  adbc, t1, t2, x1, x2, time=0., fset, acc_temp1, acc_temp2,
  p_max, p_min, f1_old=0., f2_old=0., te, dz,
  period1, period2, h1, h2, damp[2][2],
  k2, p_dy, p_py, b_dy, b_py, ks, beta, fai1, fai2,
  acc_next1, acc_next2, dis_next1, dis_next2, vel_next1, vel_next2,
  acc1_max=0., acc2_max=0., dis1_max=0., dis2_max=0.,
  frc1_max=0., frc2_max=0.;
int n, number, ncomp, total, i, iteration, division, pier, 
    iter, bearing ;
char buffer[8] ;
const double EPS = 1.e-5 ;
const int linear=0, clough=1, bilin=2;
bool iflag;
FILE *fpin, *fpout, *fplog;

  Form2->Hide();
  Form3->Hide();
  
// パラメーターが全部入力されているか
  if ( check_1() ) { return ; }

// ファイルのオープン
// ファイルのフル・パスは，○○->Hintに保管してある

  if (( fpin=fopen(Edit_eq->Hint.c_str(),"rt")) == NULL ){
    Application->MessageBox("地震波ファイルを設定し直して下さい。",
      "入力ファイルエラー",MB_ICONEXCLAMATION | MB_OK ) ;
    return ;
    }

  if (( fpout=fopen(Edit_resp->Hint.c_str(),"wt")) == NULL ){
    Application->MessageBox("結果の出力ファイルを設定し直して下さい。",
      "出力ファイルエラー",MB_ICONEXCLAMATION | MB_OK ) ;
    return ;
  }

  if (( fplog=fopen(Edit_log->Hint.c_str(),"wt")) == NULL ){
    Application->MessageBox("ログファイルを設定し直して下さい。",
      "出力ファイルエラー",MB_ICONEXCLAMATION | MB_OK ) ;
    return ;
  }

  if ( pier_model->ItemIndex == 0 ) {
    pier = linear ;
    kp = StrToFloat( Edit_lin_kp->Text );           // 橋脚の初期剛性
    p_min = 0.;                                     // 最小橋脚変位
    p_max = 0.;                                     // 最大橋脚変位
  }
  else if ( pier_model->ItemIndex == 1 ) {
    pier = clough ;
    p_dy = StrToFloat( Edit_clo_dy->Text );         // 橋脚降伏変位
    p_py = StrToFloat( Edit_clo_py->Text );         // 橋脚降伏耐力
    kp = p_py / p_dy;                               // 橋脚の初期剛性
    p_min = -p_dy;                                  // 最小橋脚変位
    p_max =  p_dy;                                  // 最大橋脚変位
  }
  else {
    pier = bilin ;
    p_dy = StrToFloat( Edit_bil_dy->Text );         // 橋脚降伏変位
    p_py = StrToFloat( Edit_bil_py->Text );         // 橋脚降伏耐力
    kp = p_py / p_dy;                               // 橋脚の初期剛性
    ks = StrToFloat( Edit_bil_k2->Text );           // 橋脚の二次剛性
    p_min = -p_dy;                                  // 最小橋脚変位
    p_max =  p_dy;                                  // 最大橋脚変位
  }

  if ( bearing_model->ItemIndex == 0 ) {
    bearing = linear ;
    k1 = StrToFloat( Edit_b_lin_k1->Text );         // 支承の初期剛性
  }
  else {
    bearing = bilin ;
    b_dy = StrToFloat( Edit_b_bil_dy->Text );       // 支承の降伏変位
    b_py = StrToFloat( Edit_b_bil_py->Text );       // 支承の降伏耐力
    k1 = b_py / b_dy;                               // 支承の初期剛性
    k2 = StrToFloat( Edit_b_bil_k2->Text );         // 支承の２次剛性
  }

  m1 = StrToFloat( Edit_m1->Text ) ;           // 橋脚質量
  m2 = StrToFloat( Edit_m2->Text ) ;           // 桁質量  
  h1 = StrToFloat( Edit_h1->Text );            // 橋脚の減衰定数
  h2 = StrToFloat( Edit_h2->Text );            // 支承の減衰定数
  dt = StrToFloat( Edit_dt->Text );            // 時間刻み
  total = StrToInt( Edit_number->Text );       // 地震波個数
  number = int( total / 100 );                 // 計算の途中経過表示用
  beta = 1. / StrToFloat( Edit_beta->Text );   // Newmarkβのβ値
  iteration=100;                               // Newmarkβ陰解法の最大反復数
  division=StrToInt( Edit_div->Text );         // 時間刻みの細分化数

  natural( m1, m2, kp, k1, &period1, &period2, &fai1, &fai2, h1, h2, damp ) ;
                                               // 固有周期と歪みエネルギー比例減衰の計算
// 計算結果のログファイルへの出力
  fprintf(fplog,"%s\t%lf\n","m1=",m1);
  fprintf(fplog,"%s\t%lf\n","m2=",m2);
  fprintf(fplog,"%s\t%lf\n","kp=",kp);
  fprintf(fplog,"%s\t%lf\n","k1=",k1);
  fprintf(fplog,"%s\t%lf\n","h1=",h1);
  fprintf(fplog,"%s\t%lf\n","h2=",h2);
  fprintf(fplog,"%s\t%lf\n","period 1st mode=",period1);
  fprintf(fplog,"%s\t%lf\n","period 2nd mode=",period2);
  fprintf(fplog,"%s\t%lf\n","fai 11=",fai1);
  fprintf(fplog,"%s\t%lf\n","fai 12=",fai2);
  fprintf(fplog,"%s\n","damping matrix:");
  fprintf(fplog,"%lf\t%lf\n",damp[0][0],damp[0][1]);
  fprintf(fplog,"%lf\t%lf\n\n",damp[1][0],damp[1][1]);

  dt_org = dt ;
  dt = dt / (double)division ;               // division回にわけて計算
  dt2 = dt * dt;

// 出力用ファイルに各列の説明を出力、nは時間ステップ
  n = 0;
  ncomp = number;
  fprintf(fpout,"%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n",
    "時間","地震加速度(m/s/s)","橋脚変位(m)","支承変位(m)","橋脚速度(m/s)",
    "桁速度(m/s)","橋脚加速度(m/s/s)","桁加速度(m/s/s)",
    "橋脚復元力(kN)","支承復元力(kN)");

// 入力地震波データがある限り繰り返す

  z_old = 0. ;
  while( fscanf(fpin,"%lf",&z_new ) != EOF ){

// 時間をすすめ、時間経過バーを表示
    time = time + dt_org ;
    n = n + 1;
    if ( n > total ) {
      break ;
    }
    if ( n == ncomp ) {
      ProgressBar1->StepIt();
      ncomp = ncomp + number;
    }

// Newmarkβ陰解法による数値積分  division回にわけて計算

    z_new /= 100.;              // gal->m/s/s
    z = -z_old ;
    dz=(-z_new-z)/(double)division ;

    for ( i=0 ; i<division ; i++ ) {
      z+=dz;
      acc_next1=accr1;
      acc_next2=accr2;

      for ( iter=0; iter<iteration; iter++ ) {
        dis_next1=dis1+dt*vel1+(0.5-beta)*dt2*accr1+beta*dt2*acc_next1;
        dis_next2=dis2+dt*vel2+(0.5-beta)*dt2*accr2+beta*dt2*acc_next2;
        vel_next1=vel1+dt*accr1/2.+dt*acc_next1/2.;
        vel_next2=vel2+dt*accr2/2.+dt*acc_next2/2.;

// 橋脚は、剛性劣化(Clough)型かバイリニア型か線形

        if ( pier==linear ) {
          f1 = kp*dis_next1;
        }
        else if ( pier==clough ) {
          f1 = degrading( dis_next1, d1_old, f1_old, p_py, kp, p_max, p_min ) ;
        }
        else {
          f1 = bilinear( dis_next1, d1_old, f1_old, p_dy, p_py, kp, ks ) ;
        }

// 支承はバイリニアまたは線形
        if ( bearing==linear ) {
          f2 = k1* ( dis_next2-dis_next1 ) ;
        }
        else {
          f2 = bilinear( dis_next2-dis_next1, d2_old, f2_old, b_dy, b_py, k1, k2 ) ;
        }

// 運動方程式から加速度を計算
        acc_temp1=z-damp[0][0]*vel_next1/m1-damp[0][1]*vel_next2/m1-f1/m1+f2/m1;
        acc_temp2=z-damp[1][0]*vel_next1/m2-damp[1][1]*vel_next2/m2-f2/m2;

// 収束判定　iflagが真のままなら収束
        iflag=true;
        if (fabs(acc_temp1-acc_next1)>EPS) iflag=false;
        if (fabs(acc_temp2-acc_next2)>EPS) iflag=false;
        if ( iflag ) {                                   // 収束した場合
          dis1 = dis_next1; dis2 = dis_next2;
          vel1 = vel_next1; vel2 = vel_next2;
          accr1= acc_next1; accr2= acc_next2;
          d1_old = dis1;
          d2_old = dis2 - dis1;
          f1_old = f1;
          f2_old = f2;
          if ( dis1 > p_max ) {
            p_max = dis1;
          }
          if ( dis1 < p_min ) {
            p_min = dis1;
          }
          break;                                         // 収束計算のループを抜ける
        }
        else {                                           // 収束しなかった場合再計算
          acc_next1=acc_temp1;
          acc_next2=acc_temp2;
        }
      }                                    // ここまで最大iteration回繰り返す
      if (iter==iteration) fprintf(fplog,"%s%d%s%d\n","Not converged at ",n,"-",i);
                                           // 収束しなかったメッセージを出力
//      else fprintf(fplog,"%d%s%d%s%d\n",n,"-",i,":",iter);

    } // ここまでをdivision回繰り返す

    z_old = z_new ;

// ファイルに出力
    fprintf(fpout,"%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n",
        time, z_new, dis1, d2_old, vel1, vel2-vel1, accr1+z_new, accr2+z_new, f1, f2);

// 最大応答値の更新
    if ( acc1_max < fabs(accr1+z) )   acc1_max = fabs(accr1+z);
    if ( acc2_max < fabs(accr2+z) )   acc2_max = fabs(accr2+z);
    if ( dis1_max < fabs(dis1) )      dis1_max = fabs(dis1);
    if ( dis2_max < fabs(d2_old) )    dis2_max = fabs(d2_old);
    if ( frc1_max < fabs(f1) )        frc1_max = fabs(f1);
    if ( frc2_max < fabs(f2) )      frc2_max = fabs(fset);

  }

  fclose(fpin);         // ファイルを閉じる
  fclose(fpout);
  fclose(fplog);

  sprintf( buffer , "%8.2f" , dis1_max ) ;   // 結果表示画面の準備
  Form2->Label_dis1->Caption = buffer ;
  sprintf( buffer , "%8.2f" , acc1_max ) ;
  Form2->Label_acc1->Caption = buffer ;
  sprintf( buffer , "%8.2f" , frc1_max ) ;
  Form2->Label_frc1->Caption = buffer ;
  sprintf( buffer , "%8.2f" , dis2_max ) ;
  Form2->Label_dis2->Caption = buffer ;
  sprintf( buffer , "%8.2f" , acc2_max ) ;
  Form2->Label_acc2->Caption = buffer ;
  sprintf( buffer , "%8.2f" , frc2_max ) ;
  Form2->Label_frc2->Caption = buffer ;

  ProgressBar1->Position = 0 ;
  Form2->Show();                             // 結果表示画面を出して終了

}
//---------------------------------------------------------------------------
bool TForm1::check_1()
// 入力されていない枠がないかどうかチェック
{
bool yet ;

  yet = false ;
  if ( pier_model->ItemIndex == 0 ) {
    if ( Edit_lin_kp->Text == "" ) yet = true ;
  }
  else if ( pier_model->ItemIndex == 1 ) {
    if ( Edit_clo_dy->Text == "" ) yet = true ;
    if ( Edit_clo_py->Text == "" ) yet = true ;
  }
  else {
    if ( Edit_bil_dy->Text == "" ) yet = true ;
    if ( Edit_bil_py->Text == "" ) yet = true ;
    if ( Edit_bil_k2->Text == "" ) yet = true ;
  }

  if ( bearing_model->ItemIndex == 0 ) {
    if ( Edit_b_lin_k1->Text == "" ) yet = true ;
  }
  else {
    if ( Edit_b_bil_dy->Text == "" ) yet = true ;
    if ( Edit_b_bil_py->Text == "" ) yet = true ;
    if ( Edit_b_bil_k2->Text == "" ) yet = true ;
  }

  if ( Edit_m1->Text == "" ) yet = true ;
  if ( Edit_m2->Text == "" ) yet = true ;
  if ( Edit_h1->Text == "" ) yet = true ;
  if ( Edit_h2->Text == "" ) yet = true ;
  if ( Edit_dt->Text == "" ) yet = true ;
  if ( Edit_number->Text == "" ) yet = true ;
  if ( Edit_beta->Text == "" ) yet = true ;
  if ( Edit_div->Text == "" ) yet = true ;
  if ( Edit_eq->Text == "" ) yet = true ;
  if ( Edit_resp->Text == "" ) yet = true ;
  if ( Edit_log->Text == "" ) yet = true ;
  if ( yet == true ) {
    Application->MessageBox("パラメーターが入力されていない箇所があります。",
      "入力エラー",MB_ICONEXCLAMATION | MB_OK ) ;  }
  return ( yet ) ;
}
//---------------------------------------------------------------------------
void __fastcall TForm1::BitBtn_eqClick(TObject *Sender)
{
// 地震波ファイルの選択
// 選択されたファイルを，ファイル名だけ画面に表示
// フル・パスは，->Hintでヒント枠に記憶

AnsiString moji, fname ;
int j, mojisu ;

  if ( OpenDialog1->Execute() ){
    moji = OpenDialog1->FileName ;
    mojisu = moji.Length() ;
    for ( j=mojisu; j>0; j-- ) {
      if ( IsPathDelimiter( moji,j ) ) break;
      fname.Insert( moji[j], 1 ) ;
    }
    Edit_eq->Text = fname ;
    Edit_eq->Hint = moji.c_str();
  }

}
//---------------------------------------------------------------------------
void __fastcall TForm1::BitBtn_respClick(TObject *Sender)
{
// 結果出力ファイルの選択

AnsiString moji, fname ;
int j, mojisu ;

  if ( SaveDialog1->Execute() ){
    moji = SaveDialog1->FileName ;
    mojisu = moji.Length() ;
    for ( j=mojisu; j>0; j-- ) {
      if ( IsPathDelimiter( moji,j ) ) break;
      fname.Insert( moji[j], 1 ) ;
    }
    Edit_resp->Text = fname ;
    Edit_resp->Hint = moji.c_str();
  }

}
//---------------------------------------------------------------------------
void __fastcall TForm1::BitBtn_logClick(TObject *Sender)
{
// ログファイルの選択

AnsiString moji, fname ;
int j, mojisu ;

  if ( SaveDialog1->Execute() ){
    moji = SaveDialog1->FileName ;
    mojisu = moji.Length() ;
    for ( j=mojisu; j>0; j-- ) {
      if ( IsPathDelimiter( moji,j ) ) break;
      fname.Insert( moji[j], 1 ) ;
    }
    Edit_log->Text = fname ;
    Edit_log->Hint = moji.c_str();
  }

}
//---------------------------------------------------------------------------
double TForm1::bilinear( double x, double x_old, double y,
            double xy, double yy, double k1, double k2 )
// バイリニア型
{
double y_linear, y_plus, y_minus;

  y_linear = k1 * ( x - x_old ) + y;
  y_plus   = k2 * ( x - xy ) + yy;
  y_minus  = k2 * ( x + xy ) - yy;

  if ( y_linear > y_plus ) {
    return(y_plus);
  }
  else if ( y_linear < y_minus ) {
    return(y_minus);
  }
  else {
    return(y_linear);
  }
}
//---------------------------------------------------------------------------
double TForm1::degrading( double x, double x_old, double y, double yy,
            double k, double x_max,double x_min )
// 二次剛性がゼロの剛性劣化(Clough)型
//   x     : このステップの変位
//   x_old : 前ステップの変位     y     : 前ステップの復元力
//   yy    : 降伏耐力             k     : 初期剛性
//   x_max : これまでの最大変位   x_min : これまでの最小変位

{
double x0, y1, y2, y3, y4, y5, y_new;

  x0 = x_old - y / k;
  y1 = k * ( x - x0 );

  if ( x >= x_old ) {
    y2 = ( x - x0 ) * yy / ( x_max - x0 );
    if ( x_max == x_old ) {
      y4 = yy;
    }
    else {
      y4 = y + ( yy - y ) * ( x - x_old ) / ( x_max - x_old );
    }
    if ( (y1>=y2) && (y2>=y4) ){
      if ( y2 < yy ) y_new = y2;
                        else y_new = yy;
    }
    else if ( (y1<=y2) && (y2<=y4) ){
      if ( y2 < yy ) y_new = y2;
                        else y_new = yy;
    }
    else if ( (y2>=y4) && (y4>=y1) ){
      if ( y4 < yy ) y_new = y4;
                        else y_new = yy;
    }
    else if ( (y2<=y4) && (y4<=y1) ){
      if ( y4 < yy ) y_new = y4;
                        else y_new = yy;
    }
    else {
      if ( y1 < yy ) y_new = y1;
                        else y_new = yy;
    }
  }
  else {
    y3 = ( x0 - x ) * yy / ( x_min - x0 );
    if ( x_old == x_min ){
      y5 = -yy;
    }
    else {
      y5 = y + ( yy + y ) * ( x - x_old ) / ( x_old - x_min );
    }
    if ( (y1>=y3) && (y3>=y5) ){
      if ( y3 > -yy ) y_new = y3;
                        else y_new = -yy;
    }
    else if ( (y1<=y3) && (y3<=y5) ){
      if ( y3 > -yy ) y_new = y3;
                        else y_new = -yy;
    }
    else if ( (y3>=y5) && (y5>=y1) ){
      if ( y5 > -yy ) y_new = y5;
                        else y_new = -yy;
    }
    else if ( (y3<=y5) && (y5<=y1) ){
      if ( y5 > -yy ) y_new = y5;
                        else y_new = -yy;
    }
    else {
      if ( y1 > -yy ) y_new = y1;
                        else y_new = -yy;
    }
  }

  return ( y_new );

}

//---------------------------------------------------------------------------
void TForm1::natural( double m1, double m2, double k1, double k2,
     double *t1, double *t2, double *fai1, double *fai2,
     double h1, double h2, double damp[2][2] )

// 固有値解析と歪みエネルギー比例減衰の計算

{
double a, b, c, x, y, z, temp1, temp2, h11, h12, h21, h22, hm1, hm2, det ;
double mat1[2][2], mat2[2][2], mat3[2][2];
const double pi = 3.1415926535897932384626433832795 ;

  a = m1 * m2 ;                         // 2自由度系の固有値解析の開始
  b = m1 * k2 + m2 * k1 + m2 * k2 ;
  c = k1 * k2 ;
  x = sqrt( b * b - 4 * a * c ) ;
  y = ( b - x ) / 2. / a ;              // yは1次の固有円振動数の2乗
  z = ( b + x ) / 2. / a ;              // zは2次の固有円振動数の2乗
  *t1 = 2. * pi / sqrt ( y ) ;          // t1とt2が固有周期
  *t2 = 2. * pi / sqrt ( z ) ;

  temp1=k1*(k2-m2*y)*(k2-m2*y);         // 歪みエネルギー比例減衰の計算開始
  temp2=k2*m2*m2*y*y;
  h11=temp1/(temp1+temp2);              // h11が部材1（橋脚）の減衰定数に対する係数
  h12=temp2/(temp1+temp2);              // h11が部材2（支承）の減衰定数に対する係数
  hm1=2.*(h11*h1+h12*h2)*sqrt(y)*m1;    // 2×1次の減衰定数×ω1（対角化された減衰行列の要素）
  temp1=k1*(k2-m2*z)*(k2-m2*z);
  temp2=k2*m2*m2*z*z;
  h21=temp1/(temp1+temp2);
  h22=temp2/(temp1+temp2);
  hm2=2.*(h21*h1+h22*h2)*sqrt(z)*m2;    // 2×2次の減衰定数×ω2
  
  temp1=(k2-m2*y)/k2;
  *fai1=temp1;                          // 1次固有モードベクトル（部材2を1とする）
  temp2=(k2-m2*z)/k2;
  *fai2=temp2;                          // 2次固有モードベクトル

  mat1[0][0]=temp1; mat1[0][1]=temp2; mat1[1][0]=1.; mat1[1][1]=1.; // 減衰行列の計算
  mat2[0][0]=hm1;  mat2[0][1]=0.;   mat2[1][0]=0.; mat2[1][1]=hm2;
  matmult( mat1, mat2, mat3 );
  det=temp1-temp2;
  mat1[0][0]=1./det; mat1[0][1]=-temp2/det;
  mat1[1][0]=-1./det; mat1[1][1]=temp1/det;
  matmult( mat3, mat1, damp );              // モーダル行列×対角化減衰行列×モーダル行列の逆行列

}
//---------------------------------------------------------------------------
void TForm1::matmult( double mat1[2][2], double mat2[2][2],
                                 double mat3[2][2])
// 行列どうしのかけ算をする関数
// mat3=mat1×mat2
{
  int i,j,k;
  double temp;
  for (i=0;i<2;i++) {
    for (j=0;j<2;j++) {
      temp=0.;
      for (k=0;k<2;k++) temp += mat1[i][k]*mat2[k][j];
      mat3[i][j]=temp;
    }
  }
}
//---------------------------------------------------------------------------
void __fastcall TForm1::pier_modelChange(TObject *Sender)

// 橋脚の履歴復元力特性に合わせて入力パネルの表示を変更
{
  switch ( pier_model->ItemIndex ) {
    case 0:
      Panel_bilinear->Visible = false;
      Panel_clough->Visible = false;
      Panel_linear->Visible = true;
      break;
    case 1:
      Panel_bilinear->Visible = false;
      Panel_linear->Visible = false;
      Panel_clough->Visible = true;
      break;
    case 2:
      Panel_linear->Visible = false;
      Panel_clough->Visible = false;
      Panel_bilinear->Visible = true;
  }
}
//---------------------------------------------------------------------------
void __fastcall TForm1::bearing_modelChange(TObject *Sender)

// 支承の履歴復元力特性に合わせて入力パネルの表示を変更
{
  switch ( bearing_model->ItemIndex ) {
    case 0:
      Panel_b_bilinear->Visible = false;
      Panel_b_linear->Visible = true;
      break;
    case 1:
      Panel_b_linear->Visible = false;
      Panel_b_bilinear->Visible = true;
  }
}
//---------------------------------------------------------------------------
void __fastcall TForm1::BitBtn_helpClick(TObject *Sender)

// 解説画面の表示
{
  Form3->Show();
}
//---------------------------------------------------------------------------