//---------------------------------------------------------------------------
// Ver.0.1    '04.7.29
// Ver.0.2    '04.8.10
// Ver.0.3    '04.8.26
// Ver.0.4    '08.6.9
// Ver.1.0    '10.7.18

#include <vcl.h>
#pragma hdrstop
#include <cmath>
#include <complex>
#include <except.h>

#include "main.h"
#include "about.h"
using namespace std;
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
	: TForm(Owner)
{
}
//---------------------------------------------------------------------------
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_execClick(TObject *Sender)
{
float x,y,z;
double s, dt, a, temp, *xin, *yin, *zin, *w ;
complex<double> *comp;
int n, total, i, nn, m, j;
char buffer[8];
string kai;

	shindo->Visible=False;
	shindokai->Visible=False;
	Form2->Hide() ;
	StatusBar1->SimpleText = "計算開始" ;
	StatusBar1->Refresh() ;

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

// ファイルのオープン

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

	dt = StrToFloat( edit_dt->Text );         // 時間刻み
	total = StrToInt( edit_number->Text );    // 地震波個数

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

	i = 0;
	while( fscanf(fpin,"%f %f %f",&x,&y,&z ) != EOF ){
		xin[i] = (double) x;
		yin[i] = (double) y;
		zin[i] = (double) z;
		i++ ;
		if ( i >= total ) break;
	}
	fclose(fpin);
	n = i;
	a = log( (double)(n-1) ) / log( 2. ) + 1.0 ;
	nn = pow ( 2., (int) a ) ;

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

	fourier ( xin, comp, dt, n, nn );
	fourier ( yin, comp, dt, n, nn );
	fourier ( zin, comp, dt, n, nn );

	for ( i=0; i<n; i++ ) {
		w[i] = sqrt( xin[i]*xin[i]+yin[i]*yin[i]+zin[i]*zin[i] );
	}

	for ( m=n-1; m>0; m-- )
		for ( i=0; i<m; i++ ) {
			if ( w[i] < w[i+1] ) {
				temp = w[i];
				w[i]=w[i+1];
				w[i+1]=temp;
			}
		}

	j=floor(0.3000001/dt)-1;           // 基準以上が0.3秒続く大きさ
	s=2.*log10( w[j] )+0.94;       // 計測震度の式
	s = floor( 100.*s + 0.5 );     // 小数第3位を四捨五入
	s = floor( s/10 ) / 10.;       // 小数第2位を切り捨て

	if (s<0.5) {
		kai="  0";
	} else if (s<1.5) {
		kai="  1";
	} else if (s<2.5) {
		kai="  2";
	} else if (s<3.5) {
		kai="  3";
	} else if (s<4.5) {
		kai="  4";
	} else if (s<5.0) {
		kai="  5弱";
	} else if (s<5.5) {
		kai="  5強";
	} else if (s<6.0) {
		kai="  6弱";
	} else if (s<6.5) {
		kai="  6強";
	} else {
		kai="  7";
	}

	delete[] comp ;
	delete[] xin ;
	delete[] yin ;
	delete[] zin ;
	delete[] w;

	sprintf( buffer , "%5.1f" , s ) ;
	shindo->Text = buffer ;
	shindokai->Text = kai.c_str() ;
	shindo->Visible=True;
	shindokai->Visible=True;

	StatusBar1->SimpleText = "計算終了" ;
	StatusBar1->Refresh() ;
}
//---------------------------------------------------------------------------
void TForm1::fourier ( double x[], complex<double> comp[], double dt, int n, int nn )
{
const complex<double> zero(0.,0.);
int i, nyquist;
double f,y,f1,f2,f3;
	                              // FFTの準備
    for ( i=0; i<n; i++ )
        comp[i] = complex<double>( x[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++ ) comp[i] /= (double)nn;
                              // フィルターをかける
    StatusBar1->SimpleText = "フィルター操作" ;
    StatusBar1->Refresh() ;

    nyquist=nn/2;
    comp[0]=zero;

    for ( i=1; i<=nyquist; i++ ) {
        f=(double)i/(double)nn/dt;
        y=f/10.;
        f1=sqrt(1./f);
        f2=1./sqrt(1.+0.694*y*y+0.241*pow(y,4)+0.0557*pow(y,6)
                +0.009664*pow(y,8)+0.00134*pow(y,10)+0.000155*pow(y,12));
        f3=sqrt(1.-exp(-(pow(f/0.5,3))));
        comp[i]=f1*f2*f3*comp[i];
    }
    for ( i=nyquist+1; i<nn; i++ ) {
        f=(double)(nn-i)/(double)nn/dt;
        y=f/10.;
        f1=sqrt(1./f);
        f2=1./sqrt(1.+0.694*y*y+0.241*pow(y,4)+0.0557*pow(y,6)
                +0.009664*pow(y,8)+0.00134*pow(y,10)+0.000155*pow(y,12));
        f3=sqrt(1.-exp(-(pow(f/0.5,3))));
        comp[i]=f1*f2*f3*comp[i];
    }
                              // フーリエ逆変換で時刻歴波形に戻す
    StatusBar1->SimpleText = "フーリエ逆変換" ;
    StatusBar1->Refresh() ;
    fast ( comp, nn, +1 ) ;
    for ( i=0; i<n; i++ ) {
        x[i] = real( comp[i] ) ;
    }
}
//---------------------------------------------------------------------------
// 高速フーリエ変換
//      参考文献：大崎順彦「新・地震動のスペクトル解析入門」鹿島出版会

void TForm1::fast ( complex<double> x[], int nn, int ind )
{
complex<double> temp, theta;
int i, j, k, m, kmax, istep, npower;
static const double PI = 6*asin( 0.5 ) ;

    npower = log( (double) nn + 1 ) / log( 2. ) ;
    for ( i=1; i<nn-1; i++ ) {
        m = 1 ;
        j = 0 ;
        for ( k=0; k<npower; k++ ) {
            j += (( i & m ) >> k) << ( npower - k - 1 ) ;
            m *= 2 ;
        }
        if ( i < j ) {
            temp = x[j] ;
            x[j] = x[i] ;
            x[i] = temp ;
        }
    }
    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 ;
    }
}
//---------------------------------------------------------------------------
// パラメータの入力忘れチェック

bool TForm1::check_1()
{
bool nogood = false ;
	if ( edit_number->Text == "" )  nogood = true ;
	if ( edit_dt->Text == "" )  nogood = true ;
	if ( edit_eq->Text == "" )  nogood = true ;
	if ( nogood ) {
		Application->MessageBox(L"パラメーターが入力されていない箇所があります。",
		L"エラー",MB_ICONEXCLAMATION | MB_OK ) ;
	}
	return ( nogood ) ;
}
//---------------------------------------------------------------------------
void __fastcall TForm1::BitBtn_helpClick(TObject *Sender)
{
	Form2->Show();
}
//---------------------------------------------------------------------------