use NetCDF;
require ("timelocal.pl");

$strikeNum = 0;
$minute = 60;
$PI = 3.14159265358979;
# approximate radius at 40N
$radiusEarth = 6365;

##################################################################
# get most recent filename
##################################################################
# path
$netcdfPath = "/data/fxa/point/binLightning/netcdf/";

# files in directory
opendir(DIR, $netcdfPath);
@fileList = reverse sort readdir(DIR);
closedir(DIR);

# most recent file
$netcdfFile = $fileList[0];
$netcdfFilename = $netcdfPath . $netcdfFile;
##################################################################


##################################################################
# get lat and lon from command line
##################################################################
if($ARGV[1]) { 
    $latLoc = $ARGV[0];
    $lonLoc = $ARGV[1];
} else {
    die "usage: $0 lat(decimal degrees) lon(decimal degrees)\n";
}

$latA = deg2rad($latLoc);
$lonA = deg2rad($lonLoc);

$distanceToRDA = 111.12;
$closestDistance = $distanceToRDA;
##################################################################


##################################################################
# open netcdf file
##################################################################
my $ncid = NetCDF::open($netcdfFilename, NetCDF::NOWRITE);

# get number of records
my $numRecords = NetCDF::dimid($ncid, "strikeNum");

# get dimensions using dimension inquiry
NetCDF::diminq($ncid, $numRecords, \$max, \$lastRecord);

##################################################################
# cdf variables
##################################################################

# read lat, lon, and time of strike each into a list
$varid = NetCDF::varid($ncid, "lat");
NetCDF::varget($ncid, $varid, 0, $lastRecord, \@lats);
    
$varid = NetCDF::varid($ncid, "lon");
NetCDF::varget($ncid, $varid, 0, $lastRecord, \@lons);

$varid = NetCDF::varid($ncid, "time");
NetCDF::varget($ncid, $varid, 0, $lastRecord, \@times);

##################################################################
# close netcdf file
##################################################################
NetCDF::close($ncid);


##################################################################
# iterate through list to find closest 
# strike to supplied lat and lon coords
##################################################################
for $index (0..$lastRecord) {

    $latB = deg2rad($lats[$index]);
    $lonB = deg2rad(abs($lons[$index]));

    $distance = $radiusEarth * arcCos(sin($latA) * sin($latB) + cos($latA) * cos($latB) * cos(abs($lonA - $lonB)));
    
    if($distance < $closestDistance) {
        $closestDistance = $distance;
        $mylat = $latB;
        $mylon = $lonB;
        $strikeTime = $times[$index];
    }
}


##################################################################
# print out results
##################################################################
if($closestDistance != $distanceToRDA) {
    # format time into human readable form
    @strikeTime2 = gmtime($strikeTime);
    ($sec, $min, $hour, $mday, $mon, $year, $wday, $ampm, $isdst) = @strikeTime2;

    # format time for  output
    $mon++;
    $year += 1900;
    $mon = substr("0" . $mon, -2);
    $hour = substr("0". $hour, -2);
    $min = substr("0" . $min, -2);
    $strikeTime = $hour . $min . "Z " . $mon . "/" . $mday . "/" . $year;

    printf "Lightning strike within 60NM of location at: %2.2fN: %3.2fW at %s\n", rad2deg($mylat), rad2deg($mylon), $strikeTime;
        
} else {
    print "No lightning detected within 60NM of location\n";
}
##################################################################



##################################################################
# subroutines begin here
##################################################################
sub deg2rad {
    my $degrees = shift;
    return ($degrees / 180) * $PI;
}

sub rad2deg {
    my $radians = shift;
    return ($radians / $PI) * 180;
}

sub arcCos {
    my $x = shift;
    return atan2(sqrt(1 - $x * $x), $x);
}