#!/usr/bin/perl

# Interpolate positions for all the new RN obs from the
#  subset with positions.

# Uses linear interpolation using 2 constraints
#  1) Don't interpolate over more than 48 hours
#  2) Don't interpolate over more than 5 degrees (lat or long)
# If adjacent good positions fail these checks, intermediate obs will
#  not have a position interpolated.

use strict;
use warnings;
use IMMA;
use Time::Local;
use Date::Calc qw(Delta_DHMS);    # LocalTime no good for dates before 1904

# Read in the obs and sort into date order
my @a;
while ( my $ob = imma_read( \*STDIN ) ) {
    push @a, $ob;
}
@a = sort td_compare @a;

# Mark the obs with positions
my @WithP;
foreach (@a) {
    if ( defined( $_->{LAT} ) && defined( $_->{LON} ) ) {
        push @WithP, 1;
    }
    else {
        push @WithP, 0;
    }
}

# Interpolate the missing locations
for ( my $i = 0 ; $i < scalar(@a) ; $i++ ) {
    if ( defined( $a[$i]->{LAT} ) || defined( $a[$i]->{LON} ) ) {
        next;
    }
    my $Previous = find_previous($i);
    my $Next     = find_next($i);
    if (   defined($Previous)
        && defined($Next) )
    {
        my ( $Lat1, $Lon1 ) = interpolate( $Previous, $Next, $a[$i] );
        if (   defined($Lat1)
            && defined($Lon1) )
        {
            $a[$i]->{LAT} = $Lat1;
            $a[$i]->{LON} = $Lon1;
            $a[$i]->{LI}  = 3;       # Mark position as interpolated
        }
    }

}

# Output the new obs
foreach (@a) { $_->write( \*STDOUT ); }

# Find the last previous ob that has a date
sub find_previous {
    my $Point = shift;
    for ( my $j = $Point - 1 ; $j >= 0 ; $j-- ) {
        if ( $a[$j]->{ID} ne $a[$Point]->{ID} ) { return; }
        if ( $WithP[$j] == 1 ) { return $a[$j]; }
    }
    return;
}

# Find the first subsequent ob that has a date
sub find_next {
    my $Point = shift;
    for ( my $j = $Point + 1 ; $j < scalar(@a) ; $j++ ) {
        if ( $a[$j]->{ID} ne $a[$Point]->{ID} ) { return; }
        if ( $WithP[$j] == 1 ) { return $a[$j]; }
    }
    return;
}

# Interpolate from the positions of a previous and a subsequent ob
sub interpolate {
    my $Previous = shift;
    my $Next     = shift;
    my $Target   = shift;

    # Check that all obs have valid dates
    unless ( checkValid($Previous) && checkValid($Next) && checkValid($Target) )
    {
        return ( undef, undef );
    }

    # Give up if the gap is longer than 10 days hours
    if ( Delta_Seconds( $Previous, $Next ) > 10 * 86400 ) {
        return ( undef, undef );
    }

    # Deal with any logitude wrap-arounds
    my $Next_long = $Next->{LON};
    if ( $Next_long - $Previous->{LON} > 180 ) { $Next_long -= 360; }
    if ( $Next_long - $Previous->{LON} < -180 ) { $Next_long += 360; }

    # Give up if the separation is too great
    if (   abs( $Next_long - $Previous->{LON} ) > 20
        || abs( $Next->{LAT} - $Previous->{LAT} ) > 10 )
    {
        return ( undef, undef );
    }

    # Do the interpolation
    if ( Delta_Seconds( $Target, $Next ) < 0 ) { return ( undef, undef ); }
    my $Weight =
      Delta_Seconds( $Target, $Next ) / Delta_Seconds( $Previous, $Next );
    if ( $Weight < 0 || $Weight > 1 ) { return ( undef, undef ); }
    my $Target_long = $Next_long * ( 1 - $Weight ) + $Previous->{LON} * $Weight;
    if ( $Target_long < -180 ) { $Target_long += 360; }
    if ( $Target_long > 180 ) { $Target_long -= 360; }
    my $Target_lat =
      $Next->{LAT} * ( 1 - $Weight ) + $Previous->{LAT} * $Weight;
    return ( $Target_lat, $Target_long );
}

# Compare two records, sort by callsign then date
sub td_compare {
    my ( $Idb, $Yrb, $Mob, $Dyb, $Hrb ) =
      ( $a->{ID}, $a->{YR}, $a->{MO}, $a->{DY}, $a->{HR} );
    unless ( defined($Idb) ) { $Idb = 'SHIP     '; }
    unless ( defined($Yrb) ) { $Yrb = 3000; }
    unless ( defined($Mob) ) { $Mob = 12; }
    unless ( defined($Dyb) ) { $Dyb = 31; }
    unless ( defined($Hrb) ) { $Hrb = 24; }
    my ( $Ida, $Yra, $Moa, $Dya, $Hra ) =
      ( $a->{ID}, $a->{YR}, $a->{MO}, $a->{DY}, $a->{HR} );
    unless ( defined($Ida) ) { $Ida = 'SHIP     '; }
    unless ( defined($Yra) ) { $Yra = 3000; }
    unless ( defined($Moa) ) { $Moa = 12; }
    unless ( defined($Dya) ) { $Dya = 31; }
    unless ( defined($Hra) ) { $Hra = 24; }
    return $Ida cmp $Idb ||    # Compare callsign
      $Yra <=> $Yrb ||         # Compare Year
      $Moa <=> $Mob ||         # Compare Month
      $Dya <=> $Dyb ||         # Compare Day
      $Hra <=> $Hrb;           # Compare Hour
}

sub month_lengths {
    my $Month = shift;
    my $Year  = shift;
    unless ( defined($Month) && defined($Year) ) { return; }
    my @Lengths = qw(31 28 31 30 31 30 31 31 30 31 30 31);
    if ( $Year % 400 == 0 || ( $Year % 4 == 0 && $Year % 100 != 0 ) ) {
        $Lengths[1] = 29;
    }
    return $Lengths[ $Month - 1 ];
}

# Difference between 2 dates in seconds
sub Delta_Seconds {
    my $First = shift;
    my $Last  = shift;
    my ( $Dd, $Dh, $Dm, $Ds ) = Delta_DHMS(
        $First->{YR},
        $First->{MO},
        $First->{DY},
        int( $First->{HR} ),
        int( ( $First->{HR} - int( $First->{HR} ) ) * 60 ),
        0,
        $Last->{YR},
        $Last->{MO},
        $Last->{DY},
        int( $Last->{HR} ),
        int( ( $Last->{HR} - int( $Last->{HR} ) ) * 60 ),
        0
    );
    return $Dd * 86400 + $Dh * 3600 + $Dm * 60 + $Ds;
}

# Check an ob has a valid data
sub checkValid {
    my $Ob = shift;
    if (   !defined( $Ob->{HR} )
        || $Ob->{HR} < 0
        || $Ob->{HR} > 24
        || !defined( $Ob->{DY} )
        || $Ob->{DY} < 1
        || !defined( $Ob->{MO} )
        || $Ob->{MO} < 1
        || $Ob->{MO} > 12
        || $Ob->{DY} > month_lengths( $Ob->{MO}, $Ob->{YR} ) )
    {
        return;
    }
    else {
        return 1;
    }
}