#!/usr/local/bin/perl
# Pair correlation analysis for two data files(the latter refers the formaer).
# Self correlation analysis for single data file.
# Data should be numerical and two column separated by \t.
# Mon Jun 25 2001 ;  by Isoji MIYAGI
# Thu Jun 28 2001 ; modified by Isoji MIYAGI
# Tue Jul 10 2001 ; modified by Isoji MIYAGI

@plotRanges = ("42", "365");
@thresholds = ("0", "10", "50");

# data file reading. 
if(@ARGV==2){
    open(RAIN,"<$ARGV[0]") or die "cannot open rain data\n";
    @otherselves=<RAIN>;
    close(RAIN);
    open(ERUPT,"<$ARGV[1]") or die "cannot open eruption data\n";
    @myselves=<ERUPT>;
    close(ERUPT);
    chomp (@myselves, @otherselves);
    &catdog();
} elsif (@ARGV==1){
    $ARGV[1]=$ARGV[0];
    open(FILE,"<$ARGV[0]") or die "cannot open data\n";
    @myselves=<FILE>;
    close(FILE);
    chomp @myselves;
    &selfish();
} else {
    print STDERR "data file open failed\n";
    exit 0;
}


# a subroutine for pair correlation 
sub catdog(){
    my($eruptData, $rainData, $theDay);
    foreach $rainData (@otherselves) {
        ($theDay, $fall) = split("\t", $rainData);
        next if $rainData =~ /\#|X/;
        $rainFall{$theDay} = $fall;
    }
    foreach $eruptData (@myselves) {
        ($theDay, $magnitude) = split("\t", $eruptData);
        next if $eruptData =~ /\#|X/;
        $eruptMag{$theDay} = $magnitude;
    }
    
    foreach $theThreshold (@thresholds) {
        %N=0;
        open(OUT,">./deleteMe.txt");
        foreach $theRainyday (keys(%rainFall)) {
            next unless ($rainFall{$theRainyday} > $theThreshold);
            print STDOUT "$rainFall{$theRainyday}\t$theThreshold\n";
            print OUT "#\t$theRainyday";
            foreach $theEruptionday (keys(%eruptMag)) {
                $delDay = sprintf("%d", $theRainyday - $theEruptionday);
                print OUT "\t$delDay";
                $N{$delDay}++;
            }
            print OUT "\n";
        }
        foreach $theDelDay (sort{$a <=> $b}keys(%N)) {
            print OUT "$theDelDay\t$N{$theDelDay}\n";
        }
        close(OUT);
        
        &printGraph($theThreshold, @plotRanges);
    }
}


# a subroutine for self correlation 
sub selfish(){
    my($rainData, $theDay);
    foreach $rainData (@myselves) {
        ($theDay, $fall) = split("\t", $rainData);
        next if $rainData =~ /\#|X/;
        $rainFall{$theDay} = $fall;
    }
    foreach $theThreshold (@thresholds) {
        %N=0;
        open(OUT,">./deleteMe.txt");
        foreach $theRainyday (keys(%rainFall)) {
            next unless ($rainFall{$theRainyday} > $theThreshold);
#            print STDOUT "$rainFall{$theRainyday}\t$theThreshold\n";
            print OUT "#\t$theRainyday";
            foreach $theOtherRainyday (keys(%rainFall)) {
                next unless ($rainFall{$theOtherRainyday} > $theThreshold);
                $delDay = sprintf("%d", $theRainyday - $theOtherRainyday);
                print OUT "\t$delDay";
                $N{$delDay}++;
            }
            print OUT "\n";
        }
        foreach $theDelDay (sort{$a <=> $b}keys(%N)) {
            print OUT "$theDelDay\t$N{$theDelDay}\n";
        }
        close(OUT);
        
        &printGraph($theThreshold, @plotRanges);
    }
}


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