#!/usr/bin/env perl
# relative_error : calculates the relative error defined as
# abs(y1-y2)/((abs(y1)+abs(y2)/2)). The y1, y2 values are
# taken from a chosen column in two files as a function
# of the (common) first column.
# Part of "NRG Ljubljana".
# Rok Zitko, rok.zitko@ijs.si, Sep 2009
# $Id$

# CHANGE LOG
# 18.9.2009 - initial version
# 10.3.2011 - verbosity

use warnings;
use strict;

my $verbose = 0;

my $col =   (@ARGV >= 1 ? $ARGV[0] : "2");
my $fn1 =   (@ARGV >= 2 ? $ARGV[1] : "file1");
my $fn2 =   (@ARGV >= 3 ? $ARGV[2] : "file2");
my $fnout = (@ARGV >= 4 ? $ARGV[3] : "diff$col.dat");

$verbose && print "Column: $col\n";

open (F1, "<$fn1") or die "Can't open $fn1: $!\n";
open (F2, "<$fn2") or die "Can't open $fn2: $!\n";

my $lastline = ""; # last header line
my $data = 0; # flag - have we encountered numeric data yet

my ($l1, $l2);

while (<F1>) {
    chomp;
    if (/^#/) {
        $lastline = $_;
    } else {
        $l1 = $_;
        last;
    }
}
$verbose && print "1: $lastline\n";

while (<F2>) {
    chomp;
    if (/^#/) {
        $lastline = $_;
    } else {
        $l2 = $_;
        last;
    }
}
$verbose && print "2: $lastline\n";

open(OUT, ">$fnout") or die "Can't open $fnout for writing: $!\n";

$col -= 1; # offset 1

while (1) {
    my @data1 = split(' ', $l1);
    my @data2 = split(' ', $l2);
    my $x = $data1[0];
    ($x eq $data2[0]) or die "Mismatch!\n";
    my $y1 = $data1[$col];
    my $y2 = $data2[$col];

    my $absdiff = abs($y1-$y2);
    my $avgabs = (abs($y1)+abs($y2))/2;
    my $ratio = $absdiff/$avgabs;
    
    my $s = "$x $ratio $y1 $y2\n";
    $verbose && print $s;
    print OUT $s;

    $l1 = <F1>;
    $l1 or last;
    chomp($l1);
    $l2 = <F2>;
    $l2 or die "Failed reading file 2.\n";
    chomp($l2);
}

close(OUT);

