#!/usr/local/bin/perl -w

# Pair correlation analysis for two data files(the latter refers the formaer).
# Self correlation analysis for single data file.
# Correlation is limited for those closest, prior day to the reference.
# Data should be numerical and one column.
# Tue Jul 10 2001 ;  by Isoji MIYAGI


# data file reading. 
if(@ARGV==2){
    open(RAIN,"<$ARGV[0]") or die "cannot open rain data\n";
    @rains=<RAIN>;
    close(RAIN);
    open(ERUPT,"<$ARGV[1]") or die "cannot open eruption data\n";
    @eruptions=<ERUPT>;
    close(ERUPT);
    chomp (@eruptions, @rains);
    @eruptions = reverse(sort {$a <=> $b} @eruptions);
    @rains = reverse(sort {$a <=> $b} @rains);
    &catdog();
    &printGraph("pair correlation", "p", "");
} elsif (@ARGV==1){
    open(FILE,"<$ARGV[0]") or die "cannot open data\n";
    @myselves=<FILE>;
    close(FILE);
    chomp @myselves;
    @myselves = reverse(sort {$a <=> $b} @myselves);
    &selfish();
    &printGraph("self correlation", "s", "");
} else {
    print STDERR "data file open failed\n";
    exit 0;
}


# a subroutine for pair correlation 
sub catdog(){
    open(OUT,">./deleteMe.txt");
  eruptionDay:
    foreach $theEruptionday (@eruptions) {
      rainDay:
        foreach $theRainday (@rains) {
            $delDay = sprintf("%d", $theRainday - $theEruptionday);
            if($delDay <= 0) {
                $N{$delDay}++;
                print STDOUT "$theRainday $theEruptionday $delDay\n";
                next eruptionDay;
            } else {
                next rainDay;
            }
        }
    }
    foreach $theDelDay (sort{$a <=> $b}keys(%N)) {
        print OUT "$theDelDay\t$N{$theDelDay}\n";
    }
    close(OUT);
}


# a subroutine for self correlation 
sub selfish () {
    open(OUT,">./deleteMe.txt");
  rainReferenceDay:
    foreach $theDay (@myselves) {
      rainDay:
        foreach $theOtherDay (@myselves) {
            $delDay = sprintf("%d", $theOtherDay - $theDay);
            if($delDay < 0) {
                print STDOUT "$theOtherDay $theDay $delDay\n";
                $N{$delDay}++;
                $N{0}++;
                next rainReferenceDay;
            } else {
                next rainDay;
            }
        }
    }
    foreach $theDelDay (sort{$a <=> $b}keys(%N)) {
        print OUT "$theDelDay\t$N{$theDelDay}\n";
    }
    close(OUT);
}


# a subroutine for graph making
sub printGraph(){
    my($title, $type, @xrange)=@_;
    foreach $theXRange (@xrange) {
    open (GNUP, "|/usr/bin/gnuplot") || die "could not open GNUPLOT\n";
    print GNUP <<PERIOD;
set terminal png small color
set title \"$title\"
set output \"$ARGV[0]-$ARGV[1]$type$theXRange.png\"
set xlabel \"rainy day after the eruption\"
set ylabel \"number\"
plot  \[$theXRange\]\[0:\] \'deleteMe.txt\' title  \"$ARGV[0]-$ARGV[1]\"  with impulses 

PERIOD
close (GNUP);
    }
}
exit 0;