//---------------------------------------------------------------------------
// Ver.1.1    '99.2.10
//     基線補正をはずすオプション
// Ver.1.02   '98.2.14
//     infr の不備を修正
// Ver.2.0    '99.9.9
//     ローパスとハイパスのフィルターを追加
// Ver.2.1    '99.9.15
//     スペクトル値の間違いを修正
// Ver.3.0    '07.2.4
//     台形フィルターにコサイン型の漸減曲線を追加
// Ver.3.1    '09.2.18
//     フィルターの説明図を改良

#include <vcl.h>
#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <except.h>

#pragma hdrstop

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

void __fastcall TForm1::bt_execClick(TObject *Sender)
{
float z;
double fll, flu, ful, fuu, dt, a,
    *time, *freq, *accr, *velr, *disr,
    *afft, *vfft, *dfft ;
double acc_max, vel_max, dis_max ;
struct complex<double> *comp, *temp ;
const struct complex<double> zero(0.,0.);
const double PI = acos(-1.0) ;
int n, total, i, nn ;
char buffer[8] ;

    TCursor Save_Cursor = Screen->Cursor;

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

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

    Form2->Hide() ;
    Form3->Hide() ;
    StatusBar1->SimpleText = "計算開始" ;
    StatusBar1->Refresh() ;

// パラメーターが全部入力されているか
    if ( check_1() ) {
        StatusBar1->SimpleText = "" ;
        StatusBar1->Refresh() ;
        Screen->Cursor = Save_Cursor;
        return ;
    }

// ファイルのオープン

	if (( fpin=_wfopen(OpenDialog1->FileName.c_str(),L"rt")) == NULL ){
		Application->MessageBox(L"加速度波形ファイルを設定し直して下さい。",
		L"入力ファイルエラー",MB_ICONEXCLAMATION | MB_OK ) ;
        StatusBar1->SimpleText = "" ;
        StatusBar1->Refresh() ;
        Screen->Cursor = Save_Cursor;
        return ;
    }

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

    fll = 0. ;  flu = 0. ;  flu = 0. ; fuu = 0. ;
    if ( RadioButton1->Checked ) {
        fll = StrToFloat( edit_fll->Text ) ;
        flu = StrToFloat( edit_flu->Text ) ;
        ful = StrToFloat( edit_ful->Text ) ;
        fuu = StrToFloat( edit_fuu->Text ) ;
    } else if ( RadioButton2->Checked ) {
        fll = 1. / StrToFloat( edit_fll->Text ) ;
        flu = 1. / StrToFloat( edit_flu->Text ) ;
    } else {
        ful = StrToFloat( edit_ful->Text ) ;
        fuu = StrToFloat( edit_fuu->Text ) ;
    }
    dt = StrToFloat( edit_dt->Text );         // 時間刻み
    total = StrToInt( edit_number->Text );    // 地震波個数

    try {    // 例外のテスト
      accr = new double[total];
    }
    catch (bad_alloc) {
		Application->MessageBox(L"入力地震波を記憶するメモリが足りません。終了します。",
		L"エラー",MB_ICONEXCLAMATION | MB_OK ) ;
        StatusBar1->SimpleText = "" ;
        StatusBar1->Refresh() ;
        Screen->Cursor = Save_Cursor;
        exit (-1);
    }

    acc_max = 0. ;
    vel_max = 0. ;
    dis_max = 0. ;
    i = 0;
    while( fscanf(fpin,"%f",&z ) != EOF ){
        accr[i] = (double) z;
        i++ ;
        if ( i >= total ) break;
    }
    fclose(fpin);
    n = i;
    a = log( (double)(n-1) ) / log( 2. ) + 1.0 ;
    nn = std::pow ( 2., (int) a ) ;

    try {                      // 例外のテスト
        time = new double[nn];
        freq = new double[nn];
        velr = new double[nn];
        disr = new double[nn];
        afft = new double[nn];
        vfft = new double[nn];
        dfft = new double[nn];
        comp = new struct complex<double>[nn];
        temp = new struct complex<double>[nn];
    }
    catch (bad_alloc) { // bad_alloc が送出された場合にのみこのブロックに入る
		Application->MessageBox(L"メモリが足りません。終了します。",
		L"エラー",MB_ICONEXCLAMATION | MB_OK ) ;
        StatusBar1->SimpleText = "" ;
        StatusBar1->Refresh() ;
        Screen->Cursor = Save_Cursor;
        exit (-1);
    }

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

                              // 線形回帰によって基線補正
    if ( CheckBox_hosei->Checked )
        senkai( time, accr, n ) ;

                              // FFTの準備
    for ( i=0; i<n; i++ )
        comp[i] = complex<double>( accr[i], 0. ) ;
    for ( i=n; i<nn; i++ )
        comp[i] = zero ;
                              // フーリエ変換で加速度のフーリエスペクトル
    StatusBar1->SimpleText = "フーリエ変換" ;
    StatusBar1->Refresh() ;
    fast  ( comp, nn, -1 ) ;
    for ( i=0; i<nn; i++ )
        afft[i] = abs( comp[i] ) * dt ;
                              // 台形フィルターをかけ、積分
    StatusBar1->SimpleText = "フィルター操作" ;
    StatusBar1->Refresh() ;
    filter ( comp, dt, fll, flu, ful, fuu, nn ) ;
    StatusBar1->SimpleText = "積分：加速度→速度" ;
    StatusBar1->Refresh() ;
    infr ( comp, dt, nn, 0.0 ) ;
    for ( i=0; i<nn; i++ ) {
        vfft[i] = abs( comp[i] ) * dt * dt * nn / 2. / PI ;
        temp[i] = comp[i] ;
    }
                              // フーリエ逆変換で速度波形
    StatusBar1->SimpleText = "フーリエ逆変換：速度" ;
    StatusBar1->Refresh() ;
    fast ( temp, nn, +1 ) ;
    for ( i=0; i<nn; i++ ) {
        velr[i] = real( temp[i] ) / 2. / PI * dt ;
    }
                              // フィルターをかけ、さらに積分
    StatusBar1->SimpleText = "フィルター操作" ;
    StatusBar1->Refresh() ;
    filter ( comp, dt, fll, flu, ful, fuu, nn ) ;
    StatusBar1->SimpleText = "積分：速度→変位" ;
    StatusBar1->Refresh() ;
    infr ( comp, dt, nn, 0.0 ) ;
    for ( i=0; i<nn; i++ ) {
        dfft[i] = abs( comp[i] ) * dt * dt * dt * nn * nn / 4. / PI / PI ;
    }
                              // フーリエ逆変換で変位波形
    StatusBar1->SimpleText = "フーリエ逆変換：変位" ;
    StatusBar1->Refresh() ;
    fast ( comp, nn, +1 ) ;
    for ( i=0; i<nn; i++ ) {
        disr[i] = real( comp[i] ) * nn * dt * dt / 4. / PI / PI ;
    }
                              // 最大値を求める
    StatusBar1->SimpleText = "最大値の計算" ;
    StatusBar1->Refresh() ;
    for ( i=0; i<n; i++ )
        if ( fabs( accr[i] ) > acc_max ) acc_max = fabs(accr[i]);

    for ( i=0; i<nn; i++ ) {
        if ( fabs( velr[i] ) > vel_max ) vel_max = fabs(velr[i]);
        if ( fabs( disr[i] ) > dis_max ) dis_max = fabs(disr[i]);
    }

// ファイルに出力

    StatusBar1->SimpleText = "ファイルへ出力中" ;
    StatusBar1->Refresh() ;
    fprintf(fpout,"%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n",
        "時間(sec)","変位(cm)","速度(kine)","加速度(gal)",
        "周波数(Hz)","変位スペクトル(cm.sec)",
        "速度スペクトル(kine.sec)","加速度スペクトル(gal.sec)");
    for ( i=0; i<n; i++ )
        fprintf(fpout,"%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
            time[i], disr[i], velr[i], accr[i],
            freq[i], dfft[i], vfft[i], afft[i] );

    fclose(fpout);

    delete[] time ;
    delete[] freq ;
    delete[] accr ;
    delete[] velr ;
    delete[] disr ;
    delete[] afft ;
    delete[] vfft ;
    delete[] dfft ;

    sprintf( buffer , "%8.2f" , dis_max ) ;
    Form2->Label_dis1->Caption = buffer ;
    sprintf( buffer , "%8.2f" , vel_max ) ;
    Form2->Label_vel1->Caption = buffer ;
    sprintf( buffer , "%8.2f" , acc_max ) ;
    Form2->Label_acc1->Caption = buffer ;

    StatusBar1->SimpleText = "" ;
    StatusBar1->Refresh() ;
    Screen->Cursor = Save_Cursor;

    Form2->Show();

}
//---------------------------------------------------------------------------
// 入力加速度波形の選択

void __fastcall TForm1::bt_eqClick(TObject *Sender)
{
    if ( OpenDialog1->Execute() ){
        edit_eq->Text = OpenDialog1->FileName;
        edit_eq->Hint = edit_eq->Text;
    }

}
//---------------------------------------------------------------------------
// 結果保存ファイルの選択

void __fastcall TForm1::bt_outClick(TObject *Sender)
{
       if ( SaveDialog1->Execute() ){
        edit_output->Text = SaveDialog1->FileName;
        edit_output->Hint = edit_output->Text ;
    }
}
//---------------------------------------------------------------------------
// 説明文の表示

void __fastcall TForm1::BitBtn_helpClick(TObject *Sender)
{
    Form3->Show() ;

}
//---------------------------------------------------------------------------
// 高速フーリエ変換
//      参考文献：大崎順彦「新・地震動のスペクトル解析入門」鹿島出版会

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 = acos(-1.0) ;

    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::infr ( struct complex<double> c[], double dt, int nn, double v0 )
{
int nfold, k;
const double PI = acos(-1.0) ;
double c1, cr, pn;

    nfold = nn / 2 ;
    pn = PI / nn ;
    c1 = real ( c[0] ) ;
    cr = ( v0 / dt * nn + (nn-1) / 2. * c1 ) * pn;
	for ( k=1; k<nfold; k++ ) {
        cr = cr - imag( c[k] ) / (double)k;
        c[k] = complex<double>(-1.,cos( k*pn )/sin( k*pn )) * c1 * pn
               - complex<double>(0., 1.) * c[k] / (double)k;
        c[nn-k] = conj( c[k] ) ;
    }
    c[0] = complex<double>( cr*2. , 0. ) ;
    c[nfold] = complex<double>( -c1 * pn , 0. ) ;
}
//---------------------------------------------------------------------------
// フィルタ　（台形フィルタ） Band-pass, High-pass, Low-pass

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;
const double PI = acos(-1.0) ;
int kll, klu, kul, kuu, kll1, kuu1, i;

	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 ;
	if ( RadioButton1->Checked || RadioButton2->Checked ) {
		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 = 0.5-0.5*cos(((double)i/(double)nn/dt-fll)*PI/(flu-fll)) ;
			x[i] *= factor;
			x[nn-i] *= factor;
		}
	}
	if ( RadioButton1->Checked || RadioButton3->Checked ) {
		for ( i=kul+1; i<kuu; i++ ) {
			factor = 0.5+0.5*cos(((double)i/(double)nn/dt-ful)*PI/(fuu-ful)) ;
			x[i] *= factor;
			x[nn-i] *= factor;
		}
		for ( i=kuu; i<=kuu1; i++ ) x[i] = zero ;
	}
}
//---------------------------------------------------------------------------
// 基線補正　線形回帰

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] ;
}
//---------------------------------------------------------------------------
// パラメータの入力忘れチェック

bool TForm1::check_1()
{
bool nogood = false ;

    if ( edit_fll->Text == "" ) nogood = true ;
    if ( edit_flu->Text == "" ) nogood = true ;
    if ( edit_ful->Text == "" ) nogood = true ;
    if ( edit_fuu->Text == "" ) nogood = true ;
    if ( edit_number->Text == "" )  nogood = true ;
    if ( edit_dt->Text == "" )  nogood = true ;
    if ( edit_eq->Text == "" )  nogood = true ;
    if ( edit_output->Text == "" ) nogood = true ;
    if ( nogood ) {
		Application->MessageBox(L"パラメーターが入力されていない箇所があります。",
		L"エラー",MB_ICONEXCLAMATION | MB_OK ) ;
    }
    return ( nogood ) ;
}
//---------------------------------------------------------------------------
// 初期画面の大きさ設定

void __fastcall TForm1::FormCreate(TObject *Sender)
{
double ratio ;

    ratio = PixelsPerInch / 96. ;
    Form1->ClientHeight = 433 * ratio ;
    Form1->ClientWidth  = 480 * ratio ;
    ImageHigh->Visible = false ;
    ImageLow->Visible  = false ;

}
//---------------------------------------------------------------------------
// Band Pass Filter
void __fastcall TForm1::RadioButton1Click(TObject *Sender)
{
    LabelFll->Visible = true ;
    LabelFlu->Visible = true ;
    LabelFul->Visible = true ;
    LabelFuu->Visible = true ;
    Label1->Visible = true ;
    Label2->Visible = true ;
    Label3->Visible = true ;
    Label4->Visible = true ;
    edit_fll->Visible = true ;
    edit_flu->Visible = true ;
    edit_ful->Visible = true ;
    edit_fuu->Visible = true ;
    ImageBand->Visible = true ;
    ImageHigh->Visible = false ;
    ImageLow->Visible  = false ;
}
//---------------------------------------------------------------------------
// High Pass Filter
void __fastcall TForm1::RadioButton2Click(TObject *Sender)
{
    LabelFll->Visible = true ;
    LabelFlu->Visible = true ;
    LabelFul->Visible = false ;
    LabelFuu->Visible = false ;
    Label1->Visible = true ;
    Label2->Visible = false ;
    Label3->Visible = true ;
    Label4->Visible = false ;
    edit_fll->Visible = true ;
    edit_flu->Visible = true ;
    edit_ful->Visible = false ;
    edit_fuu->Visible = false ;
    ImageBand->Visible = false ;
    ImageHigh->Visible = true ;
    ImageLow->Visible  = false ;
}
//---------------------------------------------------------------------------
// Low Pass Filter
void __fastcall TForm1::RadioButton3Click(TObject *Sender)
{
    LabelFll->Visible = false ;
    LabelFlu->Visible = false ;
    LabelFul->Visible = true ;
    LabelFuu->Visible = true ;
    Label1->Visible = false ;
    Label2->Visible = true ;
    Label3->Visible = false ;
    Label4->Visible = true ;
    edit_fll->Visible = false ;
    edit_flu->Visible = false ;
    edit_ful->Visible = true ;
    edit_fuu->Visible = true ;
    ImageBand->Visible = false ;
    ImageHigh->Visible = false ;
    ImageLow->Visible  = true ;
}
//---------------------------------------------------------------------------