//---------------------------------------------------------------------------
#include <vcl.h>
#pragma hdrstop

#include "main.h"
#include "input.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
    : TForm(Owner)
{
}
//---------------------------------------------------------------------------

void __fastcall TForm1::BitBtn_inputClick(TObject *Sender)
{
     Form2->Show() ;
}
//---------------------------------------------------------------------------
void __fastcall TForm1::BitBtn_firstClick(TObject *Sender)
{
double a, ti, kakudo ;
float z ;
int i, n, total, jiban, type, hajime ;
const double PI = 3.14159265358979 ;

    StatusBar1->SimpleText = "計算中" ;
    StatusBar1->Refresh() ;

    // ファイルのオープン

    if (( fpin=fopen(AnsiString(Form2->Edit_input->Text).c_str(),"rt")) == NULL ){
		Application->MessageBox(L"加速度波形ファイルを設定し直して下さい。",
		L"入力ファイルエラー",MB_ICONEXCLAMATION | MB_OK ) ;
        return ;
    }

    dt  = StrToFloat( Form2->Edit_dt->Text );      // 時間刻み
    total = StrToInt( Form2->Edit_number->Text );  // 地震波個数
    h = StrToFloat( Form2->Edit_gensui->Text );      // 構造物の減衰定数
    type  = Form2->CSpinEdit_type->Value ;         // 地震波タイプ
    jiban = Form2->CSpinEdit_jiban->Value ;        // 地盤種別

    a = log( (double)(total-1) ) / log( 2. ) + 1.0 ;
    nn = pow ( 2., (int) a ) ;

    try {                      // 例外のテスト
      accr = new double[nn];
      time2 = new double[nn];
      freq = new double[nn];
      r_freq = new double[nn];
      sa = new double[nn];
      res = new double[nn];
      comp = new struct complex<double>[nn];
    }
    catch (bad_alloc) { // bad_alloc が送出された場合にのみこのブロックに入る
        Application->MessageBox
		(L"メモリが足りません。データ個数を減らして下さい。",
		 L"エラー",MB_ICONEXCLAMATION | MB_OK ) ;
        exit (-1);
    }

      TCursor Save_Cursor = Screen->Cursor;

      Screen->Cursor = crHourGlass;    // 砂時計カーソルを表示する

      try
      {
        // 砂時計カーソルを表示する
      }
      catch (...)
      {
            Screen->Cursor = Save_Cursor; // 必ず通常のカーソルに戻る
            throw;
      }

    i = 0;
    while( fscanf(fpin,"%f",&z ) != EOF ){
        accr[i] = (double) z;
        i++ ;
        if ( i >= total ) break;
    }
    fclose(fpin);
    n = i;

    for ( i=0; i<nn; i++ ) {
        time2[i] = dt * (double)i ;
        freq[i] = (double)i / (double)nn / dt ;
    }

    senkai( time2, accr, n ) ;
    hajime = n / 10 ;
    for ( i=0; i<hajime; i++ ) {
        kakudo = PI / 2. / (double)hajime * (double)i ;
        accr[i] = accr[i] * ( 1. - cos( kakudo ) ) ;
        accr[n-i-1] = accr[n-i-1] * ( 1. - cos( kakudo ) ) ;
    }

    for ( i=n; i<nn; i++ ) accr[i] = 0. ;

    kpoint = 0 ;
    for ( i=0; i<nn; i++ ) {
        if ( freq[i] > 14.5 ) break ;
        if ( freq[i] > 0.25 ) {
            r_freq[kpoint] = freq[i] ;
            ti = 1. / freq[i] ;
            sa[kpoint] = sekkei( ti, type, jiban ) ;
            kpoint++ ;
        }
    }

                                // 応答スペクトルの計算
    resspec( r_freq, accr, res, dt, h, nn, kpoint ) ;

    if ( Series1->FirstValueIndex != -1 ) {
        Series1->Clear() ;
        Series2->Clear() ;
    }

    for ( i=0; i<kpoint; i++ ) {
        ti = 1. / r_freq[i] ;
        Series1->AddXY( ti,res[i],"",clRed);
        Series2->AddXY( ti,sa[i], "",clGreen);
    }

    StatusBar1->SimpleText = "計算終了" ;
    Screen->Cursor = Save_Cursor;                              // フーリエ変換で加速度のフーリエスペクトル

}
//---------------------------------------------------------------------------
// 設計スペクトル

double TForm1::sekkei( double ti, int type, int jiban )
{
double sa ;

    if ( type == 1 ) {
        switch ( jiban ) {
            case 1 :
                if ( ti <= 1.4 ) sa = 700. ;
                else sa = 980 / ti ;
                break ;
            case 2 :
                if ( ti < 0.18 ) {
                    sa = 1505. * pow( ti, 1./3. ) ;
                    if ( sa < 700. ) sa = 700. ;
                }
                else {
                    if ( ti > 1.6 ) sa = 1360. / ti ;
                    else sa = 850. ;
                }
                break ;
            case 3 :
                if ( ti < 0.29 ) {
                    sa = 1511. * pow( ti, 1./3. ) ;
                    if ( sa < 700. ) sa = 700. ;
                }
                else {
                    if ( ti > 2.0 ) sa = 2000. / ti ;
                    else sa = 1000. ;
                }
        }
    }
    else if ( type == 2 ) {
        switch ( jiban ) {
            case 1 :
                if ( ti < 0.3 ) sa = 4463. * pow( ti, 2./3. ) ;
                else {
                    if ( ti > 0.7 ) sa = 1104. / pow( ti, 5./3. ) ;
                    else sa = 2000. ;
                }
                break ;
            case 2 :
                if ( ti < 0.4 ) sa = 3224. * pow( ti, 2./3. ) ;
                else {
                    if ( ti > 1.2 ) sa = 2371. / pow( ti, 5./3. ) ;
                    else sa = 1750. ;
                }
                break ;
            case 3 :
                if ( ti < 0.5 ) sa = 2381. * pow( ti, 2./3. ) ;
                else {
                    if ( ti > 1.5 ) sa = 2948. / pow( ti, 5./3. ) ;
                    else sa = 1500. ;
                }
        }
    }
    else {
        switch ( jiban ) {
            case 1 :
                if ( ti < 0.2 ) sa = 17600. * pow( ti, 2./3. ) ;
                else {
                    if ( ti > 0.5 ) sa = 1890. / pow( ti, 5./3. ) ;
                    else sa = 6000. ;
                }
                break ;
            case 2 :
                if ( ti < 0.4 ) sa = 9300. * pow( ti, 2./3. ) ;
                else {
                    if ( ti > 0.8 ) sa = 3450. / pow( ti, 5./3. ) ;
                    else sa = 5000. ;
                }
                break ;
            case 3 :
                if ( ti < 0.5 ) sa = 4750. * pow( ti, 2./3. ) ;
                else {
                    if ( ti > 1.5 ) sa = 5900. / pow( ti, 5./3. ) ;
                    else sa = 3000. ;
                }
        }
    }
    return ( sa ) ;
}
//---------------------------------------------------------------------------
// 高速フーリエ変換
//      参考文献：大崎順彦「新・地震動のスペクトル解析入門」鹿島出版会

void TForm1::fast ( struct complex<double> x[], int nn, int ind )
{
struct complex<double> temp, theta;
int i, j, k, m, kmax, istep;
const double PI = 3.14159265358979 ;

    j = 0 ;
    for ( i=0; i<nn; i++ ) {
        if ( i<j ) {
            temp = x[j] ;
            x[j] = x[i] ;
            x[i] = temp ;
        }
        m = nn / 2 ;
        while ( j+1>m ) {
            j = j - m ;
            m = m / 2 ;
            if ( m<2 ) break;
        }
        j = j + m ;
    }
    kmax = 1 ;
    while ( kmax<nn ) {
        istep = kmax * 2 ;
        for ( k=0; k<kmax; k++ ) {
            theta = complex<double>( 0., PI * ind * k / kmax ) ;
            for ( i=k; i<nn; i=i+istep ) {
                j = i + kmax ;
                temp = x[j] * exp( theta ) ;
                x[j] = x[i] - temp ;
                x[i] = x[i] + temp ;
            }
        }
        kmax = istep ;
    }
}
//---------------------------------------------------------------------------
// 基線補正 線形回帰

void TForm1::senkai ( double x[], double y[], int n )
{
double ax=0., ay=0., b=0., c=0., alph, beta;
int i;

    for ( i=0; i<n; i++ ) {
        ax += x[i] ;
        ay += y[i] ;
    }
    ax /= n ;
    ay /= n ;

    for ( i=0; i<n; i++ ) {
        b += ( x[i] - ax ) * ( y[i] - ay ) ;
        c += ( x[i] - ax ) * ( x[i] - ax ) ;
    }
    beta = b / c ;
    alph = ay - beta * ax ;
    for ( i=0; i<n; i++ )
        y[i] = y[i] - alph - beta * x[i] ;
}
//---------------------------------------------------------------------------
void TForm1::resspec ( double freq[], double z[],
        double res[], double dt, double h, int total, int kpoint )
{

double ainv, dis, vel, acc, dt2, omega,
       c, k, x, acc_max ;
const double PI = 3.14159265358979 ;
int i, j ;

	dt2 = dt * dt;

    for ( i=0 ; i<kpoint ; i++ ) {
    	omega = 2. * PI * freq[i];
        k = omega * omega ;
        c = 2. * h * omega ;
        acc = 0.;
        vel = 0.;
        dis = 0.;
        acc_max = 0.;

//  計算途中で使う係数を計算

		ainv = 1. + dt * c / 2. + dt2 * k / 4. ;

// Operator Splitting法による数値積分

        for ( j=0; j<total; j++ ) {
			dis = dis + dt * vel + dt2 / 4. * acc;
        	x = ( -z[j] - k * dis - c * ( vel + dt * acc / 2. )) / ainv;

        	dis = dis + dt2 / 4 * x;
        	vel = vel + dt / 2. * ( acc + x );
        	acc = x;

// 最大応答値の更新

            if ( fabs(acc+z[j]) > acc_max ) acc_max = fabs(acc+z[j]);

    	}

        res[i] = acc_max;
	}
}


//---------------------------------------------------------------------------
void __fastcall TForm1::BitBtn_modifyClick(TObject *Sender)
{
int i, k, hajime, nfold ;
double ck, kakudo, ti, fll, flu, ful, fuu ;
const complex<double> zero(0.,0.);
const double PI = 3.14159265358979 ;

      StatusBar1->SimpleText = "計算中" ;
      StatusBar1->Refresh() ;
      TCursor Save_Cursor = Screen->Cursor;

      Screen->Cursor = crHourGlass;    // 砂時計カーソルを表示する

      try
      {
        // 砂時計カーソルを表示する
      }
      catch (...)
      {
            Screen->Cursor = Save_Cursor; // 必ず通常のカーソルに戻る
            throw;
      }

        for ( i=0; i<nn; i++ )
            comp[i] = complex<double>( accr[i], 0. ) ;
        fast  ( comp, nn, -1 ) ;
        for ( i=0; i<nn; i++ )
            comp[i] /= (double)nn ;

                                // フーリエスペクトルの修正
        k = 0 ;
        nfold = nn / 2 ;
        comp[0] = zero ;
        for ( i=1; i<nfold; i++ ) {
            if ( freq[i] > 14.5 ) break ;
            if ( freq[i] > 0.25 ) {
                ck = sa[k] / res[k] ;
                comp[i] = comp[i] * ck ;
                comp[nn-i] = comp[nn-i] * ck ;
                k++ ;
            }
        }
        for ( k=i; k<nfold; k++ ) {
            comp[k] = comp[k] * ck ;
            comp[nn-k] = comp[nn-k] * ck ;
        }
        fll = 1. / 10. ;
        flu = 1. / 4. ;
        ful = 14.5 ;
        fuu = 20. ;
        filter ( comp, dt, fll, flu, ful, fuu, nn ) ;
                                // フーリエ逆変換で加速度に
        fast ( comp, nn, +1 ) ;
        for ( i=0; i<nn; i++ )
            accr[i] = real( comp[i] ) ;

        hajime = nn / 10 ;
        for ( i=0; i<hajime; i++ ) {
            kakudo = PI / 2. / (double)hajime * (double)i ;
            accr[i] = accr[i] * ( 1. - cos( kakudo ) ) ;
            accr[nn-i-1] = accr[nn-i-1] * ( 1. - cos( kakudo ) ) ;
        }
                                // 応答スペクトルの計算
        resspec( r_freq, accr, res, dt, h, nn, kpoint ) ;

        for ( i=kpoint-1; i>=0; i-- )
            Series1->Delete( i ) ;

        for ( i=0; i<kpoint; i++ ) {
            ti = 1. / r_freq[i] ;
            Series1->AddXY( ti,res[i],"",clRed);
        }

        saveflag = 1 ;
        StatusBar1->SimpleText = "計算終了" ;
        Screen->Cursor = Save_Cursor;

}
//---------------------------------------------------------------------------

void __fastcall TForm1::BitBtn_saveClick(TObject *Sender)
{
int i;
double ti;

      StatusBar1->SimpleText = "保存中" ;
      StatusBar1->Refresh() ;
      TCursor Save_Cursor = Screen->Cursor;

      Screen->Cursor = crHourGlass;    // 砂時計カーソルを表示する

      try
      {
        // 砂時計カーソルを表示する
      }
      catch (...)
      {
            Screen->Cursor = Save_Cursor; // 必ず通常のカーソルに戻る
            throw;
      }

// ファイルに出力

	if (( fpout=fopen(AnsiString(Form2->Edit_output->Text).c_str(),"wt")) == NULL ){
		Application->MessageBox(L"結果の出力ファイルを設定し直して下さい。",
		L"出力ファイルエラー",MB_ICONEXCLAMATION | MB_OK ) ;
        return ;
    }
    fprintf(fpout,"%s\t%s\n",
        "時間(sec)", "加速度(gal)" );
    for ( i=0; i<nn; i++ )
        fprintf(fpout,"%lf\t%lf\n",
            time2[i], accr[i] );

    fclose(fpout);

// スペクトルの保存

    if (( fpout=fopen(AnsiString(Form2->Edit_spec->Text).c_str(),"wt")) != NULL ){
        fprintf(fpout, "%s\t%s\t%s\n",
            "周期(sec)","応答スペクトル(gal)","設計スペクトル(gal)" );
        for ( i=kpoint-1; i>=0; i-- ) {
            ti = 1. / r_freq[i] ;
            fprintf(fpout,"%lf\t%lf\t%lf\n",
                ti, res[i], sa[i] );
        }
        fclose(fpout) ;
    }

    saveflag = 0 ;
    Screen->Cursor = Save_Cursor;                              // フーリエ変換で加速度のフーリエスペクトル
    StatusBar1->SimpleText = "データ保存終了" ;

}
//---------------------------------------------------------------------------

void __fastcall TForm1::BitBtn_quitClick(TObject *Sender)
{
    if ( saveflag == 1 ) {
        int RetValue = Application->MessageBox(
		L"データが保存されていません。\nこのまま終了しますか？",
		L"確認", MB_ICONEXCLAMATION | MB_YESNO ) ;
        if ( RetValue == IDNO ) return ;
    }
        delete[] time2 ;
        delete[] freq ;
        delete[] accr ;
        delete[] r_freq ;
        delete[] comp ;
        delete[] sa ;
        delete[] res ;
    Close() ;

}
//---------------------------------------------------------------------------

void __fastcall TForm1::BitBtn1Click(TObject *Sender)
{
     Application->HelpFile = "modify.hlp";
     Application->HelpCommand(HELP_FINDER, 0);
}
//---------------------------------------------------------------------------

void __fastcall TForm1::FormCreate(TObject *Sender)
{
    	ClientHeight = 335 ;
    	ClientWidth  = 570 ;
        saveflag = 0 ;
}
//---------------------------------------------------------------------------
// バンドパスフィルタ　（台形フィルタ）

void TForm1::filter ( struct complex<double> x[], double dt, double fll, double flu,
    double ful, double fuu, int nn )
{
const struct complex<double> zero(0.,0.);
double xll, xlu, xul, xuu, factor;
int kll, klu, kul, kuu, i, kll1, kuu1;

	xll = nn * dt * fll ;
	xlu = nn * dt * flu ;
	xul = nn * dt * ful ;
	xuu = nn * dt * fuu ;
	kll = Ceil(xll) ;
	klu = Ceil(xlu) ;
	kul = Ceil(xul) ;
	kuu = Ceil(xuu) ;
	if ( kll == klu ) kll = klu - 1 ;
	if ( kul == kuu ) kuu = kul + 1 ;
	kll1 = nn - kll ;
	kuu1 = nn - kuu ;
	for ( i=0; i<kll; i++ ) x[i] = zero ;
	for ( i=kll1+1; i<nn; i++) x[i] = zero ;
	for ( i=kll; i<klu; i++ ) {
		factor = ((double)i/(double)nn/dt-fll) / (flu-fll) ;
		x[i] *= factor;
		x[nn-i] *= factor;
	}
	for ( i=kul+1; i<kuu; i++ ) {
		factor = (fuu-(double)i/(double)nn/dt) / (fuu-ful) ;
		x[i] *= factor;
		x[nn-i] *= factor;
	}
	for ( i=kuu; i<=kuu1; i++ ) x[i] = zero ;
}
//---------------------------------------------------------------------------