#   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.
#

# PERL programs to compute the probability of a tropical cyclone hitting
# a latitude/longitude.
#
#   StrikeProb latitude longitude [filename] [location string]
#
#   latitude:   dd.mm  + = North
#                      - = South
#   Longitude:  dd.mm  + = West
#                      - = East
#
# 1996-03-19 TRM Added interpolation of forecast data for finer resolution
# 1996-03-22 TRM Added cumulative probabilities and changed the output to
#                use HTML tables.
# 1996-04-02 TRM Modified the calculation of the probabilities.  The 
#                calculation should be much more robust now.
# 1996-04-12 TRM A few more mods to the math.  Moved the city list
#                to a subroutine and alphabatized it.
# 1996-04-19 TRM Changed SLOPE from 9.5625 to 6.1
# 1996-04-23 TRM Fixed a bug in the stopping condition in &SPROB.  Also
#                use logs in SPROB calcs to make it more robust.
#                Moved constants out of &StrikeProb so that SPROB can be
#                called independently.
# 1996-05-07 TRM Fixed a bug in the interpolation subroutine.  Also added
#                check for distance between forecasts before interpolation.
# 1996-05-08 TRM Removed all code that assumed $STORMSIZE has 3 elements.
#                $STORMSIZE can now have arbitrary length.
# 1996-05-16 TRM Changed AVGSPROB to compute the average exactly.
#                Changes SPROB so deal with R=0 and S=0 properly.
#                Use incomplete gamma function in SPROB and AVGSPROB.
# 1996-05-20 TRM Changed the CUMMULATE subroutine.  More robust now.
# 1996-05-23 TRM Changed the CUMMULATE subroutine.  More robust now.
# 1996-05-28 TRM Changed the CUMMULATE subroutine.  Added The SPROBGIVEN
#                subroutine.  CUMMULATE/SPROBGIVEN combination is almost 
#                exact now.
#
# 1996-05-31 TRM VERSION 1.00
#
# 1997-01-08 TRM Fixed a bug in the YEAR2DATE subroutine which did leap
#                years not quite right.
# 1997-08-31 TRM Changed the way AVGSPROB is used by applying Bayes' theorem.
# 1997-09-19 TRM New AVGSPROB approximation + some sanity checks in AVGSPROB
#                Lot's of minor changes as well.
#
# 1997-09-19 TRM VERSION 1.02
#
# 2000-01-12 TRM Corrected another leap year bug in YEAR2DATE
# 2000-08-30 TRM Added smallscreen
#
# Version 1.4
#
# 2004-09-03 TRM Changed IDIST coeff to 1.0 from 3.0.  Changed interpolation
#                so that it will not interpolate if the strike prob is less
#                than 2%.
# 2004-09-07 TRM Added wml.
# 2004-09-12 TRM Check AVGSPROB when deciding whether or not to interpolate.
#
# Version 1.5
#
# 2004-09-16 TRM Better algorithm for AVGSPROB integral with variable step size.
#                Check probs from STORMSIZE[1] when deciding whether to interpolate.
# 2004-09-17 TRM Added $gofast approximation.
# 2004-09-23 TRM Added sanity check on cumprobs.
# 2004-09-27 TRM Added closest approach check to interpolation.
# 2004-09-30 TRM Added $INTERCEPT.  20 nmi.

package StrikeProb;
#require "/data1/WWW/Tropical/Bin/locations.pl";
#require "/data1/WWW/Tropical/Bin/RungeKutta.pl";

$VERSION = 1.50;

$SLOPE = 3.10;    # The time rate of change of the forecast error.  This
                  # is the central number which defines how the probabilites
                  # are computed.  Units:  nmi/hr.
                  # 1996-07-11 6.1 --> 4.5
                  # 1997-07-17 4.5 --> 4.0
                  # 2002-08-27 4.0 --> 3.2
                  # 2004-01-08 3.2 --> 3.0
                  # 2004-09-10 3.0 --> 3.1
$INTERCEPT = 20.0; # The error at time 0 in nmi, i.e. the error in the actual storm location
#$EPOWER = 1.18;   # Error = PSLOPE*(T**EPOWER)
$EPOWER = 1.00;  # 2004-09-21
$PSLOPE = (36**(1-$EPOWER))*$SLOPE;   # match slope calculation at 36 hours

# 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;
$sqrtpi = sqrt($pi);
$lnsqrtpi = log($sqrtpi);
$ln2 = log(2);

# These refer to the column in the tropical storm data base.

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

################### Main Subroutine ######################

sub StrikeProb {

# compute the probability of a tropical cyclone hitting
# a latitude/longitude.
#
#   StrikeProb latitude longitude [filename]
#
#   latitude:   dd.mm  + = North
#                      - = South
#   Longitude:  dd.mm  + = West
#                      - = East
#

   # Some Parameters

   @STORMSIZE = sort numerically (60,120,240);  # Radii of storm ... change at will

   $nbisect = 1;   # number of bisect passes for interpolation.  Compute time
                   # may increase rapidly with this number.
                   # 0 = No interpolation
                   # 1 = Interpolate at midpoints
                   # 2 = Interpolates at midpoints and midpoints of midpoints, etc.

   $IDIST = $STORMSIZE[0]*1.0;  # The maximum distance between forecasts (nmi) before interpolation.  If the distance
                              # is less than $IDIST, there is no interpolation; if greater, a new forecast
                              # is created by bisection, but the number of bisections is limited
                              # by $nbisect above.

   # End of parameters, start of code

   $- = 0;  # force a top-of-page
   $| = 1;  # force STDOUT to be unbuffered

   @argv = @_;

   $storm_near = 0;

   ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = gmtime;
   if ($year>=96) {$year+=1900;} else {$year+=2000;}
   $now = $year + ($yday+1 + ($hour/24.0) + ($min/1440.))/366.0;

   # Get the name of the file which holds the forecast data

   #print join(" ",@argv);
   #return;

   if ($#argv >= 2) {
      $tropfile = $argv[2];
   } else {
      $tropfile = '-';
   }

   if ($#argv >= 3) {
      $locstring = $argv[3];
   } else {
      $locstring = "";
   }

   if ($#argv >= 4) {
      $smallscreen = $argv[4];
   } else {
      $smallscreen = 0; 
   }

   if ($#argv >= 5) {
      $wml = $argv[5];
   } else {
      $wml = 0;
   }

   if ($#argv >= 6) {
      $gofast = $argv[6];
   } else {
      $gofast = 0;
   }

   if ($#argv >= 7) {
      $nowoverride = $argv[7];
   } else {
      $nowoverride = 0;
   }



   if ($#argv < 0) {print "<p>Latitude and Longitude must be specified</p>\n"; exit;}

   # Define lat and lon for various cities.

   $lat = $argv[0];
   if (!&ISNUMBER($lat)) {
      ($lat,$lon,$locstring) = &CITY($argv[0]);
      if (!&ISNUMBER($lat) || !&ISNUMBER($lon)) {
         print "<p>Unreckognized location code: $argv[0]</p>\n";
         exit;
      }
   } else {
      $lon = $argv[1];
      if (!&ISNUMBER($lon)) {print "<p>Longitude must be a number $lon</p>\n"; exit;}
   }
   $lat = int($lat) + ($lat-int($lat))*1.66666667;
   $lon = int($lon) + ($lon-int($lon))*1.66666667;

   # parse the forecast data by storm name

   open(TROPICAL,$tropfile) || die "Couldn't open $tropfile";

   ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = gmtime;
   if ($year>=96) {$year+=1900;} else {$year+=2000;}
   $now = $year + ($yday+1 + ($hour/24.0))/366.0;

   if ($nowoverride) {$now = $nowoverride;}

   $ptr = 0;
   $actdate = 99999;

   while(<TROPICAL>) {
      s/  / /g;  # Make sure there are no multiple spaces
      @line = split(/ /);
      $joinline = join(" ",@line);
      if ($#line >= 23) {    # This is a valid line
         if ($line[$DATE] > ($now-$oneday)) {
            if ($line[$OBSTYPE] eq 'FOR') {    # This is a current forecast
               if ($line[$ACTDATE] > ($now-1.0*$oneday)) {
                  if ($lastline{$line[$TNAME]} != "" && 1) {   # interpolate
                     @interpline = ();
                     @newinterpline = ();
                     $interpline[0] = $lastline{$line[$TNAME]};
                     $interpline[1] = $joinline;
                     for ($bisect=1; $bisect<=$nbisect+1; ++$bisect) {  # +1 to do pca check
                        my ($pcacheck);  # check for point of closest approach?
                        if ($bisect <= $nbisect) {$pcacheck=0;} else {$pcacheck=1;}
                        $ilineptr = 0;
                        for ($biline=0; $biline<$#interpline; ++$biline) {
                           my($lat1,$lon1,$lat2,$lon2,$tim1,$tim2,
                                 $good1,$good2,$out1,$out2,$K12,$strike1,$strike2,$nstrike);
                           $newinterpline[$ilineptr++] = $interpline[$biline];
                           $newinterpline[$ilineptr++] = &INTERPOLATE($interpline[$biline],$interpline[$biline+1]);
                           # Check the distance\time between forecasts.  
                           # If it is less than $IDIST nmi, or less than 12 hours,
                           # don't use the interpolation.  Also don't use if the strike prob is <5%
                           ($lat1,$lon1,$tim1) = &GETLATLON($interpline[$biline]);
                           ($lat2,$lon2,$tim2) = &GETLATLON($interpline[$biline+1]);
                           my $testdist = &GCDISTANCE($lat1,$lon1,$lat2,$lon2);
                           # Compute relevant interpolated forecast data.
                           ($good,$out)=&DOSPROB($lat,$lon,split(/ /,$newinterpline[$ilineptr-1]));      # sets actdate too
                           if ($good) {
                              my($testname,$testyear,$testmonth,$testday,$testhour,$testminute,$testydate,$testforwind,$testforlat,$testforlon,$testfortdiff,@teststrike)=split(/ /,$out);
                              my($lasttestname,$lasttestyear,$lasttestmonth,$lasttestday,$lasttesthour,$lasttestminute,$lasttestydate,$lasttestforwind,$lasttestforlat,$lasttestforlon,$lasttestfortdiff,@lastteststrike)=split(/ /,$outline[$ptr-1]);
                              # interpolate if we have found a minimum in the distance
                              my $testdist = &GCDISTANCE($lat,$lon,$testforlat,$testforlon);
                              my $testdist1 = &GCDISTANCE($lat,$lon,$lat1,$lon1);
                              my $testdist2 = &GCDISTANCE($lat,$lon,$lat2,$lon2);
                              #print "<!-- INTERP: $testday $testhour: $testdist $testdist1 $testdist2 -->\n";
                              my $testoverlap = &AVGSPROB($STORMSIZE[0],0,$testdist,($tim2-$tim1)/$onehour,0.5,0.5);
                              my $testplace = 1;
                              if ($testplace > $#STORMSIZE || $gofast) {$testplace = $#STORMSIZE;}
                              if ( ( !$pcacheck &&
                                     ( ( ($teststrike[$testplace] >= 0.02 || $lastteststrike[$testplace] >=0.02) && 
                                         (($testdist>=$IDIST) || ($testoverlap<0.6)) 
                                       ) || ((($tim2-$tim1)/$onehour)>12.1)
                                     ) 
                                   ) || ( # check for point of closest approach
                                          ($testdist <= $testdist1 && $testdist <= $testdist2)
                                        )
                                 ) {
                                 # Save relevant interpolated forecast data for output.
                                 $outline[$ptr] = $out;
                                 $outline[$ptr] =~ s/  / /g;
                                 ++$ptr;
                              }
                           }
                           $nxtbiline = $interpline[$biline+1];
                        }
                        $newinterpline[$ilineptr++] = $nxtbiline;
                        @interpline = @newinterpline;
                     }
                  }
                  # Compute and save relevant forecast data for output
                  ($good,$out)=&DOSPROB($lat,$lon,@line);      # sets actdate too
                  if ($good) {
                     $outline[$ptr] = $out;
                     $outline[$ptr] =~ s/  / /g;
                     ++$ptr;
                  }
               }
            } else {  # Actual data: Is there a storm nearby now?
               # This also saves the last actual data.  I assume that the storm
               # data base is sorted!  (It should be ...)
               $sname = substr($line[$TNAME],0,length($line[$TNAME])-3); 
               $actlat{$sname} = $line[$LATITUDE];
               $actlon{$sname} = $line[$LONGITUDE];
               if ($line[$NS] eq 'S') {$actlat{$sname} = -$actlat{$sname};}
               if ($line[$EW] eq 'E') {$actlon{$sname} = -$actlon{$sname};}
               $actdistance{$sname} = &GCDISTANCE($lat,$lon,$actlat{$sname},$actlon{$sname});
               $actcourse{$sname} = &GCCOURSE($lat,$lon,$actlat{$sname},$actlon{$sname});
               $actdate{$sname} = $line[$DATE];
               my $tyear = $line[$TYEAR];
               my $tmonth = $line[$TMONTH];
               my $tday = $line[$TDAY];
               my $thour = int($line[$TTIME]);
               my $tminute = int(0.5+100*($line[$TTIME]-int($line[$TTIME])));
               if (length($tmonth)<2) {$tmonth = "0".$tmonth;}
               if (length($tday)<2) {$tday = "0".$tday;}
               if (length($thour)<2) {$thour = "0".$thour;}
               if (length($tminute)<2) {$tminute = "0".$tminute;}
               $actdatetext{$sname} = "$tyear-$tmonth-$tday $thour:$tminute UT";
               $acttdiff = ($now - $line[$DATE])/$oneday;  # in days
               if ($acttdiff <= 1 && $actdistance{$sname} <= $STORMSIZE[$#STORMSIZE]) {
                  $storm_near = 1;
                  $storm_near{$sname} = 1;
               }
            }
         }
         $lastline{$line[$TNAME]} = $joinline;
      }
   }

   if ($lat < 0) {$ns='S'; $lat=-$lat;} else {$ns='N';}
   if ($lon < 0) {$ew='E'; $lon=-$lon;} else {$ew='W';}
   $latmin = int(($lat-int($lat))*60);
   $lonmin = int(($lon-int($lon))*60);
   if (length($latmin) < 2) {$latmin = "0".$latmin;}
   if (length($lonmin) < 2) {$lonmin = "0".$lonmin;}
   $latstring = int($lat).'&#176;'.$latmin.'&#180;';
   $lonstring = int($lon).'&#176;'.$lonmin.'&#180;';
   if ($locstring eq "") {
      $locstring = $latstring.$ns." ".$lonstring.$ew;
      $locstringtext = $locstring;
   } else {
      $locstringtext = $locstring;
      $locstring = $locstring.": ".$latstring.$ns." ".$lonstring.$ew;
   }

   # Now do the HTML table output and compute the cumulative probabilities

   $= = 99999;  # Only one page
   if ($ptr <= 0) {
      if ($storm_near) {
         foreach $key (keys(%storm_near)) {
            $name = $key;
            print "<p>No forecasts available for storm $name near $locstring</p>\n";

            my $locstringhere = int($actdistance{$name})." nmi ".&DEG2DIRECTION($actcourse{$name});
            my $timestring = $actdatetext{$name};
            my $ftimediff = int(($now - $actdate{$name})/$onehour+0.5);
            my $ftimediffplural;
            if ($ftimediff != 1) {$ftimediffplural = "s";}
            my $ftimediffstring = " (roughly $ftimediff hour$ftimediffplural ago)";
            if (!($timestring =~ /\d\d\d\d\-\d\d\-\d\d\s+\d\d\:\d\d\s+UT/i)) {
               $timestring = "last update";
            } else {
               $timestring = $timestring . $ftimediffstring;
            }
            print "<p>At $timestring, $name was $locstringhere of $locstringtext</p>\n";

            if ($smallscreen) {print "<hr><p><a href=\"file:Hurricane.pqa/hurricane.html\">Home</a><p>";}
            if ($wml) {
               print "<p><anchor accesskey=\"1\">Prev<prev/></anchor></p>\n";
               print "<p><anchor accesskey=\"2\">Home<go href=\"/Tropical/HTML/WML/index.wml#home\"/></anchor></p>\n";
               print "<p><anchor accesskey=\"3\">Strike Probabilities<go href=\"/Tropical/HTML/WML/strike.wml\"/></anchor></p>\n";
            }
         }
      } else {
         print "<p>No tropical cyclones near $locstring</p>\n";
         if ($smallscreen) {print "<hr><p><a href=\"file:Hurricane.pqa/hurricane.html\">Home</a><p>";}
         if ($wml) {
            print "<p><anchor accesskey=\"1\">Prev<prev/></anchor></p>\n";
            print "<p><anchor accesskey=\"2\">Home<go href=\"/Tropical/HTML/WML/index.wml#home\"/></anchor></p>\n";
            print "<p><anchor accesskey=\"3\">Strike Probabilities<go href=\"/Tropical/HTML/WML/strike.wml\"/></anchor></p>\n";
         }
      }
   } else {
      $ssizestring = "";
      $ssizerowhead1 = "";
      $ssizerowhead2 = "";
      for ($s=0; $s<=$#STORMSIZE; $s++) {
         $ssize = int($STORMSIZE[$s]);
         $ssizestring = $ssizestring."$ssize nmi";
         if ($s < $#STORMSIZE) {$ssizestring = $ssizestring.", ";}
         if ($s+1 == $#STORMSIZE) {$ssizestring = $ssizestring."and ";}
         if ($smallscreen) {
            $ssizerowhead1 = $ssizerowhead1."<td align=center>$ssize</td>";
         } else {
            $ssizerowhead1 = $ssizerowhead1."<TD ALIGN=CENTER COLSPAN=3 NOWRAP><B>$ssize nmi</B></TD><TD NOWRAP></TD>";
            $ssizerowhead2 = $ssizerowhead2."<TD ALIGN=CENTER NOWRAP>T</TD><TD ALIGN=CENTER NOWRAP>I</TD><TD ALIGN=CENTER NOWRAP>C</TD><TD NOWRAP></TD>";
         }
      }
      if ($smallscreen) {
         print "Strike Probabilities for $locstring\n";
         print "<table border=1 width=153>\n";
         print "<tr><td align=center rowspan=2>Name</td><td align=center rowspan=2>Time</td><td align=center rowspan=1 colspan=$s>Distance (nmi)<td align=center rowspan=2>Wind</td></tr>\n";
         print "<tr>$ssizerowhead1</tr>\n";
      } elsif ($wml) {
         my $numcols = $s+3;
         print "<p>Strike Probabilities for $locstring.</p>\n";
      } else {
         print "<P ALIGN=CENTER><H3>Tropical Cyclone Strike Probabilities for $locstring</H3></P>\n";
      print "<P>Probability that the storm center will be within $ssizestring of $locstring.  Probabilities at intermediate times may be higher than those shown.  The probabilities are given in percent (\"&lt;1\" means \"less than one percent probability\") and are updated whenever a new forecast is received.<P>\n";
         print "<P><DIV ALIGN=CENTER><TABLE BORDER=1>\n";
         print "<TR><TD ALIGN=CENTER ROWSPAN=2 NOWRAP><B>Name</B></TD><TD ALIGN=CENTER ROWSPAN=2 NOWRAP><B>Date</B></TD><TD ALIGN=CENTER ROWSPAN=2 NOWRAP><B>Time</B></TD><TD ALIGN=CENTER ROWSPAN=1 COLSPAN=2WRAP><B>Relative Time</B></TD><TD NOWRAP></TD>$ssizerowhead1<TD ALIGN=CENTER ROWSPAN=2 NOWRAP><B>Wind</B></TD></TR>\n";
         print "<TR><TD ALIGN=CENTER NOWRAP>Now</TD><TD ALIGN=CENTER NOWRAP>Fcst</TD><TD NOWRAP></TD>$ssizerowhead2</TR>\n";
      }
      $outline[$ptr]=" ";
      @outline = @outline[0..($ptr-1)];   # Get rid of possible garbage
      @outline = sort @outline;           # Sort output by name and date
      $lastname = "";
      for ($i=0; $i<$ptr; $i++) {
         ($name,$year,$month,$day,$hour,$minute,$ydate,$forwind,$forlat,$forlon,$fortdiff,@strike)=split(/ /,$outline[$i]);
         ($nextname,$nextyear,$nextmonth,$nextday,$nexthour,$nextminute,$nextydate,$nextforwind,$nextforlat,$nextforlon,$nextfortdiff,@nextstrike)=split(/ /,$outline[$i+1]);
         if ($#strike == $#STORMSIZE) {

            if ($lastname ne $name) {   # start of a new storm
               $lastreltime2 = 0;
               if ($smallscreen) {
                  print "<tr><td colspan=6></td></tr>\n";
               } elsif ($wml) {
               } else {
                  print "<TR><TD COLSPAN=17 NOWRAP></TD></TR>\n";
               }
               @llarchive = ();
               $lastforlat = $actlat{$name};
               $lastforlon = $actlon{$name};
               $lastdate = $actdate{$name};
               $lastfortdiff = 0;
               for ($s=0; $s<=$#STORMSIZE; $s++) {
                  if (defined($actdistance{$name}) && $actdistance{$name} <= $STORMSIZE[$s]) {
                     $laststrike[$s] = 1;
                  } else {
                     $laststrike[$s] = 0;
                  }
                  $lastcumprob[$s] = 0;
                  $lastintervalprob[$s] = 0;
               }
               $llarchive[0] = join(" ",($lastdate,$lastforlat,$lastforlon,$lastfortdiff,@laststrike)); 
            } 
            # Store info for calculation of cumulative probability
            $fortdiff = ($ydate-$actdate{$name})/$onehour;
            $llarchive[$#llarchive+1] = join(" ",($ydate,$forlat,$forlon,$fortdiff,@strike)); 
            $tdiff = ($ydate-$lastdate)/$onehour;
            for ($s=0; $s<=$#STORMSIZE; $s++) {  # get cumulative probs

               # Cumulative probability is
               # P(A+B|M) = P(A|M) + P(B|M) - P(AB|M)
               # where P(AB|M) = P(B|M)P(A|BM)
               # P(A+B|M) = P(A|M) + P(B|M) - P(B|M)P(A|BM) 
               # Here, A is the cumulative prob and B is the prob at the
               # next time.
               # $overlap and $cumoverlap are P(A|BM), the cumulative prob up 
               # to the current time given that the storm is within STORMSIZE nmi
               # of the next forecast.  This requires going through all the 
               # forecasts again which is why they were saved in @llarchive.

               # $overlap is the probability that the storm is within
               # $STORMSIZE[$s] nmi at the time of the current forecast
               # given that the storm was within $STORMSIZE[$s] nmi at
               # the time of the previous forecast.

               # integral average of K(r,s,e)

               #$overlap[$s] = &AVGSPROB($STORMSIZE[$s],$fortdiff,&GCDISTANCE($forlat,$forlon,$lastforlat,$lastforlon),$tdiff,$strike[$s],$laststrike[$s]);
               
               my(@given) = ($#llarchive-1);
               $overlap[$s] = &SPROBGIVEN($s,$#llarchive,\@given,\@llarchive);
               my ($upper,$lower);
               if ($laststrike[$s] > 0) {
                  $upper = $strike[$s]/$laststrike[$s];
                  $lower = ($strike[$s]+$laststrike[$s]-1)/$laststrike[$s];
               } else {
                  $upper = 1;
                  $lower = 0;
               }
               if ($overlap[$s]>$upper) {$overlap[$s]=$upper;}
               if ($overlap[$s]<$lower) {$overlap[$s]=$lower;}
               #if ($overlap[$s]<$lower) {$overlap[$s]=$strike[$s];}
               #if ($overlap[$s]>$upper || $overlap[$s]<$lower) {$overlap[$s]=$strike[$s];}
               $intervalprob[$s] = $strike[$s] 
                                 + $laststrike[$s] 
                                 - $overlap[$s]*($laststrike[$s]);
               #print "<!-- interval: $strike[$s] $laststrike[$s] $overlap[$s] $intervalprob[$s] -->\n";

               if ($intervalprob[$s] > 1) {$intervalprob[$s] = 1;}
               if ($intervalprob[$s] < 0) {$intervalprob[$s] = 0;}
               #if ($intervalprob[$s] < $laststrike[$s]) {$intervalprob[$s] =$laststrike[$s];}
               #if ($intervalprob[$s] < $strike[$s]) {$intervalprob[$s] =$strike[$s];}

               if ($lastname ne $name) {   # Start of new storm?
                  $cumprob[$s] = $laststrike[$s];
               }
               my $gotviolation = 0;
               if ($cumprob[$s] < 1) {  # if cumprob is 1, we can't add any more
                  my(@given) = ($#llarchive);
                  if ($gofast) {  # Fast, but approximate
                     $cumoverlap[$s] = &SPROBGIVEN($s,$#llarchive-1,\@given,\@llarchive);
                  } else {
                     $cumoverlap[$s] = &CUMMULATE($#llarchive-1,$s,\@given,\@llarchive);
                  }
                  my ($upper,$lower);
                  if ($strike[$s] > 0) {
                     $upper = $cumprob[$s]/$strike[$s];
                     $lower = ($strike[$s]+$cumprob[$s]-1)/$strike[$s];
                  } else {
                     $upper = 1;
                     $lower = 0;
                  }
                  if ($cumoverlap[$s]>$upper) {$cumoverlap[$s]=$upper;}
                  if ($cumoverlap[$s]<$lower) {$cumoverlap[$s]=$lower; $gotviolation=1;}
                  #if ($cumoverlap[$s]<$lower) {$cumoverlap[$s]=$cumprob[$s]; $gotviolation=1;}
                  $cumprob[$s] += $strike[$s] - $cumoverlap[$s]*$strike[$s];
               }

               # The next block does a crude correction for any inaccuracies
               # in CUMMULATE.  It does not make a very large correction.
               # The first condition is (S_{i}-S{i-1})>(I_{i}-I_{i-1}). 
               # Since S_{i} = S_{i-1}+A_{i}-A_{i}(S_{i-1}|A_{i}) and
               # I_{i} = A_{i-1}+A_{i}-A_{i}(A_{i-1}|A_{i}),
               # (S_{i}-S{i-1})>(I_{i}-I_{i-1}) is the same as
               # A_{i} [(S_{i-1}|A_{i}) - (A_{i-1}|A_{i})] < I_{i-1} - A_{i-1}.
               # Since I_{i-1} - A_{i-1} > 0, the condition is the same as
               # (S_{i-1}|A_{i}) < (A_{i-1}|A_{i}).   In the gofast
               # approximation, (S_{i-1}|A_{i}) = (A_{i-1}|A_{i}) so this
               # condition should always be satisfied under the gofast approximation.
               # It is not satisfied otherwise, but the approximation is good, so use
               # the check if there was a violation above.  Otherwise use the condition
               # that S{i} > I{i} which is always true.
               # S represents the cummulative probabilities and I represents the interval probs.
               my $testcumprob = $intervalprob[$s]+$lastcumprob[$s]-$lastintervalprob[$s];
               if ($cumprob[$s] < $testcumprob && $testcumprob < 1 && $testcumprob >0 && ($gofast || $gotviolation)) {
                  $cumprob[$s] = $testcumprob;
               } else {
                  if ($cumprob[$s] < $intervalprob[$s]) {$cumprob[$s] = $intervalprob[$s]}
               }
               if ($cumprob[$s] > 1) {$cumprob[$s] = 1;}
               if ($cumprob[$s] < 0) {$cumprob[$s] = 0;}

               $lastcumprob[$s] = $cumprob[$s];    # save previous values
               $lastintervalprob[$s] = $intervalprob[$s];
            }

            # Prepare and print all the output strings.

            $storm_near{$name} = 0;  # There was forecast data for this storm
            if (length($month) < 2) {$month = "0".$month;}
            if (length($day) < 2) {$day = "0".$day;}
            if (length($hour) < 2) {$hour = "0".$hour;}
            if (length($minute) < 2) {$minute = "0".$minute;}
            if (length($minute) < 2) {$minute = "0".$minute;}
            while (length($forwind) < 3) {$forwind = "0".$forwind;}
            $date = $year."-".$month."-".$day;
            $time = $hour.":".$minute." UT";
            $reltime1 = int(($ydate - $now)/$onehour+0.5);
            $reltime2 = int(($ydate - $actdate{$name})/$onehour+0.5);
            if ($smallscreen != 2) {
               if ($reltime1 > 0) {$reltime1 = "+"."$reltime1";}
               else {$reltime1=-$reltime1; $reltime1 = "-"."$reltime1";}
               if ($reltime2 > 0) {$reltime2 = "+"."$reltime2";}
               else {$reltime2=-$reltime2; $reltime2 = "-"."$reltime2";}
            }
            for ($si=0; $si<=$#STORMSIZE; $si++) {
               $sstrike[$si] = int($strike[$si]*100); if ($sstrike[$si]<1) {$sstrike[$si]="&lt;1";}
               $istrike[$si] = int($intervalprob[$si]*100); if ($istrike[$si]<1) {$istrike[$si]="&lt;1";}
               $cstrike[$si] = int($cumprob[$si]*100); if ($cstrike[$si]<1) {$cstrike[$si]="&lt;1";}
               if ($sstrike[$si]>99) {$sstrike[$si]=99;} 
               if ($istrike[$si]>99) {$istrike[$si]=99;} 
               if ($cstrike[$si]>99) {$cstrike[$si]=99;} 
            } 
            if ($smallscreen) {
               my $shortname = substr($name,0,5);
               print "<tr>\n";
               print "   <td align=center>$shortname</td>\n";
               #print "   <td align=center>$reltime2 h</td>\n";
               if ($smallscreen == 2) {
                  print "   <td align=center>$lastreltime2-$reltime2 h</td>\n";
               } else {
                  print "   <td align=center>$reltime2 h</td>\n";
               }
               for ($si=0; $si<=$#STORMSIZE; $si++) {
                  my($sstrike,$istrike,$cstrike) = ($sstrike[$si],$istrike[$si],$cstrike[$si]);
                  if ($smallscreen == 2) {
                     print "   <td align=right>$istrike</td>\n";
                  } else {
                     print "   <td align=right>$cstrike</td>\n";
                  }
               }
               print "   <td align=center>$forwind</td>\n";
               print "</tr>\n";
            } elsif ($wml) {
               if ($name ne $nextname) {
                  #my $shortname = substr($name,0,5);
                  print "<p>$name strike probabilities over the next $reltime1 hours:</p>\n";
                  for ($si=0; $si<=$#STORMSIZE; $si++) {
                     my($sstrike,$istrike,$cstrike) = ($sstrike[$si],$istrike[$si],$cstrike[$si]);
                     print "<p>$STORMSIZE[$si] nmi: $cstrike\%</p>\n";
                  }
               }
            } else {
               print "<TR>\n";
               print "   <TD ALIGN=CENTER NOWRAP>$name</TD>\n";
               print "   <TD ALIGN=CENTER NOWRAP>$date</TD>\n";
               print "   <TD ALIGN=CENTER NOWRAP>$time</TD>\n";
               print "   <TD ALIGN=CENTER NOWRAP>$reltime1 h</TD>\n";
               print "   <TD ALIGN=CENTER NOWRAP>$reltime2 h</TD>\n";
               print "   <TD ALIGN=CENTER NOWRAP></TD>\n";
               for ($si=0; $si<=$#STORMSIZE; $si++) {
                  my($sstrike,$istrike,$cstrike) = ($sstrike[$si],$istrike[$si],$cstrike[$si]);
                  print "   <TD ALIGN=RIGHT NOWRAP>$sstrike</TD>\n";
                  print "   <TD ALIGN=RIGHT NOWRAP>$istrike</TD>\n";
                  if ((split(" ",$outline[$i+1]))[0] ne $name) {
                     print "   <TD ALIGN=RIGHT NOWRAP><B>$cstrike</B></TD>\n";
                  } else {
                     print "   <TD ALIGN=RIGHT NOWRAP>$cstrike</TD>\n";
                  }
                  print "   <TD ALIGN=CENTER NOWRAP></TD>\n";
               }
               print "   <TD ALIGN=CENTER NOWRAP>$forwind</TD>\n";
               print "</TR>\n";
            }
            $lastreltime2 = $reltime2;

            # Save the last forecast for use during the next forecast.

            ($lastname,$lastyear,$lastmonth,$lastday,$lasthour,$lastminute,$lastdate,$lastforwind,$lastforlat,$lastforlon,@laststrike)=($name,$year,$month,$day,$hour,$minute,$ydate,$forwind,$forlat,$forlon,@strike)
         }
      }
      if ($smallscreen ) {
         print "</table>\n";
      } elsif ($wml) {
      } else {
         print "</TABLE></DIV>\n";
      }
      foreach $key (keys(%storm_near)) {
         if ($storm_near{$key}) {  # There was no forecast data for this storm
            $name = $key;
            print "<p>No forecasts available for storm $name near $locstring.</p><p></p>\n";
         }
      }
      print "<p></p>\n";
      foreach $name (keys(%actdate)) {
         $isnear = grep(/^$name$/i,keys(%storm_near)); # Is this $name nearby?
         if ($actdate{$name} ne "" && $isnear>0) {
            ($foryear,$formon,$formday,$forhour,$formin) = &YEAR2DATE($actdate{$name});
            $fordate = $foryear."-".$formon."-".$formday." ".$forhour.":".$formin." UT";
            $year2 = substr($foryear,2,2);
            print "<p>\n";
            if ($smallscreen || $wml)  { 
               print "$name probabilities based on forecast data from $fordate. \n";
            } else {
               print "<A HREF=\"\/Tropical\/Data\/$name-$year2.html\">$name</A> probabilities based on forecast data from $fordate.<br>\n";
            }
            if (defined($actdistance{$name})) {
               my $locstringhere = int($actdistance{$name})." nmi ".&DEG2DIRECTION($actcourse{$name});
               my $timestring = $actdatetext{$name};
               my $ftimediff = int(($now - $actdate{$name})/$onehour+0.5);
               my $ftimediffplural;
               if ($ftimediff != 1) {$ftimediffplural = "s";}
               my $ftimediffstring = " (roughly $ftimediff hour$ftimediffplural ago)";
               if (!($timestring =~ /\d\d\d\d\-\d\d\-\d\d\s+\d\d\:\d\d\s+UT/i)) {
                  $timestring = "last update";
               } else {
                  $timestring = $timestring . $ftimediffstring;
               }
               print "At $timestring, $name was $locstringhere of $locstringtext\n"
            }
            print "</p><p></p>\n";
         }
      }
      if ($smallscreen) {
         print "<p><a href=\"file:Hurricane.pqa/helptable.html\" button>Help</a><p>";
         print "<p><hr><p><a href=\"file:Hurricane.pqa/hurricane.html\">Home</a> | ";
         foreach $name (sort keys(%storm_near)) {
            if ($name) {
               print "<a href=\"/Tropical/HTML/Clip/$name.html\">$name Summary</a> | ";
               print "<a href=\"/Tropical/HTML/Clip/${name}_history.html\">History~2k</a><p>";
            }
         }
      } elsif ($wml) {
         my $accesskey = 1;
         my $accessstring;
         foreach $name (sort keys(%storm_near)) {
           if ($name) {
              if ($accesskey <= 9) {$accessstring = " accesskey=\"$accesskey\"";} else {$accessstring = "";}
              print "<p><anchor$accessstring>$name Summary<go href=\"/Tropical/HTML/WML/index.wml#s$name\"/></anchor></p>\n";
              ++$accesskey;
              if ($accesskey <= 9) {$accessstring = " accesskey=\"$accesskey\"";} else {$accessstring = "";}
              print "<p><anchor$accessstring>$name History<go href=\"/Tropical/HTML/WML/${name}_history.wml\"/></anchor></p>\n";
              ++$accesskey
           }
         }
         if ($accesskey <= 9) {$accessstring = " accesskey=\"$accesskey\"";} else {$accessstring = "";}
         print "<p><anchor$accessstring>Home<go href=\"/Tropical/HTML/WML/index.wml#home\"/></anchor></p>\n";
         ++$accesskey;
         if ($accesskey <= 9) {$accessstring = " accesskey=\"$accesskey\"";} else {$accessstring = "";}
         print "<p><anchor$accessstring>Strike Probabilities<go href=\"/Tropical/HTML/WML/strike.wml\"/></anchor></p>\n";
         ++$accesskey;
         if ($accesskey <= 9) {$accessstring = " accesskey=\"$accesskey\"";} else {$accessstring = "";}
         print "<p><anchor$accessstring>Prev<prev/></anchor></p>\n";
         ++$accesskey;
      } else {
         print '<P><BIG><A HREF="/Tropical/StrikeProb/instructions.html">Instructions</A></BIG><P>';
      }
      my $algorithm = "Slow";
      if ($gofast) {$algorithm="Fast";}
      print "\n<!-- StrikeProb.pl Version=$VERSION Algorithm=$algorithm -->\n";
      print "\n";
   }
}

################### SUBROUTINES ###############################

sub numerically { $a <=> $b;}

sub INTERPOLATE {

# Midpoint interpolation of storm data line

# TRM 1995-04-24 Made @l1,@l2,@l3 local

   my($line1,$line2) = @_;
   my(@l1) = split(/ /,$line1);
   my(@l2) = split(/ /,$line2);
   my(@l3) = @l2;
   my($s1,$s2);
   my($year,$mon,$mday,$hour,$min);
   

   $l3[$DATE] = ($l1[$DATE]+$l2[$DATE])/2;
   if ($l2[$ACTDATE] < $l1[$ACTDATE]) {$l3[$ACTDATE] = $l2[$ACTDATE]}
   else {$l3[$ACTDATE] = $l1[$ACTDATE]}
   $l3[$OBSTYPE] = "FOR";
   if ($l1[$NS] eq "S") {$s1 = -1;} else {$s1 = +1;}
   if ($l2[$NS] eq "S") {$s2 = -1;} else {$s2 = +1;}
   $l3[$LATITUDE] =  &INTERP($l1[$DATE],$l2[$DATE],$l3[$DATE],$l1[$LATITUDE]*$s1,$l2[$LATITUDE]*$s2);
   if ($l1[$EW] eq "E") {$s1 = -1;} else {$s1 = +1;}
   if ($l2[$EW] eq "E") {$s2 = -1;} else {$s2 = +1;}
   if ($l1[$EW] ne $l2[$EW]) {
      $l2[$LONGITUDE] -= $s2*360;
      $s2 *= -1; 
   }
   $l3[$LONGITUDE] =  &INTERP($l1[$DATE],$l2[$DATE],$l3[$DATE],$l1[$LONGITUDE]*$s1,$l2[$LONGITUDE]*$s2);
   if ($l3[$LATITUDE] >= 0) {$l3[$NS] = "N"} else {$l3[$NS] = "S"}
   if ($l3[$LONGITUDE] >= 0) {$l3[$EW] = "W"} else {$l3[$EW] = "E"}
   if ($l3[$LATITUDE] < 0) {$l3[$LATITUDE] = -$l3[$LATITUDE];}
   if ($l3[$LONGITUDE] < 0) {$l3[$LONGITUDE] = -$l3[$LONGITUDE];}
   $l3[$TNAME] = $l1[$TNAME];
   ($year,$mon,$mday,$hour,$min) = &YEAR2DATE($l3[$DATE]);
   $l3[$TYEAR] = $year;
   $l3[$TMONTH] = $mon;
   $l3[$TDAY] = $mday;
   $l3[$TTIME] = $hour.".".$min;
   $l3[$TWIND] = int(&INTERP($l1[$DATE],$l2[$DATE],$l3[$DATE],$l1[$TWIND],$l2[$TWIND])); 
   $l3[$TGUST] = int(&INTERP($l1[$DATE],$l2[$DATE],$l3[$DATE],$l1[$TGUST],$l2[$TGUST])); 

   join(" ",@l3);

}

sub INTERP {

   # Do a linear interpolation

   my($t1,$t2,$t,$v1,$v2) = @_;

   $v1+(($t-$t1)/($t2-$t1))*($v2-$v1);

}

@gammacof = (76.18009173,
             -86.50532033,
             24.01409822,
             -1.231739516,
             0.120858003e-2,
             -0.536382e-5);

sub GAMMALN {

   # ln of the gamma function.

   my($z) = @_;
   my($answer,$x,$y,$tmp,$j,$ser);
   $y = $z;
   $x = $z;
   $tmp = $x+5.5;
   $tmp -= ($x+0.5)*log($tmp);
   $ser = 1.000000000190015;
   for ($j=0; $j<=5; $j++) {$ser += $gammacof[$j]/(++$y);}
   -$tmp+log(2.5066282746310005*$ser/$x);
}

sub GAMMA {

   # The gamma function

   my($z) = @_;
   exp(&GAMMALN($z));
}


sub LNFACTORIAL {

   # log factorial function
   # Positive integers only

   my($x) = @_;

   if ($x > 12) {&GAMMALN($x+1);}
   else {
      log(&FACTORIAL($x));
   }
}

sub FACTORIAL {

   # factorial function
   # Positive integers only

   my($x) = @_;

   if ($x > 12) {&GAMMA($x+1);}
   else {
      if ($x <= 1) { 1;} else {$x*&FACTORIAL($x-1);}
   }
}

sub GAMMP {

   # Incomplete Gamma Function, Gamma(a)*P(a,x)
   # $ x must be positive or zero, $a must be positve

   my($a,$x) = @_;
   my($gammp);

   if ($x < ($a + 1)) {
      $gammp = &GSER($a,$x);
   } else {
      $gammp = 1 - &GCF($a,$x);
   }

   $gammp;
   
}

$ITMAX = 100;
$EPS = 3.0e-7;

sub GSER {

   # Incomplete Gamma Function, P(a,x)
   # $ x must be positive or zero, $a must be positve

   my($a,$x) = @_;
   my($gln,$gser,$ap,$del,$sum,$n,$absdel,$abssum);

   $gln = &GAMMALN($a);
   if ($x <= 0.0) {
      $gser = 0.0;    # Input Error
   } else {
      $ap = $a;
      $del = $sum = 1.0/$a;
      for ($n=1; $n<=$ITMAX; $n++) {
         ++$ap;
         $del *= $x/$ap;
         $sum += $del;
         $absdel = $del;
         if ($absdel < 0.0) {$absdel = -$absdel;}
         $abssum = $sum;
         if ($abssum < 0.0) {$abssum = -$abssum;}
         if ($absdel < $abssum*$EPS) {
            last;
         }
      }
      $gser = $sum * exp(-$x+$a*log($x)-$gln)
   }
   $gser;
}

$FPMIN = 1.0e-30;

sub GCF {

   # Incomplete Gamma Function, Q(a,x)
   # $ x must be positive or zero, $a must be positve

   my($a,$x) = @_;
   my($gln,$b,$c,$d,$h,$i,$an,$absd,$del,$absdel1);

   $gln = &GAMMALN($a);
   $b = $x + 1.0 - $a;
   $c = 1.0/$FPMIN;
   $d = 1.0/$b;
   $h = $d;
   for ($i=1; $i<=$ITMAX; $i++) {
      $an = -$i*($i-$a);
      $b += 2.0;
      $d = $an*$d+$b;
      $absd = $d;
      if ($absd < 0.0) {$absd = -$absd;}
      if ($absd < $FPMIN) {$d = $FPMIN;}
      $c = $b + $an/$c;
      $absc = $c;
      if ($absc < 0.0) {$absc = -$absc;}
      if ($absc < $FPMIN) {$c = $FPMIN;}
      $d = 1.0/$d;
      $del = $d*$c;
      $h *= $del;
      $absdel1 = $del - 1.0;
      if ($absdel1 < 0) {$absdel1 = -$absdel1;}
      if ($absdel1 < $EPS) {last;}
   }
   exp(-$x+$a*log($x)-$gln)*$h;

}

sub ACOS {

# Compute the arc cosine
# Input: Cosine (forced to -1 <= cosine <= 1)
# Output: Angle in radians (0<angle<pi)

   my($cosine) = @_;
   my($argument);
   $argument = 1.0-$cosine**2;
   if ($argument<=0) {
      if ($cosine > 0) {0;} else {$pi;}
   } else {
      $piover2 - atan2($cosine/(sqrt($argument)),1);
   }
}

sub tan { sin($_[0]) / cos($_[0]) }

sub GCDISTANCE {

# Compute Great Circle distance between positions
# Input:  lat1, lon1, lat2, lon, all in decimal degrees
# Output: Distance in Nautical Miles

   my($latitude1,$longitude1,$latitude2,$longitude2) = @_;

   &ACOS(sin($latitude1*$dtor)*sin($latitude2*$dtor) + 
         cos($latitude1*$dtor)*cos($latitude2*$dtor) *
         cos(($longitude1-$longitude2)*$dtor))*60/$dtor;
}

sub GCCOURSE {

   # Compute great circle course

   my($latitude1,$longitude1,$latitude2,$longitude2) = @_;

   $latitude1 *= $dtor;
   $latitude2 *= $dtor;
   $longitude1 *= $dtor;
   $longitude2 *= $dtor;

   my $dlon = $longitude1-$longitude2;
   atan2(sin($dlon),
         cos($latitude1)*tan($latitude2)-sin($latitude1)*cos($dlon))/$dtor

}

sub DEG2DIRECTION {
   my ($c) = @_;
   my $tcourse;

   while ($c > 180) {$c -= 360;}
   while ($c < -180) {$c += 360;}

   if ($c > -11.25 && $c <= +11.25) {$tcourse = "N";}
   if ($c > +11.25 && $c <= +33.75) {$tcourse = "NNE";}
   if ($c > +33.75 && $c <= +56.25) {$tcourse = "NE";}
   if ($c > +56.25 && $c <= +78.75) {$tcourse = "ENE";}
   if ($c > +78.75 && $c <= +101.75) {$tcourse = "E";}
   if ($c > +101.75 && $c <= +123.75) {$tcourse = "ESE";}
   if ($c > +123.75 && $c <= +146.25) {$tcourse = "SE";}
   if ($c > +146.25 && $c <= +168.75) {$tcourse = "SSE";}
   if ($c > +168.75 || $c <= -168.75) {$tcourse = "S";}
   if ($c <= -146.25 && $c > -168.75) {$tcourse = "SSW";}
   if ($c <= -123.75 && $c > -146.25) {$tcourse = "SW";}
   if ($c <= -101.25 && $c > -123.75) {$tcourse = "WSW";}
   if ($c <= -78.75 && $c > -101.25) {$tcourse = "W";}
   if ($c <= -56.25 && $c > -78.75) {$tcourse = "WNW";}
   if ($c <= -33.75 && $c > -56.25) {$tcourse = "NW";}
   if ($c <= -11.25 && $c > -33.75) {$tcourse = "NNW";}
   #print "$c $tcourse\n";

   $tcourse;
}


sub ISNUMBER {

   # Is the argument a number, e.g. 1.0, .1, 1, 1e4, 1.0e4, etc.?
   # Returns true or false.

   $_[0] =~ /^\s*(\+|-|)(\d+|\d+\.\d*|\d*\.\d+)(e\d+|)\s*$/i;
}

sub DOSPROB {

   # compute the strike probability for a lat and lon given the forecast
   # data in @line.  Returns with a string of data, outline, including
   # name, date, time, wind, lat, lon, and strike probabilities.

   my($lat,$lon,@line) = @_;
   my($forlat,$forlon,$distance,$tdiff,$forwind,$name,$year,$month,$day,$hour,$minute,$ForecastError,$good,$dist,$k1,$k2,$strike,$outline);

   $forlat = $line[$LATITUDE];
   $forlon = $line[$LONGITUDE];
   if ($line[$NS] eq 'S') {$forlat = -$forlat;}
   if ($line[$EW] eq 'E') {$forlon = -$forlon;}
   $distance = &GCDISTANCE($lat,$lon,$forlat,$forlon);
   $tdiff = ($line[$DATE] - $line[$ACTDATE])/$onehour;  # in hours
   $forwind = $line[$TWIND];
   $good = 0;
   $outline = "";
   if ($tdiff > 0 && $distance < ($tdiff*60+$STORMSIZE[$#STORMSIZE])) {
      $name   = substr($line[$TNAME],0,length($line[$TNAME])-3);
      $year   = $line[$TYEAR];
      $month  = $line[$TMONTH];
      if (length($month) == 1) {$month = "0".$month;}
      $day = $line[$TDAY];
      if (length($day) == 1) {$day = "0".$day;}
      $line[$TTIME] =~ /(\d+)\.(\d*)/i;
      $hour   = $1;
      if (length($hour) == 1) {$hour = "0".$hour;}
      $minute = $2;
      if (length($minute) == 1) {$minute = "0".$minute;}
      $date = $line[$DATE];
      $outline = "$name $year $month $day $hour $minute $date $forwind $forlat $forlon $tdiff";
      #$ForecastError = &FORECASTERROR($tdiff,$INTERCEPT);
      if ($actdate{$name} eq "" || $actdate{$name} > $line[$ACTDATE]) {$actdate{$name} = $line[$ACTDATE];}
      $good = 0;
      for ($i=0; $i<=$#STORMSIZE; ++$i) {
         $STORMSIZE = $STORMSIZE[$i];
         $strike = &SPROB($distance,$STORMSIZE,$tdiff);
         $outline = $outline." $strike";
         if ($strike >= 0.005) {$good=1;}  # Don't show storms with 0 probability
      }
   }
   ($good,$outline);
}

sub SPROB {
   # computes the probability that a storm will be within $STORMSIZE nmi
   # of a point $distance nmi from the forecast position given that the
   # forecast is $tdiff hours old.  Output in percent: range 0-1.

# Compute K(R,S,E[,error]). The optional fourth argument controls
# termination of the summation (default: 1.0e-6).
# Extreme cases:
# R-S <= -4 => K=1
# R-S >=  4 => K=0
# were determined numerically with termination control set to 1.0e-6.
# Note: this subroutine is equivalent to the old SPROB.
# Last revision: July 18, 2006 (ftm)
  # Variables:
  my ($R,$S,$tdiff,$eps) = @_;
  if (!$eps) { $eps = 1e-6; } # Default convergence test is 1e-6
  if ($tdiff < 0) {$tdiff = -$tdiff;}   # absolute value
  if ($R < 0) {$R = 0;}
  $E = &FORECASTERROR($tdiff,$INTERCEPT);
  my $R_E = $R/$E;
  my $S_E = $S/$E;
  # Test for extreme cases ("4" was determined using 1e-6 accuracy):
  my $R_E_S_E = $R_E-$S_E;
  if ($R_E_S_E <= -4) { return 1; }
  elsif ($R_E_S_E >= 4) { return 0; }
  # Initialize:
  my ($n,$x,$y,$Tn_1,$Tn,$Un_1,$Un,$sum,$x_n,$y_n);
  $n = 1;
  $x = $R_E*$R_E;
  $y = $S_E*$S_E;
  $Tn_1 = 1;
  $Un_1 = exp($y)-1;
  $sum = $Un_1;
  # Recursion:
  do {
    $x_n = $x/$n;
    $y_n = $y/$n;
    $Tn = $x_n*$y_n*$Tn_1;
    $Un = $x_n*$Un_1-$Tn;
    $sum += $Un;
    $Tn_1 = $Tn;
    $Un_1 = $Un;
    $n++;
  } until (abs($Un) < $eps*$sum);
  return exp(-$x-$y)*$sum;
}

sub AVGSPROBALT1 {

   # Computes the integral from 0 to S of SPROB(r,S,T) 2 pi r dr divided
   # by pi s^2.

   # This assumes that the storm does not evolve from B given that B is
   # true

   my($STORMSIZE,$tdiff,$distance,$reltdiff,$aprob,$bprob) = @_;
   my($ForecastError,$strike,$s_over_e,$i,
         $term,$a1,$a1ln,$term1,$term2,$j,$s_over_e2);
   my $answer;

   # compute average strike probability

   if ($tdiff < 0) {$reltdiff = -$reltdiff;}   # absolute value
   $ForecastError = &FORECASTERROR($reltdiff,0);
   if ($ForecastError < $STORMSIZE/30.0) {
      # if S/E gets too big the series converges *very* slowly, so limit this.
      #if ($distance <= $STORMSIZE-2*$ForecastError) {return 1;}
      #if ($distance >= $STORMSIZE+2*$ForecastError) {return 0;}
      $ForecastError = $STORMSIZE/30.0;
   }
   $s_over_e = $STORMSIZE/$ForecastError;
   $s_over_e2 = $s_over_e**2;
   if ($s_over_e <= 0) {
      $answer = 0.0;
   } else {
      $strike = 0.0;
      $i=0;
      do {
         $term = ((&GAMMP($i+1,$s_over_e2))**2);
         $strike += $term;
         ++$i;
      } until (($term*100000 < $strike) || ($i > 300));

      $answer = $strike/($s_over_e2);
   }

#   #if ($aprob < $bprob) {$answer = $answer * $aprob/$bprob;} # Bayes theorem: symmetry condition
#
#   # the condition that P(A|BF)<P(A)/P(B) and
#   # the condition that P(A|BF)> (P(A)+P(B)-1)/P(B) are violated when the forecasts are
#   # too far apart.  To deal with this, use P(A|BF)=P(A) when the conditions are violated.
#   # To make this smooth set P(A|BF)=P(A) when P(A)+P(B)>1 and K<P(A).
#
#   if ( (($aprob+$bprob) > 1) && ($answer < $aprob) ) {$answer = $aprob;} # do NOT use >=1 !!!
#
#   # Check that the approximation does not violate obvious conditions. 
#   if ($bprob > 0) {
#      $testab = $aprob/$bprob;
#      $testab2 = ($aprob+$bprob-1)/$bprob;
#   } else {
#      $testab=1;
#      $testab2=0;
#   }
#   if ($testab > 1) {$testab = 1;} elsif ($testab < 0) {$testab = 0;}
#   if ($testab2 > 1) {$testab2 = 1;} elsif ($testab2 < 0) {$testab2 = 0;}
#   if ($answer > $testab) {
#      #print "<!-- OVER: $answer $testab $aprob $bprob -->\n";
#      if ($testab > 0.99) {
#         $answer = $testab;   # correct roundoff error
#      } else {
#         $answer = $aprob;    # method failed revert to P(A|B)=P(A)
#      }
#   }
#   if ($answer < $testab2) {
#      #print "<!-- UNDER: $answer $testab2 $aprob $bprob -->\n";
#      if ($testab2 < 0.01) {
#         $answer = $testab2;  # correct roundoff error
#      } else {
#         $answer = $aprob;    # method failed revert to P(A|B)=P(A)
#      }
#   }

   $answer;
}

sub AVGINTRK4 {

   # Do the integral using Runge-Kutta with auto step scaling.

   my($i,$r_over_e,$r_over_e2,$s_over_e,$s_over_e2) = @_;

   my $lowbound = $r_over_e - $s_over_e;
   if ($lowbound < 0) {$lowbound = 0;}
   my $upbound = ($r_over_e+$s_over_e);
   @RK4PARAMS = ($i,&LNFACTORIAL($i),$r_over_e,$s_over_e,$r_over_e2,$s_over_e2);   # must have global scope
   my $t = $lowbound;
   my $nsteps = 20;
   my $dt = ($upbound-$lowbound)/$nsteps;
   my (@y,@errors,$step);
   $errors[0] = 0.1;  # Should be smaller, but it is too slow
   $y[0] = 0.0;

   while ($t+$dt < $upbound) {
      ($t, $dt, @y) = rk4_auto(\@y,\&dydt,$t,$dt,\@errors);
   }
   $dt = $upbound - $t;
   if ($dt) {($t, @y) = &rk4( \@y, \&dydt, $t, $dt );}

   $y[0];
}


sub dydt {
   my ($t,@y) = @_;
   my ($rho,$rho2,$term);
   $rho = $t;
   $rho2 = $rho*$rho;
   my ($i,$LNFACTI,$r_over_e,$s_over_e,$r_over_e2,$s_over_e2) = @RK4PARAMS;

   if ($rho <= 0) {
      $term = 0;
   } else {
      $term = exp(log($rho)*(2*$i+1) - $LNFACTI - $rho2);
   }
   if (($rho+$r_over_e) <= $s_over_e || $r_over_e == 0) {
      $term *= $pi;
   } else {
      $term *= &ACOS(($r_over_e2-$s_over_e2+$rho2)/(2*$r_over_e*$rho));
   }
   return ($term);
}

sub AVGINT {

    # Trapezoidal rule integration

    my($i,$r_over_e,$r_over_e2,$s_over_e,$s_over_e2) = @_;
    my($rho,$rho2,$h,$integral,$step,$lowbound,$term,$nterms,$lnfacti,$facti);
    
    $integral = 0;
    $nterms = 20;
    $lowbound = $r_over_e - $s_over_e;
    if ($lowbound < 0) {$lowbound = 0;}
    $step = ($r_over_e+$s_over_e-$lowbound)/$nterms;
    $lnfacti = &LNFACTORIAL($i);
    for ($h=0; $h<=$nterms; ++$h) {
       $rho = $step*$h + $lowbound;
       $rho2 = $rho*$rho;
       if ($rho <= 0) {
          $term = 0;
       } else {
          $term = exp(log($rho)*(2*$i+1) - $lnfacti - $rho2);
       }
       if (($rho+$r_over_e) <= $s_over_e) {
	  $term *= $pi;
       } else {
          $term *= &ACOS(($r_over_e2-$s_over_e2+$rho2)/(2*$r_over_e*$rho));
       }
       if ($h == 0 || $h == $nterms) {
          $integral += $term/2;
       } else {
          $integral += $term;
       }
    }
    $integral * $step;

}

sub AVGINT2 {

    # Trapezoidal rule integration with variable spacing

    my($i,$r_over_e,$r_over_e2,$s_over_e,$s_over_e2) = @_;
    my($rho,$rho2,$h,$integral,$step,$lowbound,$upbound,$term,$nterms,$lnfacti,$facti);
    my($lastterm,$lastrho);
    
    $integral = 0;
    $nterms = 30;  # if you change this, change $step below too
    $lowbound = $r_over_e - $s_over_e;
    $upbound = $r_over_e + $s_over_e;
    if ($lowbound < 0) {$lowbound = 0;}
    # use a variable step size which is small for
    # small rho but gets bigger for larger rho
    # 465 is 1+2+3+4+...+29+30
    $step = ($upbound-$lowbound)/465.0;  # 465 only works for nterms = 30
    $lnfacti = &LNFACTORIAL($i);
    $lastrho = $lowbound;
    $lastterm = 0;
    for ($h=0; $h<=$nterms; ++$h) {
       $rho = $step*$h + $lastrho;
       $rho2 = $rho*$rho;
       if ($rho <= 0) {
          $term = 0;
       } else {
          $term = exp(log($rho)*(2*$i+1) - $lnfacti - $rho2);
       }
       if (($rho+$r_over_e) <= $s_over_e || $r_over_e == 0) {
	  $term *= $pi;
       } else {
          $term *= &ACOS(($r_over_e2-$s_over_e2+$rho2)/(2*$r_over_e*$rho));
       }
       if ($h) {$integral += (0.5*($term+$lastterm)*($rho-$lastrho));}
       $lastrho = $rho;
       $lastterm = $term;
    }
    #print "check these are equal $rho $upbound\n";
    $integral;

}


sub AVGSPROB {


   my($STORMSIZE,$tdiff,$distance,$reltdiff,$aprob,$bprob) = @_;
   my($ForecastError,$strike,$s_over_e,$r_over_e,$r_over_e2,$i,$testab,$testab2,$s_over_ee,$s_over_ee2,
         $term,$a1,$a1ln,$term1,$term2,$j,$s_over_e2,$integral,$answer);

   # compute average strike probability over a cicle of radius $STORMSIZE
   # displaced from the forecast position by $distance.

   if ($reltdiff < 0) {$reltdiff = -$reltdiff;}   # absolute value
   $ForecastError = &FORECASTERROR($reltdiff,0);
   if ($ForecastError < $STORMSIZE/30.0) {
      # if S/E gets too big the series converges *very* slowly, so limit this.
      #if ($distance <= $STORMSIZE-2*$ForecastError) {return 1;}
      #if ($distance >= $STORMSIZE+2*$ForecastError) {return 0;}
      $ForecastError = $STORMSIZE/30.0;
   }

   $s_over_e = $STORMSIZE/$ForecastError;
   $s_over_e2 = $s_over_e**2;
   $r_over_e = $distance/$ForecastError;
   $r_over_e2 = $r_over_e**2;
   if ($s_over_e <= 0) {
      $answer = 0.0;
   } else {
      $strike = 0.0; 
      $i=0;
      do { 
         if ($r_over_e <= 0) {
            $term = $piover2*(&GAMMP($i+1,$s_over_e2)**2); # piover2 to cancel below factor
         } else {
            #$term = &AVGINT($i,$r_over_e,$r_over_e2,$s_over_e,$s_over_e2) * &GAMMP($i+1,$s_over_e2);
            $term = &AVGINT2($i,$r_over_e,$r_over_e2,$s_over_e,$s_over_e2) * &GAMMP($i+1,$s_over_e2);
            # to use AVGINTRK4, you must uncomment the "require RungeKutta.pl" at the top of this file
            #$term = &AVGINTRK4($i,$r_over_e,$r_over_e2,$s_over_e,$s_over_e2) * &GAMMP($i+1,$s_over_e2);
         }
         $strike += $term;
         ++$i;
      } until (($term*1000000 < $strike) || ($i > 1000) );

      $answer = $strike/($s_over_e2*$piover2);
   }

#   #if ($aprob < $bprob) {$answer = $answer * $aprob/$bprob;} # Bayes theorem: symmetry condition
#
#   # the condition that P(A|BF)<P(A)/P(B) and
#   # the condition that P(A|BF)> (P(A)+P(B)-1)/P(B) are violated when the forecasts are
#   # too far apart.  To deal with this, use P(A|BF)=P(A) when the conditions are violated.
#   # To make this smooth set P(A|BF)=P(A) when P(A)+P(B)>1 and K<P(A).
#
#   if ( (($aprob+$bprob) > 1) && ($answer < $aprob) ) {$answer = $aprob;} # do NOT use >=1 !!!
#
#   # Check that the approximation does not violate obvious conditions.
#   if ($bprob > 0) {
#      $testab = $aprob/$bprob;
#      $testab2 = ($aprob+$bprob-1)/$bprob;
#   } else {
#      $testab=1;
#      $testab2=0;
#   }
#   if ($testab > 1) {$testab = 1;} elsif ($testab < 0) {$testab = 0;}
#   if ($testab2 > 1) {$testab2 = 1;} elsif ($testab2 < 0) {$testab2 = 0;}
#   if ($answer > $testab) {
#      #print "<!-- OVER: $answer $testab $aprob $bprob -->\n";
#      if ($testab > 0.99) {
#         $answer = $testab;   # correct roundoff error
#      } else {
#         $answer = $aprob;    # method failed revert to P(A|B)=P(A)
#      }
#   }
#   if ($answer < $testab2) {
#      #print "<!-- UNDER: $answer $testab2 $aprob $bprob -->\n";
#      if ($testab2 < 0.01) {
#         $answer = $testab2;  # correct roundoff error
#      } else {
#         $answer = $aprob;    # method failed revert to P(A|B)=P(A)
#      }
#   }

   $answer;

}

sub AVGSPROBALTint {

   my($STORMSIZE,$tdiff,$distance,$reltdiff,$aprob,$bprob) = @_;

   # Assumes that the forecast is shifted to B when B is given.  This
   # implementation is very approximate ... should be doing an integral.
   # This approximation is so bad, it is pretty much useless.

   #&SPROB($distance,$STORMSIZE,$tdiff);

   # This is a better approximation to the integral but it is very slow.

   my($d1) = $distance;
   my($d2) = $distance + $STORMSIZE;
   my($d3) = $distance - $STORMSIZE;
   my($d4) = sqrt($distance**2 + $STORMSIZE**2);
   if ($d3 < 0.) {$d3 = -$d3;}

   #(&SPROB($d1,$STORMSIZE,$reltdiff) +
   #&SPROB($d2,$STORMSIZE,$reltdiff) +
   #&SPROB($d3,$STORMSIZE,$reltdiff) +
   #2*&SPROB($d4,$STORMSIZE,$reltdiff))/5.0;

   my($d5) = sqrt( ($distance-0.3536*$STORMSIZE)**2 + (0.3536*$STORMSIZE)**2 );
   my($d6) = sqrt( ($distance+0.3536*$STORMSIZE)**2 + (0.3536*$STORMSIZE)**2 );

   (&SPROB($d1,$STORMSIZE,$reltdiff) +
   &SPROB($d2,$STORMSIZE,$reltdiff) +
   &SPROB($d3,$STORMSIZE,$reltdiff) +
   2*&SPROB($d4,$STORMSIZE,$reltdiff) +
   2*&SPROB($d5,$STORMSIZE,$reltdiff)+
   2*&SPROB($d6,$STORMSIZE,$reltdiff))/9.0;


}

sub AVGSPROBALT3 {

   my($STORMSIZE,$tdiff,$distance,$reltdiff,$aprob,$bprob) = @_;

   # This version assumes that B has no information on A.  This is the
   # fastest version.  It is perhaps the best as well, since a reforcast
   # would be required to really know what the effect of A true B is on A.

   $aprob;   # Assumes that B has no information for A

}

sub AVGSPROBALT4 {
     
   my($STORMSIZE,$tdiff,$distance,$reltdiff,$aprob,$bprob) = @_;
     
   # This version assumes that B has no information on A, but
   # the error on A is reduced to the time diff bewtween A
   # and B.
     
   if ($reltdiff < 0) {$reltdiff = -$reltdiff;}   # absolute value
   my $answer = &SPROB($distance,$STORMSIZE,$reltdiff);
   # Check that the approximation does not violate obvious conditions
   if ($bprob > 0) {
      $testab = $aprob/$bprob;
      $testab2 = ($aprob+$bprob-1)/$bprob;
   } else {
      $testab=1;
      $testab2=0;
   }
   if ($testab > 1) {$testab = 1;} elsif ($testab < 0) {$testab = 0;}
   if ($testab2 > 1) {$testab2 = 1;} elsif ($testab2 < 0) {$testab2 = 0;}
   if ($answer > $testab) {$answer = $testab;}
   if ($answer < $testab2) {$answer = $testab2;}
   $answer;
  
}

                                                                                                                              
sub FORECASTERROR {

   my($tdiff,$intercept) = @_;
   if ($tdiff < 0) {$tdiff = -$tdiff;} # absolute value
   $intercept + $PSLOPE*($tdiff**$EPOWER);

}

sub YEAR2DATE {

   # converts a decimal year to a date (year, month, day, hour, minute)
   # The decimal year is assumed to be computed as 
   # year + doy/366.0 where year is the year (like 1996) and doy is the
   # day of year (like 215.5).

   my($date) = @_;
   my($foryr70,$forsec,$formin,$forhour,$formday,$formon,$foryear,$forwday,$foryday,$forisdst);

   $foryr70 = int($date)-1970;
   # This will fail in 2100 since it is not a leap year
   $sec_since_1970 = $foryr70*31536000 + int((int($date)-1-1972)/4)*86400 + ($date-int($date))*31622400 + 0.5;
   ($forsec,$formin,$forhour,$formday,$formon,$foryear,$forwday,$foryday,$forisdst) = gmtime($sec_since_1970);
   $formon += 1;
   if ($foryear>=96) {$foryear+=1900;} else {$foryear+=2000;}
   if ($formon < 10) {$formon = "0".$formon;}
   if ($formday < 10) {$formday = "0".$formday;}
   if ($forhour < 10) {$forhour = "0".$forhour;}
   if ($formin < 10) {$formin = "0".$formin;}
   ($foryear,$formon,$formday,$forhour,$formin);
}

sub CITY {

   my($city) = @_;
   my($lat,$lon,$locstring);

   if    ($city =~ /^(HNL|HONOLULU)/i)     {$lat=21.21; $lon=157.56; $locstring="Honolulu, Oahu";}
   elsif ($city =~ /^(ACAPULCO)/i) {$lat=16.50; $lon=99.56; $locstring = "Acapulco, Mexico";}
   elsif ($city =~ /^(AGANA)/i) {$lat=13.29; $lon=-144.48; $locstring = "Agana, Guam";}
   elsif ($city =~ /^(ANDERSEN AFB|ANDERSENAFB)/i) {$lat=13.34; $lon=-144.55; $locstring = "Andersen AFB, Guam";}
   elsif ($city =~ /^(APALACHICOLA)/i) {$lat=29.43; $lon=85.01; $locstring = "Apalachicola, Florida";}
   elsif ($city =~ /^(ATLANTIC CITY|ATLANTICCITY)/i) {$lat=39.45; $lon=74.57; $locstring = "Atlantic City, New Jersey";}
   elsif ($city =~ /^(AUCKLAND)/i) {$lat=-37.01; $lon=-174.48; $locstring = "Auckland, New Zealand";}
   elsif ($city =~ /^(BAGUIO CITY|BAGUIOCITY)/i) {$lat=16.25; $lon=-120.36; $locstring = "Baguio City, Philippines";}
   elsif ($city =~ /^(BALTIMORE)/i) {$lat=39.11; $lon=76.40; $locstring = "Baltimore, Maryland";}
   elsif ($city =~ /^(BANGKOK)/i) {$lat=13.23; $lon=-100.36; $locstring = "Bangkok, Thailand";}
   elsif ($city =~ /^(BARAHONA)/i) {$lat=18.13; $lon=71.06; $locstring = "Barahona, Dominican Republic";}
   elsif ($city =~ /^(BATON ROUGE|BATONROUGE)/i) {$lat=30.30; $lon=91.10; $locstring = "Baton Rouge, Louisiana";}
   elsif ($city =~ /^(BEAUMONT)/i) {$lat=30.00; $lon=94.00; $locstring = "Beaumont, Texas";}
   elsif ($city =~ /^(BEEF ISLAND|BEEFISLAND)/i) {$lat=18.27; $lon=64.32; $locstring = "Beef Island, Tortola, BVI";}
   elsif ($city =~ /^(BILOXI)/i) {$lat=30.24; $lon=88.55; $locstring = "Biloxi, Mississsippi";}
   elsif ($city =~ /^(BOCA RATON|BOCARATON)/i) {$lat=26.23; $lon=80.09; $locstring = "Boca Raton, Florida";}
   elsif ($city =~ /^(BOMBAY)/i) {$lat=19.07; $lon=-72.51; $locstring = "Bombay, India";}
   elsif ($city =~ /^(BOSTON)/i) {$lat=42.00; $lon=71.00; $locstring = "Boston, Massachusetts";}
   elsif ($city =~ /^(BRIDGETOWN)/i) {$lat=13.06; $lon=59.36; $locstring = "Bridgetown, Barbados";}
   elsif ($city =~ /^(BROOME)/i) {$lat=-17.56; $lon=-122.10; $locstring = "Broome, Australia";}
   elsif ($city =~ /^(BROWNSVILLE)/i) {$lat=25.54; $lon=97.30; $locstring = "Brownsville, Texas";}
   elsif ($city =~ /^(CABO|CABO SAN LUCAS)/i) {$lat=22.52; $lon=109.54; $locstring = "Cabo San Lucas, Mexico";}
   elsif ($city =~ /^(CAIRNS)/i) {$lat=-16.53; $lon=-145.45; $locstring = "Cairns, Australia";}
   elsif ($city =~ /^(CAMPECHE)/i) {$lat=19.51; $lon=90.33; $locstring = "Campeche, Mexico";}
   elsif ($city =~ /^(CANCUN)/i) {$lat=21.09; $lon=86.45; $locstring = "Cancun, Mexico";}
   elsif ($city =~ /^(CARACAS)/i) {$lat=10.36; $lon=66.59; $locstring = "Caracas, Venezuela";}
   elsif ($city =~ /^(CEBU)/i) {$lat=10.19; $lon=-123.59; $locstring = "Cebu, Philippines";}
   elsif ($city =~ /^(CHARLSTN|CHARLESTON)/i) {$lat=32.468; $lon=79.552; $locstring = "Charleston, S. Carolina";}
   elsif ($city =~ /^(COCKATOO ISLAND|COCKATOOISLAND)/i) {$lat=-16.30; $lon=-123.36; $locstring = "Cockatoo Island, Australia";}
   elsif ($city =~ /^(COOPERCITY|COOPER CITY)/i) {$lat=26.04; $lon=80.16; $locstring = "Cooper City, Florida";}
   elsif ($city =~ /^(CORPUS CHRISTI|CORPUSCHRISTI)/i) {$lat=27.47; $lon=97.26; $locstring = "Corpus Christi, Texas";}
   elsif ($city =~ /^(DAMPIER)/i) {$lat=-20.41; $lon=-116.41; $locstring = "Dampier, Australia";}
   elsif ($city =~ /^(DARWIN)/i) {$lat=-12.24; $lon=-130.52; $locstring = "Darwin, Australia";}
   elsif ($city =~ /^(DAVAO)/i) {$lat=7.07; $lon=-125.39; $locstring = "Davao, Philippines";}
   elsif ($city =~ /^(DERBY)/i) {$lat=-17.18; $lon=-123.37; $locstring = "Derby, Australia";}
   elsif ($city =~ /^(DIEGO GARCIA|DIEGOGARCIA)/i) {$lat=-7.18; $lon=-72.24; $locstring = "Diego Garcia";}
   elsif ($city =~ /^(FREEPORT-BAH)/i) {$lat=26.33; $lon=78.42; $locstring = "Freeport, Grand Bahama";}
   elsif ($city =~ /^(FREEPORT)/i) {$lat=28.54; $lon=95.24; $locstring = "Freeport, Texas";}
   elsif ($city =~ /^(GALVESTON)/i) {$lat=29.17; $lon=94.48; $locstring = "Galveston, Texas";}
   elsif ($city =~ /^(GEORGE TOWN|GEORGETOWN)/i) {$lat=23.30; $lon=75.46; $locstring = "George Town, Exuma";}
   elsif ($city =~ /^(GRANDCAYMAN|GRAND CAYMAN)/i) {$lat=19.17; $lon=81.21; $locstring = "Airport, Grand Cayman Island";}
   elsif ($city =~ /^(HALIFAX)/i) {$lat=44.40; $lon=63.35; $locstring = "Halifax, Nova Scotia";}
   elsif ($city =~ /^(HAMILTON)$/i) {$lat=32.22; $lon=64.40; $locstring = "Hamilton, Bermuda";}
   elsif ($city =~ /^(HAMILTONNZ)/i) {$lat=-37.51; $lon=-175.20; $locstring = "Hamilton, New Zealand";}
   elsif ($city =~ /^(HARTFORD)/i) {$lat=41.56; $lon=72.41; $locstring = "Hartford, Connecticut";}
   elsif ($city =~ /^(HAVANA)/i) {$lat=23.07; $lon=82.25; $locstring = "Havana, Cuba";}
   elsif ($city =~ /^(HILO)/i) {$lat=19.43; $lon=155.04; $locstring = "Hilo, Hawaii";}
   elsif ($city =~ /^(HKONG|HONG KONG)/i) {$lat=22.30; $lon=-114.18; $locstring = "Hong Kong";}
   elsif ($city =~ /^(HOUSTON)/i) {$lat=29.45; $lon=95.25; $locstring = "Houston, Texas";}
   elsif ($city =~ /^(ILOILO)/i) {$lat=10.42; $lon=-122.34; $locstring = "Iloilo, Philippines";}
   elsif ($city =~ /^(ISHIGAKI)/i) {$lat=24.20; $lon=-124.11; $locstring = "Ishigaki, Japan";}
   elsif ($city =~ /^(IWAKUNI)/i) {$lat=34.09; $lon=-132.14; $locstring = "Iwakuni, Japan";}
   elsif ($city =~ /^(JACKSONVILLE)/i) {$lat=30.03; $lon=81.42; $locstring = "Jacksonville, Florida";}
   elsif ($city =~ /^KARRATHA/i) {$lat=-20.46; $lon=-116.47; $locstring = "Karratha, Australia";}
   elsif ($city =~ /^(KENNEDY SPACE CENTER|KSC)/i) {$lat=28.62; $lon=80.68; $locstring = "Kennedy SC, Florida";}
   elsif ($city =~ /^(KAGOSHIMA)/i) {$lat=31.33; $lon=-130.33; $locstring = "Kagoshima, Japan";}
   elsif ($city =~ /^(KADENA)/i) {$lat=26.21; $lon=-127.46; $locstring = "Kadena AB, Okinawa";}
   elsif ($city =~ /^(KAHULUI)/i) {$lat=20.54; $lon=156.26; $locstring = "Kahului, Maui";}
   elsif ($city =~ /^(KANEOHE)/i) {$lat=21.27; $lon=157.46; $locstring = "Kaneohe, Oahu";}
   elsif ($city =~ /^(KEY WEST|KEYWEST)/i) {$lat=24.33; $lon=81.45; $locstring = "Key West, Florida";}
   elsif ($city =~ /^(KINGSTON)/i) {$lat=17.58; $lon=76.48; $locstring = "Kingston, Jamaica";}
   elsif ($city =~ /^(KITTY HAWK|KITTYHAWK)/i) {$lat=36.06; $lon=75.42; $locstring = "Kitty Hawk, N. Carolina";}
   elsif ($city =~ /^(KONA)/i) {$lat=20.38; $lon=156.00; $locstring="Kona, Hawaii";}
   elsif ($city =~ /^(KWAJALEIN)/i) {$lat=8.44; $lon=-167.44; $locstring="Kwajalein, Marshall Islands";}
   elsif ($city =~ /^(LAFAYETTE)/i) {$lat=30.12; $lon=92.01; $locstring = "Lafayette, Louisiana";}
   elsif ($city =~ /^(LAKE CHARLES|LAKECHARLES)/i) {$lat=30.13; $lon=93.13; $locstring = "Lake Charles, Louisiana";}
   elsif ($city =~ /^(LAPAZ|LA PAZ)/i) {$lat=24.16; $lon=110.25; $locstring = "La Paz, Mexico";}
   elsif ($city =~ /^(LEARMONTH)/i) {$lat=-22.14; $lon=-114.05; $locstring = "Learmonth, Australia";}
   elsif ($city =~ /^(LE PORT|LEPORT)/i) {$lat=-20.56; $lon=-55.17; $locstring = "Le Port, Reunion";}
   elsif ($city =~ /^(LIHUE)/i) {$lat=21.59; $lon=159.21; $locstring = "Lihue, Kauai";}
   elsif ($city =~ /^(LONG BEACH|LONGBEACH)/i) {$lat=33.49; $lon=118.09; $locstring = "Long Beach, California";}
   elsif ($city =~ /^(MANILA)/i) {$lat=14.35; $lon=-120.59; $locstring = "Manila, Philippines";}
   elsif ($city =~ /^(MANZANILLO)/i) {$lat=19.03; $lon=104.20; $locstring = "Manzanillo, Mexico";}
   elsif ($city =~ /^(MAZATLAN)/i) {$lat=23.12; $lon=106.25; $locstring = "Mazatlan, Mexico";}
   elsif ($city =~ /^(MELBOURNE)/i) {$lat=28.04; $lon=80.36; $locstring = "Melbourne, Florida";}
   elsif ($city =~ /^(MERIDA)/i) {$lat=20.59; $lon=89.39; $locstring = "Merida, Mexico";}
   elsif ($city =~ /^(MIA|MIAMI)/i) {$lat=25.45; $lon=80.23; $locstring = "Miami, Florida";}
   elsif ($city =~ /^(MIYAKO)/i) {$lat=24.47; $lon=-125.18; $locstring = "Miyako, Japan";}
   elsif ($city =~ /^(MOBILE)/i) {$lat=30.40; $lon=88.05; $locstring = "Mobile, Alabama";}
   elsif ($city =~ /^(MOLOKAI)/i) {$lat=21.09; $lon=157.06; $locstring = "Molokai";}
   elsif ($city =~ /^(MONTERREY)/i) {$lat=25.15; $lon=99.40; $locstring = "Monterrey, Mexico";}
   elsif ($city =~ /^(MOREHEAD CITY|MOREHEADCITY)/i) {$lat=33.48; $lon=77.00; $locstring = "Morehead City, N. Carolina";}
   elsif ($city =~ /^(MYRTLE BEACH|MYRTLEBEACH)/i) {$lat=33.41; $lon=78.56; $locstring = "Myrtle Beach, S. Carolina";}
   elsif ($city =~ /^(NAHA)/i) {$lat=26.11; $lon=-127.39; $locstring = "Naha, Okinawa";}
   elsif ($city =~ /^(NANTUCKET)/i) {$lat=41.17; $lon=70.05; $locstring = "Nantucket, Massachusetts";}
   elsif ($city =~ /^(NASSAU)/i) {$lat=25.03; $lon=77.28; $locstring = "Nassau, Bahamas";}
   elsif ($city =~ /^(NEW ORLEANS|NEWORLEANS)/i) {$lat=30.00; $lon=90.03; $locstring = "New Orleans, Louisiana";}
   elsif ($city =~ /^(NEW YORK JFK|NEWYORKJFK|NEWYORK|NEW YORK)/i) {$lat=40.39; $lon=73.47; $locstring = "New York (JFK), New York";}
   elsif ($city =~ /^(NGULU)/i) {$lat=8.18; $lon=-137.29; $locstring = "Ngulu, Caroline Islands";}
   elsif ($city =~ /^(NORFOLK)/i) {$lat=36.54; $lon=76.12; $locstring = "Norfolk, Virginia";}
   elsif ($city =~ /^(OCEANCITY|OCEAN CITY)/i) {$lat=38.20; $lon=75.05; $locstring = "Ocean City, Maryland";}
   elsif ($city =~ /^(PANAMA CITY|PANAMACITY)/i) {$lat=30.10; $lon=85.41; $locstring = "Panama City, Florida";}
   elsif ($city =~ /^(PENSACOLA)/i) {$lat=30.26; $lon=87.12; $locstring = "Pensacola, Florida";}
   elsif ($city =~ /^(PERTH)/i) {$lat=-31.56; $lon=-115.57; $locstring = "Perth, Australia";}
   elsif ($city =~ /^(PLAISANCE)/i) {$lat=-20.26; $lon=-57.40; $locstring = "Plaisance, Mauritius";}
   elsif ($city =~ /^(PORT HEDLAND|PORTHEDLAND)/i) {$lat=-20.19; $lon=-118.35; $locstring = "Port Hedland, Australia";}
   elsif ($city =~ /^(PORT LOUIS|PORTLOUIS)/i) {$lat=-20; $lon=-57; $locstring = "Port Louis, Mauritius";}
   elsif ($city =~ /^(PUERTOVALLARTA|PUERTO VALLARTA)/i) {$lat=20.40; $lon=105.15; $locstring = "Puerto Vallarta, Mexico";}
   elsif ($city =~ /^(SAINT DENIS|SAINTDENIS)/i) {$lat=-20.53; $lon=-55.31; $locstring = "Saint Denis, Reunion";}
   elsif ($city =~ /^(SAINT PIERRE|SAINTPIERRE)/i) {$lat=-21.20; $lon=-55.29; $locstring = "Saint Pierre, Reunion";}
   elsif ($city =~ /^(SAIPAN)/i) {$lat=15.10; $lon=-145.40; $locstring = "Saipan, Northern Mariana Is.";}
   elsif ($city =~ /^(SALINA CRUZ|SALINACRUZ)/i) {$lat=16.10; $lon=95.12; $locstring = "Salina Cruz, Mexico";}
   elsif ($city =~ /^(SAN DIEGO|SANDIEGO)/i) {$lat=32.51; $lon=117.07; $locstring = "San Diego, California";}
   elsif ($city =~ /^(SAN JUAN|SANJUAN)/i) {$lat=18.29; $lon=66.08; $locstring = "San Juan, Puerto Rico";}
   elsif ($city =~ /^(SANTODOMINGO|SANTO DOMINGO)/i) {$lat=18.26; $lon=69.40; $locstring = "Santo Domingo, Dominican Republic";}
   elsif ($city =~ /^(SAVANNAH)/i) {$lat=32.08; $lon=81.12; $locstring = "Savannah, Georgia";}
   elsif ($city =~ /^(SPNT|SOUTH POINT)/i) {$lat=19.00; $lon=155.40; $locstring="South Point, Hawaii";}
   elsif ($city =~ /^(SUBIC BAY|SUBICBAY)/i) {$lat=14.48; $lon=-120.16; $locstring="Subic Bay, Philippines";}
   elsif ($city =~ /^(STPETERSBURG|ST PETERSBURG|ST. PETERSBURG)/i) {$lat=27.45; $lon=82.44; $locstring = "St. Petersburg, Florida";}
   elsif ($city =~ /^(STTHOMAS|ST THOMAS|ST. THOMAS)/i) {$lat=18.20; $lon=64.58; $locstring = "St. Thomas, USVI";}
   elsif ($city =~ /^(TAINAN)/i) {$lat=23.00; $lon=-120.13; $locstring = "Tainan, Taiwan";}
   elsif ($city =~ /^(TAIPEI)/i) {$lat=25.20; $lon=-121.31; $locstring = "Taipei, Taiwan";}
   elsif ($city =~ /^(TALLAHASSEE)/i) {$lat=30.10; $lon=85.41; $locstring = "Tallahassee, Florida";}
   elsif ($city =~ /^(TAMPICO)/i) {$lat=22.18; $lon=97.52; $locstring = "Tampico, Mexico";}
   elsif ($city =~ /^(TAPACHULA)/i) {$lat=14.55; $lon=92.16; $locstring = "Tapachula, Mexico";}
   elsif ($city =~ /^(TINIAN)/i) {$lat=14.55; $lon=-145.30; $locstring = "Tinian, Northern Mariana Is.";}
   elsif ($city =~ /^(TOWNSVILLE)/i) {$lat=-19.15; $lon=-146.45; $locstring = "Townsville, Australia";}
   elsif ($city =~ /^(TOKYO)/i) {$lat=35.33; $lon=-139.47; $locstring = "Tokyo, Japan";}
   elsif ($city =~ /^(TRUK)/i) {$lat=7.28; $lon=-151.51; $locstring = "Truk, Caroline Islands";}
   elsif ($city =~ /^(UPTON)/i) {$lat=40.52; $lon=72.52; $locstring = "Upton, New York";}
   elsif ($city =~ /^(VACOAS)/i) {$lat=-20.18; $lon=-57.30; $locstring = "Vacoas, Mauritius";}
   elsif ($city =~ /^(VEROBEACH)/i) {$lat=27.39; $lon=80.25; $locstring = "Vero Beach, Florida";}
   elsif ($city =~ /^(VIRGINIA BEACH|VIRGINIABEACH)/i) {$lat=36.50; $lon=75.59; $locstring = "Virginia Beach, Virginia";}
   elsif ($city =~ /^(WILMINGTON)/i) {$lat=34.10; $lon=77.50; $locstring = "Wilmington, N. Carolina";}
   elsif ($city =~ /^(WELLINGTON)/i) {$lat=-41.20; $lon=-174.48; $locstring = "Wellington, New Zealand";}
   elsif ($city =~ /^(YAP)/i) {$lat=9.29; $lon=-138.05; $locstring = "Yap, Caroline Islands";}

   ($lat,$lon,$locstring);

}

sub GETLATLON {

   # extract the latitude, longitude, and date\time from an input line.

   my($line) = @_;
   my($lat,$lon,$tim);

   my(@line) = split(/ /,$line);

   $tim = $line[$DATE];
   $lat = $line[$LATITUDE];
   $lon = $line[$LONGITUDE];
   if ($line[$NS] eq 'S') {$lat = -$lat;}
   if ($line[$EW] eq 'E') {$lon = -$lon;}
   ($lat,$lon,$tim);

}

$FASTCUMM = 0;  # If set use the faster approximate CUMMULATE
                # But it's not that much faster and it's really
                # approximate.  So don't set this.

sub CUMMULATE {

   # Compute cumulative probability given that the storm was at the
   # requested location ($lat,$lon) at indices in @given.

   # Input: $yindex,$s,*given,*llarchive
   #    $yindex = The index in @llarchive of the starting forecast
   #    $s = The index into $STORMSIZE
   #    @given = The list of indices into @llarchive of the forecasts
   #             for which a strike is true
   #    @llarchive = list of all forecasts processed
   #
   # Output: Probability 
 
   my($yindex,$s,$rgiven,$rllarchive) = @_;

   my($cumoverlap,$loverlap,$toverlap,$ltoverlap,@g1,@g2,@given,@llarchive);

   @given = @$rgiven;
   @llarchive = @$rllarchive;

   $toverlap = &SPROBGIVEN($s,$yindex,\@given,\@llarchive);
   if ($yindex == 0) {
      $cumoverlap = $toverlap;
   } else {

      $loverlap  = &CUMMULATE($yindex-1,$s,\@given,\@llarchive);

      if ($FASTCUMM) {
         @g1 = ($yindex);
         $ltoverlap = &CUMMULATE($yindex-1,$s,\@g1,\@llarchive)
                    + $loverlap*(1 - $toverlap);
         if ($ltoverlap > 1) {$ltoverlap = 1;}
      } else {
         @g1 = ($yindex,@given);
         $ltoverlap = &CUMMULATE($yindex-1,$s,\@g1,\@llarchive);
      }

      my ($upper,$lower);
      if ($toverlap > 0) {
         $upper = $loverlap/$toverlap;
         $lower = ($toverlap+$loverlap-1)/$toverlap;
      } else {
         $upper = 1;
         $lower = 0;
      }
      if ($ltoverlap>$upper) {$ltoverlap=$upper;}
      if ($ltoverlap<$lower) {$ltoverlap=$lower;}
      #if ($ltoverlap<$lower) {$ltoverlap=$loverlap;}
      #if ($ltoverlap>$upper || $ltoverlap<$lower) {$ltoverlap=$loverlap;}

      $cumoverlap = $loverlap + $toverlap - $toverlap*$ltoverlap;

   }
   if ($cumoverlap >1) {$cumoverlap=1;}
   if ($cumoverlap <0) {$cumoverlap=0;}
   $cumoverlap;

}  # END OF CUMMULATE SUBROUTINE


sub SPROBGIVEN {

   # Compute the probability of a strike given a list of strikes
   # that are known to be true.

   # Input: $s,$yindex,*given,*llarchive
   #    $s = The index into $STORMSIZE
   #    $yindex = The index in @llarchive of the starting forecast
   #    @given = The list of indices into @llarchive of the forecasts
   #             for which a strike is true
   #    @llarchive = list of all forecasts processed

   my($s,$yindex,$rgiven,$rllarchive) = @_;
   my($strike,$toverlap,$tdate,$tlat,$tlon,$tstrike,$avgstrike,
         $ydate,$ylat,$ylon,$ystrike,@tstrike,@ystrike,@g1,@g2,$ytdiff,$ttdiff,$tdiff,
         @given,@llarchive);
   @given = @$rgiven;
   @llarchive = @$rllarchive;

   if ($#given <= 0) {
      ($ydate,$ylat,$ylon,$ytdiff,@ystrike) = split(" +",$llarchive[$yindex]);
      ($tdate,$tlat,$tlon,$ttdiff,@tstrike) = split(" +",$llarchive[$given[0]]);
      $tstrike = $tstrike[$s];
      $ystrike = $ystrike[$s];
      my($d) = &GCDISTANCE($ylat,$ylon,$tlat,$tlon);
      if ($yindex < $given[0]) {
          $toverlap = &AVGSPROB($STORMSIZE[$s],$ttdiff,$d,($tdate-$ydate)/$onehour,$tstrike,$ystrike);
          ##$toverlap = &AVGSPROB($STORMSIZE[$s],$ttdiff,$d,$ttdiff,$tstrike,$ystrike);
          if ($tstrike>0) {
             $toverlap *= ($ystrike/$tstrike);  # Bayes' Theorem
          } 
      } else {
          $toverlap = &AVGSPROB($STORMSIZE[$s],$ytdiff,$d,($ydate-$tdate)/$onehour,$ystrike,$tstrike);
          ##$toverlap = &AVGSPROB($STORMSIZE[$s],$ytdiff,$d,$ytdiff,$ystrike,$tstrike);
      }
      $strike = $toverlap;
   } else {
      @g1 = $given[$#given];
      @g2 = $given[0..($#given-1)];

      # These three versions give exactly the same result ...

      #$strike = &SPROBGIVEN($s,$yindex,\@g1,\@llarchive) 
      #        + &SPROBGIVEN($s,$yindex,\@g2,\@llarchive)
      #        * (1 - &SPROBGIVEN($s,$given[$#given],\@g2,\@llarchive));

      $strike =  &SPROBGIVEN($s,$yindex,\@g2,\@llarchive)
              + &SPROBGIVEN($s,$yindex,\@g1,\@llarchive) 
              * (1 - &SPROBGIVEN($s,$given[$#given],\@g2,\@llarchive));

      #@g1 = $given[0];
      #@g2 = $given[1..($#given)];
      #$strike =  &SPROBGIVEN($s,$yindex,\@g2,\@llarchive)
      #        + &SPROBGIVEN($s,$yindex,\@g1,\@llarchive) 
      #        * (1 - &SPROBGIVEN($s,$given[0],\@g2,\@llarchive));

   }

   if ($strike > 1) {$strike = 1;}
   if ($strike < 0) {$strike = 0;}

   $strike;

}

1