#!/usr/bin/perl
use CGI;
use CGI::Carp (fatalsToBrowser);
use LWP::Simple;
print "Content-type: text/html;CHARSET=big5\n\n";

#########################################################
#  香港風力指數 version 3.0                             #
#  本程式碼提供利用Perl計算風力指數的算法，             #
#  可供大家用來參考。                                   #
#  如對以下程式有任何問題及意見，歡迎在香               #
#  港熱帶氣旋追擊站內之留言版提出。                     #
#                                                       #
#  請勿在未經同意的情況下修改風力指數計算的方法。       #
#                                                       #
#  程式及風力指數版權屬香港熱帶氣旋追擊站所有(C)        #
#  製作及編寫: JohnYeung                                #
#  2003年8月                                            #
#########################################################

@city = (25,30,20,22,30);  # 市區風速讀數陣列 (範例)
@urb  = (59,30,38,43,47);  # 離岸及高地風速讀數陣列 (範例)

# 範列執行結果: "本港最新風力指數: 40 (程度屬強勁)"

#####################   主計算程序   #####################

($u,$udv,$usize) = sdev(@city);   #市區平均
$ug = (2.28*$u)+10;               #市區陣風
$ug = int($ug*10**2+0.5)/10**2;
if ($u <= 13) {$ug = 0;}   # 區平均風速<=13時，不考慮陣風存在

($v,$vdv,$vsize) = sdev(@urb);    #空地平均
$vg = (1.58*$v) + 10;             #空地陣風
$vg = int($vg*10**2+0.5)/10**2;
if ($v <= 13) {$vg = 0;}

## 第一次風力指數計算 ##
$ur = gustfactor($u); $vr = gustfactor($v);  # 評估陣風比重
$wi = (0.56*((1-$ur)*$u+($ur*$ug))) + (0.44*((1-$vr)*$v+($vr*$vg)));

# 如果風力指數超過24km/h(強勁水平)，則移除異常值
if ($wi >= 24){     
	$cityref = \@city; $urbref = \@urb; # 製造參考 (用以把陣列傳入子程序)
	($u,$udv,$unewsize,@reudat) = remove($cityref,$u,$udv);  # 移除異常值
	($v,$vdv,$vnewsize,@revdat) = remove($urbref,$v,$vdv);

## 第二次風力指數計算 ##
$ur = gustfactor($u); $vr = gustfactor($v);
$wi = (0.56*((1-$ur)*$u+($ur*$ug))) + (0.44*((1-$vr)*$v+($vr*$vg)));

}

## 準備輸出 ##
if (($wi>= 0)and($wi< 14)) {$wd="輕微";$co="#008000";}
if (($wi>= 14)and($wi< 25)) {$wd="清勁";$co="#0000FF";}
if (($wi>= 25)and($wi< 45)) {$wd="強勁";$co="#FF0000";}
if (($wi>= 45)and($wi< 75)) {$wd="強烈";$co="#FF00FF";}
if ($wi>= 75) {$wd="極強";$co="#800040";}

$wi= int($wi*10+0.5)/10;  # RoundUp 至小數後兩個位

$report = "<strong><big>本港最新風力指數: <font color=$co> $wi (程度屬$wd)</font></big></strong>";

## 顯示結果 ##

print "$report";

#####################   計算子程序   #####################

# 計算陣風比重 
# (輸入:區平均風速 avg)
# (輸出:陣風比重 gustr)
sub gustfactor {
	my $avg = $_[0];
	my $gustr;
	
	$gustr = ((-0.00005)*($avg**2)) + (0.0081 * $avg ) - 0.0005;
	if ($avg > 81){$gustr = 0.327;}	
	return $gustr;	
}

# 計算標準差及平均 
# (輸入: 區資料陣列 data[])
# (輸出: 區平均風速 avg、區風速標準差 dev、區有效數據 size)
sub sdev { 
    my @data; 
    my ($size,$avg,$dev,$total);
    @data = @_;
    
    foreach $data(@data){
	    $total += $data;
	    if ($data > 0 ){
		    $size ++
	    }
    }
    $avg = $total/$size;
    
    foreach $data(@data){
	    if ($data > 0 ){
		    $dev += (($data - $avg)**2);
	    }	
    }
    $dev /= ($size-1);
    $dev = ($dev ** 0.5); 
    return ($avg , $dev , $size);
}

# 移除極端值  
# (輸入: 資料陣列參考 dataref、區風速標準差 dev、區平均風速 avg)
# (輸出: 新的平均 newavg、標準差 newdev、有效數據數目 newsize、被移除數據標示redat)

sub remove {
	$dataref = $_[0]; $dev = $_[2]; $avg = $_[1];
	my @data = @{$dataref};  # 解除參考
	my @redat;		
	
	for ($i=0;$i<=$#data;$i++){
		$check = $data[$i] - $avg;
		if ($check < 0){$check *= -1;}
		if ($check > $dev){
			$data[$i] = 0;
			$redat[$i] = "*"; # 被移除的數據以*號表示
		}					
	}
	($newavg,$newdev,$newsize) = sdev(@data);
	return ($newavg,$newdev,$newsize,@redat);
}