#!/usr/bin/perl

# chi2.pl
# Vypocet chi^2 pro rozdil merenych bodu od teoretickeho ciselniku.
# Calculate chi^2 for the difference between the measured points and theoretical dial.
# Miroslav Broz (miroslav.broz@email.cz), May 15th 2006

$Sigma = 0.5;	# cm, stredni chyba mereni

########################################################################

# read points to optimize

open(IN,"<BODY.dat_optimalizace");
read(IN,$in,-s IN);
close(IN);

$i = 0;
foreach $line (split("\n",$in)) {
  if ($line !~ /^#.*/) {
    @L = split(" ",$line);
    if ($#L > 1) {
      $i++;
      $bod[$i]    = $L[0];
      $x[$i]      = $L[1];
      $y[$i]      = $L[2];
      $file[$i]   = $L[3];
      $param1[$i] = $L[4];
      $param2[$i] = $L[5];

      $sigma_x[$i] = $Sigma;
      $sigma_y[$i] = $Sigma;
    }
  }
}
$n = $i;
print "# n_bod = $n\n";

########################################################################

# read date curves

open(IN,"<datum_prumet.dat");
read(IN,$in,-s IN);
close(IN);

$i = 0;
foreach $line (split("\n",$in)) {
  @L = split(" ",$line);
  if ($line =~ /# l = /) {
    $l = $L[3];
  } elsif (($line !~ /^#.*/) && ($#L > 1)) {
    $i++;
    $t_datum[$i] = $L[0];
    $x_datum[$i] = $L[3];
    $y_datum[$i] = $L[4];
    $l_datum[$i] = $l;
  }
}
$n_datum = $i;
print "# n_datum = $n_datum\n";

########################################################################

# read hour lines

open(IN,"<hodin_rysky.dat");
read(IN,$in,-s IN);
close(IN);

$i = 0;
$eps = 1.e-3;
foreach $line (split("\n",$in)) {
  @L = split(" ",$line);
  if (($line !~ /^#.*/) && ($#L > 1)) {
    $i++;
    $t_hodin[$i] = $L[0];
    $x_hodin[$i] = $L[1];
    $y_hodin[$i] = $L[2];
    if (abs($t_hodin[$i] - $t_hodin_last) < $eps) {	# letni/zimni hodinova znacka?
      $l_hodin[$i] = 2;
    } else {
      $l_hodin[$i] = 1;
    }
    $t_hodin_last = $t_hodin[$i];
  }
}
$n_hodin = $i;
print "# n_hodin = $n_hodin\n";

########################################################################

# calculate the chi^2 value

$chi2 = 0;
$eps = 1.e-3;

print "# BOD & dx [cm] & dy [cm] & chi^2\n";

for ($i = 1; $i <= $n; $i++) {

# match date curves

  if ($file[$i] eq "datum_prumet.dat") {
    $j = 1;
    while (($j <= $n_datum) and ((abs($param1[$i]-$l_datum[$j]) > $eps) || (abs($param2[$i]-$t_datum[$j]) > $eps))) {
#      print "$j $param1[$i] $l_datum[$j] $param2[$i] $t_datum[$j]\n";	# dbg
      $j++;
    }
    if ($j <= $n_datum) {
      $dx[$i] = $x[$i]-$x_datum[$j];
      $dy[$i] = $y[$i]-$y_datum[$j];
      $d[$i] = ($dx[$i]**2 + $dy[$i]**2) / $Sigma;
      $chi2 = $chi2 + $d[$i];
      printf("# bod %s %10.2f %10.2f %12.4f\n", $bod[$i], $dx[$i], $dy[$i], $d[$i]);

    } else {
      print "# Chyba: Bod `$bod[$i]' nenalezen v souboru `$file[$i]'.\n";
    }

####################################

# match hour lines

  } elsif ($file[$i] eq "hodin_rysky.dat") {
    $j = 1;
    while (($j <= $n_hodin) and ((abs($param1[$i]-$t_hodin[$j]) > $eps) || (abs($param2[$i]-$l_hodin[$j]) > $eps))) {
#      print "$j $param1[$i] $t_hodin[$j] $param2[$i] $l_hodin[$j]\n";	# dbg
      $j++;
    }
    if ($j <= $n_hodin) {
      $dx[$i] = $x[$i]-$x_hodin[$j];
      $dy[$i] = $y[$i]-$y_hodin[$j];
      $d[$i] = ($dx[$i]**2 + $dy[$i]**2) / ($sigma_x[$i]**2 + $sigma_y[$i]**2);
      $chi2 = $chi2 + $d[$i];
      printf("# bod %s %10.2f %10.2f %12.4f\n", $bod[$i], $dx[$i], $dy[$i], $d[$i]);

    } else {
      print "# Chyba: Bod `$bod[$i]' nenalezen v souboru `$file[$i]'.\n";
    }

####################################

  } else {
    print "# Chyba: Soubor `$file[$i]', odkazovany u bodu `$bod[$i]', neni nacten -> upravit skript chi2.pl?.\n";
  }
}

########################################################################

# print the results

print "sigma = $Sigma cm\n";
print "N = $n\n";
printf("chi2 = %.2f\n", $chi2);
print "\n";

exit;