#!/usr/local/bin/perl

# Check the old forecasts to see how accurate they were.
# Calling sequence:
#
# CheckForecasts Directory
#
# OR
#
# CheckForecasts -delete Directory  to delete the forecast files after use

#   Copyright (C) 1996 Thomas R. Metcalf
#
#   This software is provided "as is" and is subject to change without
#   notice.  No warranty of any kind is made with regard to this software,
#   including, but not limited to, the implied warranties of
#   merchantability and fitness for a particular purpose.  The author shall
#   not be liable for any errors or for direct, indirect, special,
#   incidental or consequential damages in connection with the furnishing,
#   performance, or use of this software:  use it at your own risk.
#
#    This program is free software; you can redistribute it and/or
#    modify it under the terms of the GNU General Public
#    License as published by the Free Software Foundation; either
#    version 2 of the License, or (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
#    General Public License for more details.
#
#    You should have received a copy of the GNU General Public
#    License along with this library; if not, write to the Free
#    Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#


$# = "%.6g";

$DeleteFiles = 0;
for ($i=0; $i<=$#ARGV; ++$i) {
   if ($ARGV[$i] =~ /^-delete$/) {
      splice(@ARGV,$i,1);
      $DeleteFiles = 1;
   }
}

if ($DeleteFiles) {$DaysUnmodified=3.5;} else {$DaysUnmodified=0;}

$ForecastDirectory = $ARGV[0];

$TDATE = 21;
$ACTDATE = 23;
$OBSTYPE = 20;
$TLAT = 5;
$TLON = 7;
$NS = 6;
$EW = 8;
$TNAME = 9;
$TYEAR = 0;
$TMONTH = 1;
$TDAY = 2;
$TTIME = 3;
$TWIND = 16;
$TGUST = 17;

# Some constants;
 
$oneminute = 0.0000018974; # Times in years
$tenminutes = 0.000018974;
$onehour = 0.000113843;
$sixhours = 0.000683061; 
$oneday = 0.00273225;
$oneweek = 0.0191257;
 
$piover2 = atan2(1,0);
$pi = $piover2*2;
$dtor = $piover2/90;
$rtod = 90/$piover2;

# Get old forecast files

opendir(ForecastDirectory,"$ForecastDirectory") || die "can't open $ForecastDirectory";
@ForecastFiles = readdir(ForecastDirectory);
closedir(ForecastDirectory);
@ForecastFiles = grep(/.*\.forecast$/,@ForecastFiles);
@ForecastFiles = grep($_ = $ForecastDirectory."/".$_,@ForecastFiles);
@ForecastFiles = grep(-M $_ >= $DaysUnmodified,@ForecastFiles);
#print join("\n",@ForecastFiles);

for ($ffile=0; $ffile<=$#ForecastFiles; ++$ffile) {

   if (open(FFile,$ForecastFiles[$ffile])) {

      # Read old forecast data

      # Sort forecast data

      sub bydate {
          $na = @aa = split(/[ ]/,$a);
          $nb = @ab = split(/[ ]/,$b);
          $aa[$TDATE] <=> $ab[$TDATE];
      }

      @ForData = sort <FFile>;

      # Get corresponding actual data

      if ($ForecastFiles[$ffile] =~ /(.*)\.forecast$/) {
         if (open(ActFile,$1)) {

            @ActData = sort bydate <ActFile>;

            for ($fline=0; $fline<=$#ForData; ++$fline) {

               # Skip duplicate entries

               if ($ForData[$fline] ne $ForData[$fline+1]) {

                  @sfline = split(" ",$ForData[$fline]);
                  # Find true lat/lon and error for each forecast

                  ALOOP: for ($aline=1; $aline<=$#ActData; ++$aline) {
                     @saline = split(" ",$ActData[$aline]);
                     if  ($saline[$TDATE] >= $sfline[$TDATE] ) {
                        #interpolate to get lat/lon
                        @saline2 = split(" ",$ActData[$aline-1]);
                        $lat1 = $saline[$TLAT];
                        if ($saline[$NS] eq "S") {$lat1 = -$lat1;}
                        $lat2 = $saline2[$TLAT];
                        if ($saline2[$NS] eq "S") {$lat2 = -$lat2;}
                        $lon1 = $saline[$TLON];
                        if ($saline[$EW] eq "E") {$lon1 = -$lon1;}
                        $lon2 = $saline2[$TLON];
                        if ($saline2[$EW] eq "E") {$lon2 = -$lon2;}
                        $lat = &INTERPOLATE($saline[$TDATE],$saline2[$TDATE],
                                            $lat1,$lat2,$sfline[$TDATE]);
                        $lon = &INTERPOLATE($saline[$TDATE],$saline2[$TDATE],
                                            $lon1,$lon2,$sfline[$TDATE]);
                        $flat = $sfline[$TLAT];
                        if ($sfline[$NS] eq "S") {$flat = -$flat;}
                        $flon = $sfline[$TLON];
                        if ($sfline[$EW] eq "E") {$flon = -$flon;}
                        ($error,$errdirection) = &GCDISTANCE($lat,$lon,$flat,$flon);
                        last ALOOP;
                     }
                  }

                  # Output checked forecast.  
                  # Adds true lat, true lon, 
                  # forcast error (nmi), forecast error (direction)

                  $ForData[$fline] =~ s/\n//;
                  print ($ForData[$fline]," ",$lat," ",$lon," ",$error," ",$errdirection,"\n");

               }

            }

         }

      }

      if ($DeleteFiles) {unlink($ForecastFiles[$ffile]);}

   }

}

exit; 


#######################################################

sub INTERPOLATE {
 
   # Do a linear interpolation
 
   local($t1,$t2,$v1,$v2,$t) = @_;
 
   $v1+(($t-$t1)/($t2-$t1))*($v2-$v1);
 
}

 
sub ACOS {
 
# Compute the arc cosine
# Input: Cosine (forced to -1 <= cosine <= 1)
# Output: Angle in radians (0<angle<pi)
 
   local($cosine) = @_;
   local($argument);
   $argument = 1.0-$cosine**2;
   if ($argument<=0) {
      if ($cosine > 0) {0;} else {$pi;}
   } else {
      $piover2 - atan2($cosine/(sqrt($argument)),1);
   }
}

sub TAN {
   # tangent of $x (in radians)
   local($x) = @_;
   sin($x)/cos($x);
}

sub GCDISTANCE {
 
# Compute Great Circle distance between positions
# Input:  lat1, lon1, lat2, lon2, all in decimal degrees
# Output: (Distance in Nautical Miles, True course)
 
   local($latitude1,$longitude1,$latitude2,$longitude2) = @_;
   local($lat1) = $latitude1*$dtor;
   local($lat2) = $latitude2*$dtor;
   local($dlon) = ($longitude1-$longitude2)*$dtor;

   local($course) = atan2(sin($dlon),cos($lat1)*&TAN($lat2)-sin($lat1)*cos($dlon))*$rtod;
   if ($course < 0) {$course += 360;}
 
   (&ACOS(sin($lat1)*sin($lat2) + cos($lat1)*cos($lat2)*cos($dlon))*60/$dtor,
    $course
   );
}