#!/usr/bin/perl
# 以下は変更する必要なし
#-----------------------------------------------------------------------------
#
# csv2srs (version  $Id: eodv,v 1.4 2007/05/22 02:22:33 Exp $) 
#
# SRSの算出
# 
# Usage:
# > csv2srs.pl hogehoge.csv 
#
# options:
#         特に無し
#
# example:
# csv2srs.pl hogehoge.csv                                #(1)
#
#-----------------------------------------------------------------------------

# Libraries
my @dt, @dg, @alpa, @beta, @A, @B, @C, @D, @yy;  
my $tt,$bs,$bb,$ymax;

$Q		= 10;
$PI		= 3.1415926;
$g 		= 9.80665;		# [m/s2]
$delt		= 1.00E-06;		#  でるｔt
$zi		= 1/(2*$Q);
$wn		= 0.0;
$fbgn		= 1;			# [Hz]
$fmax		= 10000;		# [Hz]
$df		= 0;			# Δ[Hz]
$i 		= 0;
$cnt 		= 0;
$k 		= 0;
$m		= 1;
$c		= 1;
$oct		= 100;
$aa		= 0;

#-----------------------------------------------------------------------------
# Code begin
#-----------------------------------------------------------------------------
($opt) = @ARGV;

if ($opt eq "-h") {
    &show_usage();
    exit;
}

# オプションがない場合
if ($#ARGV == -1){
	print "Error:ファイル名を指定してください\n";
}
# オプションが１つある場合
if ($#ARGV == 0){


	# 対象ファイルのオープン
	open(IN,"<$ARGV[0]") || die "Error: can't open $ARGV[0]\n";
		@line = <IN>;
	close(IN);
	
	# SRS算出
	$cnt = -1;
	foreach (@line){
		# ヘッダ情報を飛ばす
		if ($cnt == -1){$cnt++; next;}

		# 時間[sec]／加速度取得[m/s2]
		$dt[$cnt+1] = (split(/\,/,$_))[0];
		$dg[$cnt+1] = (split(/\,/,$_))[1];
		
		# インクリメント処理
		$cnt++;
	}

	# ヘッダ出力
	$disp = sprintf("f[Hz],G[g],Q=%f\n",$Q);
	print $disp;
	
	$aa = &log10(($fmax))/$oct;
	for($j=$fbgn;$j<=$oct;$j++)
	{
		# 計算
		$wn   = 2*$PI*(10**($j*$aa));
		$ymax = &cal_acc();
			
		# 結果出力
		$disp = sprintf("%f,%e\n",$wn/(2*$PI),$ymax);
		print $disp;
	}
}

sub cal_acc()
{
	# 初期化
	$C[0] = $D[0] = $ymax = 0;
	
	# ヘッダ-debug
	#$disp = sprintf("#,time,応答G,αn,βn,Cn,Dn,Cn-1,Dn-1,a'(tn),b'(tn),c'(tn),d'(tn),An,Bn,An-1,Bn-1,Mn(tn),Mn-1(tn)\n");
	#print $disp;

	for($i=1;$i<$cnt;$i++)
	{
		# 加速度計算
		#$A[$i] = $B[$i] = $C[$i] = 0;
	
		# α,β計算・例外
		$alpa[$i] = (($dg[$i+1] - $dg[$i])/($delt))*$g;
		$beta[$i] = (($dg[$i]*$dt[$i+1] - $dg[$i+1]*$dt[$i])/($delt))*$g;
			
		# 特解計算
		$C[$i] = -$alpa[$i]/($wn*$wn);
		$D[$i] = -$beta[$i]/($wn*$wn) + 2*$zi*$alpa[$i]/($wn*$wn*$wn);

		# 一般解計算
		$tt = $dt[$i];
		$bs =  &fd($tt)*(&fM($i-1,$tt)-&fM($i,$tt))-&fb($tt)*($C[$i-1] - $C[$i]);
		$bb = &fa($tt)*&fd($tt)-&fb($tt)*&fc($tt);
		$A[$i] = $bs/$bb/exp(-$zi*$wn*$tt) + $A[$i-1];

		$bs = -&fc($tt)*(&fM($i-1,$tt)-&fM($i,$tt))+&fa($tt)*($C[$i-1] - $C[$i]);
		$B[$i] = $bs/$bb/exp(-$zi*$wn*$tt) + $B[$i-1];

		# yy出力
		#$yy[$i] = &fM($tt) + $A[$i]*&fa($tt) + $B[$i]*&fb($tt);
		$yy[$i] = -(&fM($i,$tt)+($A[$i]*&fa($tt) + $B[$i]*&fb($tt))*exp(-$zi*$wn*$tt))*$wn*$wn/$g;
		
		# ymaxの算出
		if(abs($yy[$i]) > $ymax){$ymax = abs($yy[$i]);}
		
		# 出力-debug
		#$disp = sprintf("%d,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e,%e\n",$i,$dt[$i],$yy[$i],$alpa[$i],$beta[$i],$C[$i],$D[$i],$C[$i-1],$D[$i-1],&fa($tt),&fb($tt),&fc($tt),&fd($tt),$A[$i],$B[$i],$A[$i-1],$B[$i-1],&fM($i,$tt),&fM($i-1,$tt));
		#print $disp;
	}

	return $ymax;
}

# オプションが２つある場合
if ($#ARGV == 1){	
}
exit;

#-----------------------------------------------------------------------------
# Sub code begin
#-----------------------------------------------------------------------------
sub fM()
{
		return $C[$_[0]]*$_[1] + $D[$_[0]];
}
sub fa()
{
	if($_[0] == 0){
		return 1;
	}else{
		#return exp(-$zi*$wn*$_[0])*cos(sqrt(1-$zi*$zi)*$wn*$_[0]);
		return cos(sqrt(1-$zi*$zi)*$wn*$_[0]);		
	}
}
sub fb()
{
	if($_[0] == 0){
		return 0;
	}else{
		#return exp(-$zi*$wn*$_[0])*sin(sqrt(1-$zi*$zi)*$wn*$_[0]);
		return sin(sqrt(1-$zi*$zi)*$wn*$_[0]);		
	}
}
sub fc()
{
	#return $wn*exp(-$zi*$wn*$_[0])*(-$zi*cos(sqrt(1-$zi*$zi)*$wn*$_[0])-sqrt(1-$zi*$zi)*sin(sqrt(1-$zi*$zi)*$wn*$_[0]));
	return $wn*(-$zi*cos(sqrt(1-$zi*$zi)*$wn*$_[0])-sqrt(1-$zi*$zi)*sin(sqrt(1-$zi*$zi)*$wn*$_[0]));
}
sub fd()
{
	#return $wn*exp(-$zi*$wn*$_[0])*(-$zi*sin(sqrt(1-$zi*$zi)*$wn*$_[0])+sqrt(1-$zi*$zi)*cos(sqrt(1-$zi*$zi)*$wn*$_[0]));
	return $wn*(-$zi*sin(sqrt(1-$zi*$zi)*$wn*$_[0])+sqrt(1-$zi*$zi)*cos(sqrt(1-$zi*$zi)*$wn*$_[0]));	
}
sub log10 {
    my $n = shift;
    return log($n)/log(10);
}
sub show_usage()
{
  print "EoDV(放電末期電圧)を抽出・可視化する\n";
  print "\n";
  print "Usage:\n";
  print "> eodv [-h]\n";
  print "\n";
  print "options:\n";
  print "         -h : 利用方法を提示する\n";
  print "      -make : Level-2を更新する\n";
  print "\n";
  print "example:\n";
  print "eodv -h                         #(1) HELP\n";
}


