// slidespec.cpp
// 2003.2.22 伊津野和行
//
// 2自由度系の非線形応答スペクトル
// すべり摩擦支承版
// Borland C++ Builder 5 を使用
//
//---------------------------------------------------------------------------
#include <vcl\vcl.h>
#include <stdio.h>
#include <math.h>
#pragma hdrstop

#include "spec.h"
#include "about.h"
//---------------------------------------------------------------------------
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
	: TForm(Owner)
{
}
//---------------------------------------------------------------------------
void __fastcall TForm1::bt_execClick(TObject *Sender)
{
float kp, k1, m1, m2, h, dt, a, b, d, z_old, z_new,
	ainv, binv, dinv, dis1, dis2,
	vel1, vel2, accr1, accr2,
	f1, f2, z, dt2, d1_old, d2_old,
	adbc, damp, t1, t2, x1, x2,
	p_max, p_min, f1_old, f2_old,
	k2, p_dy, p_py, b_dy, b_py,fset,
	acc1_max, acc2_max, vel1_max,
	vel2_max, dis1_max, dis2_max,
	frc1_max, frc2_max, masatsu, k2_max ;
int n, number, total, i, hys ;
bool shoki ;
const float g = 9.8 ;

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

// ファイルのオープン

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

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

	shoki = CheckBox1->Checked ;
	if ( Rb_rc->Checked ) {
		h = 0.02 ;			  // ＲＣ橋脚は減衰定数2％
		hys = 0 ;
	} else if ( Rb_steel->Checked ) {
		h = 0.01 ;			  // 鋼製橋脚は減衰定数1％
		hys = 1 ;
	} else {
		h = StrToFloat( Box_h->Text ) ;
		hys = Box_model->ItemIndex ;
	}
	m1 = StrToFloat( edit_m1->Text ) ;		// 橋脚質量
	m2 = StrToFloat( edit_m2->Text ) ;		// 桁質量
	p_dy = StrToFloat( edit_pier_dy->Text )/100.;	// 橋脚降伏変位 cm->m
	p_py = StrToFloat( edit_pier_py->Text )*1000.;	// 橋脚降伏荷重 MN->kN
	dt = StrToFloat( edit_dt->Text );		// 時間刻み
	total = StrToInt( edit_number->Text );		// 地震波データ数
	if ( shoki ) {
		k1 = StrToFloat( edit_k1->Text )*1000;	// 支承の初期剛性 MN/m->kN/m
	}
	else {
		k1 = 100.;
	}
	k2_max = StrToFloat( edit_k2->Text )*1000.;     // ２次剛性の最大値 MN/m->kN/m
	masatsu = StrToFloat( edit_myu->Text );		// 摩擦係数
	kp = p_py / p_dy ;				// 橋脚の初期剛性

	damp = 2. * h * sqrt( (m1+m2) * kp );	// 粘性減衰係数
	dt = dt / 10. ;							// 10回にわけて計算
	dt2 = dt * dt;

//  計算途中で使う行列の逆行列を作成

	a = m1 + dt * damp / 2. + dt2 * (kp+k1) / 4. ;
	b = -dt2 * k1 / 4.;
	d = m2 + dt2 * k1 / 4.;
	adbc = a * d - b * b;
	ainv =  d / adbc;
	binv = -b / adbc;
	dinv =  a / adbc;

// 出力用ファイルに各列の説明を出力、nは時間ステップ
	fprintf(fpout,"%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n",
		"２次剛性(MN/m)","橋脚変位(cm)","支承変位(cm)","橋脚速度(m/s)",
		"桁速度(m/s)","橋脚加速度(gal)","桁加速度(gal)",
		"橋脚復元力(MN)","支承復元力(MN)","橋脚じん性率");

	k2 = 0. ;
	for ( number = 1 ; number <= 20 ; number++ ) {

		k2 = k2 + k2_max / 20. ;
		b_dy = m2 * masatsu * g / ( k1 - k2 ) ;		// 支承の降伏変位
		if ( shoki ) {								// 支承の降伏荷重
			b_py = k1 * b_dy ;
		} else {
			b_py = m2 * masatsu * g ;
		}
		dis1 = dis2 = vel1 = vel2 = accr1 = accr2 = 0. ;
		d1_old = d2_old = f1_old = f2_old = 0. ;
		acc1_max = acc2_max = vel1_max = vel2_max = 0. ;
		dis1_max = dis2_max = frc1_max = frc2_max = 0. ;
		p_min = -p_dy;							// 最小橋脚変位
		p_max =  p_dy;							// 最大橋脚変位
		n = 0 ;

// 入力地震波データがある限り繰り返す
		z = 0. ;
		z_old = 0. ;
		while( fscanf(fpin,"%f",&z_new ) != EOF ){
			z_new /= 100. ;
			n = n + 1 ;
			if ( n > total ) {
				break ;
			}

// Operator Splitting法による数値積分  10回にわけて計算

			for ( i=1 ; i<=10 ; i++ ) {
				z = z_old + ( z_new - z_old ) / 10. * i ;
				dis1 = dis1 + dt * vel1 + dt2 / 4. * accr1;
				dis2 = dis2 + dt * vel2 + dt2 / 4. * accr2;

// ＲＣ橋脚は剛性劣化(Clough)型、鋼製橋脚はバイリニア
				if ( hys == 0 ) {
					f1 = degrading( dis1, d1_old, f1_old, p_py, kp,
						p_max, p_min ) ;
				} else {
					f1 = bilinear( dis1, d1_old, f1_old, p_dy,
						p_py, kp, kp/100. ) ;
				}
// 支承はバイリニア、初期剛性が設定されない場合剛塑性バイリニア
				if ( shoki ) {
					f2 = bilinear( dis2-dis1, d2_old, f2_old, b_dy,
						b_py, k1, k2 ) ;
				} else {
					f2 = r_bilinear( dis2-dis1, d2_old, f2_old, b_py, k2 ) ;
				}

				t1 = -m1 * z - f1 + f2 - damp * ( vel1 + dt * accr1 / 2. );
				t2 = -m2 * z - f2 ;

				x1 = ainv * t1 + binv * t2 ;
				x2 = binv * t1 + dinv * t2 ;

				dis1 = dis1 + dt2 / 4 * x1;
				dis2 = dis2 + dt2 / 4 * x2;
				vel1 = vel1 + dt / 2. * ( accr1 + x1 );
				vel2 = vel2 + dt / 2. * ( accr2 + x2 );
				accr1 = x1;
				accr2 = x2;

				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;
				}

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

			fset = masatsu * m2 * g + ( f2 - masatsu * m2 * g ) * 2. ;
			z_old = z_new ;

// 最大応答値の更新
			acc1_max = max( (float)fabs(accr1+z) , acc1_max );
			acc2_max = max( (float)fabs(accr2+z) , acc2_max );
			vel1_max = max( (float)fabs(vel1) , vel1_max );
			vel2_max = max( (float)fabs(vel2-vel1) , vel2_max );
			dis1_max = max( (float)fabs(dis1) , dis1_max );
			dis2_max = max( (float)fabs(d2_old) , dis2_max );
			frc1_max = max( (float)fabs(f1), frc1_max ) ;
			frc2_max = max( (float)fabs(fset), frc2_max ) ;
		}

// ファイルに出力
		fprintf( fpout, "%f \t %f \t %f \t %f \t %f \t %f"
						"\t %f \t %f \t %f \t %f \n",
			k2 / 1000., dis1_max*100., dis2_max*100., vel1_max, vel2_max,
			acc1_max*100., acc2_max*100., frc1_max/1000., frc2_max/1000.,
			dis1_max / p_dy ) ;

		ProgressBar1->StepIt() ;
		rewind ( fpin ) ;
	}

	fclose(fpin);
	fclose(fpout);
	ProgressBar1->Position = 0 ;
	Application->MessageBox( "計算終了", "スペクトル計算", MB_OK ) ;

}
//---------------------------------------------------------------------------
bool TForm1::check_1()
{
bool yet ;

	yet = false ;
	if ( edit_m1->Text == "" ) { yet = true ;}
	if ( edit_m2->Text == "" ) { yet = true ;}
	if ( edit_pier_dy->Text == "" ) { yet = true ;}
	if ( edit_pier_py->Text == "" ) { yet = true ;}
	if ( edit_myu->Text == "" ) { yet = true ;}
	if ( (CheckBox1->Checked) && (edit_k1->Text == "") ) {
		yet = true ;
	}
	if ( edit_k2->Text == "" ) { yet = true ;}
	if ( edit_dt->Text == "" ) { yet = true ;}
	if ( edit_number->Text == "" ) { yet = true ;}
	if ( edit_eq->Text == "" ) { yet = true ;}
	if ( edit_output->Text == "" ) { yet = true ;}
	if ( yet == true ) {
		error_1() ;
	}
	return ( yet ) ;
}
//---------------------------------------------------------------------------
void TForm1::error_1()
{
	Application->MessageBox("パラメーターが入力されていない箇所があります。",
		"入力エラー",MB_ICONEXCLAMATION | MB_OK ) ;
}
//---------------------------------------------------------------------------
void __fastcall TForm1::bt_eqClick(TObject *Sender)

// 地震波ファイルの選択
{
	if ( OpenDialog1->Execute() ){
		edit_eq->Text = OpenDialog1->FileName;
	}
}
//---------------------------------------------------------------------------
void __fastcall TForm1::bt_outClick(TObject *Sender)

// 結果出力ファイルの選択
{
	if ( SaveDialog1->Execute() ){
		edit_output->Text = SaveDialog1->FileName;
	}
}
//---------------------------------------------------------------------------
void __fastcall TForm1::BitBtn2Click(TObject *Sender)

// 説明画面の表示
{
	Form2->Show() ;
}
//---------------------------------------------------------------------------
void __fastcall TForm1::CheckBox1Click(TObject *Sender)

// 初期剛性を設定するかどうか
{
	if ( CheckBox1->Checked ) {
		Label_k1->Enabled = true ;
		edit_k1->Enabled = true ;
		edit_k1->Visible = true ;
	}
	else {
		Label_k1->Enabled = false ;
		edit_k1->Enabled = false ;
		edit_k1->Visible = false ;
	}
}
//---------------------------------------------------------------------------
float TForm1::bilinear(float x,float x_old,float y,
						float xy,float yy,float k1,float k2)
// バイリニア型
{
float 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);
	}
}
//---------------------------------------------------------------------------
float TForm1::r_bilinear( float x, float x_old, float y,
										float yy, float k )
// 剛塑性の骨格曲線をもったバイリニア
{
	if ( x == x_old ) {
		return ( y ) ;
	}
	else if ( x > x_old ) {
		return ( k * x + yy ) ;
	}
	else {
		return ( k * x - yy ) ;
	}
}

//---------------------------------------------------------------------------
float TForm1::degrading(float x,float x_old,float y, float yy,
						float k, float x_max,float x_min)
// 剛性劣化(Clough)型
{
float 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) ){
			y_new = min( y2, yy );
		}
		else if ( (y1<=y2) && (y2<=y4) ){
			y_new = min( y2, yy );
		}
		else if ( (y2>=y4) && (y4>=y1) ){
			y_new = min( y4, yy );
		}
		else if ( (y2<=y4) && (y4<=y1) ){
			y_new = min( y4, yy );
		}
		else {
			y_new = min( y1, 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) ){
			y_new = max( y3, -yy );
		}
		else if ( (y1<=y3) && (y3<=y5) ){
			y_new = max( y3, -yy );
		}
		else if ( (y3>=y5) && (y5>=y1) ){
			y_new = max( y5, -yy );
		}
		else if ( (y3<=y5) && (y5<=y1) ){
			y_new = max( y5, -yy );
		}
		else {
			y_new = max( y1, -yy );
		}
	}

	return ( y_new );

}
//---------------------------------------------------------------------------
void __fastcall TForm1::FormCreate(TObject *Sender)
{
//	if ( Screen->Width <= 800 ) {
//		Form1->ClientHeight = 400 ;
//		Form1->ClientWidth  = 410 ;
//	} else {
//		Form1->ClientHeight = 510 ;
//		Form1->ClientWidth  = 520 ;
//	}
	Label_h->Visible = false ;
	Label_model->Visible = false ;
	Box_h->Visible = false ;
	Box_model->Visible = false ;
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Rb_etcClick(TObject *Sender)
{
	Label_h->Visible = true ;
	Label_model->Visible = true ;
	Box_h->Visible = true ;
	Box_model->Visible = true ;
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Rb_rcClick(TObject *Sender)
{
	Label_h->Visible = false ;
	Label_model->Visible = false ;
	Box_h->Visible = false ;
	Box_model->Visible = false ;
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Rb_steelClick(TObject *Sender)
{
	Label_h->Visible = false ;
	Label_model->Visible = false ;
	Box_h->Visible = false ;
	Box_model->Visible = false ;
}
//---------------------------------------------------------------------------