#!/usr/bin/perl
#
# Compute chi-squared for linear and quadratic models 
#    to some data.
#
# MWR 4/25/2012

use POSIX;

$debug = 1;

# read in data
$infile = "./electron.dat";
$ndata = 0;
open(INFILE, "$infile") || die("Can't open file $infile");
while (<INFILE>) {
  $line = $_;
  if ($line =~ /^#/) {
    next;
  }
  @words = split(/\s+/, " " . $line);
  $time_array[$ndata] = $words[1];
  $pos_array[$ndata] = $words[2];
  $posunc_array[$ndata] = $words[3];
  if ($debug > 0) {
    printf " data %3d   %5.2f  %5.2f  %5.2f \n", $ndata, 
           $time_array[$ndata],
           $pos_array[$ndata],
           $posunc_array[$ndata];
  }
  $ndata++;
}
close(INFILE);


# calc chisq for linear fit
$a = 1.54;
$b = 9.69;

$chisq = 0;
for ($i = 0; $i < $ndata; $i++) {
  $model_pos = $a + $b*$time_array[$i];
  $diff = $model_pos - $pos_array[$i];
  $top = $diff*$diff;
  $bot = $posunc_array[$i]*$posunc_array[$i];
  $chisq += $top / $bot;
  if ($debug > 0) {
    printf " linear:  data %4d  top %8.2f bot %8.2f  this %8.2f  sum %8.2f \n",
           $i, $top, $bot, $top/$bot, $chisq;
  }
}

$red_chisq = $chisq / ($ndata - 2);
printf " chisq linear %8.2f  %8.2f \n", $chisq, $red_chisq;


# calc chisq for quadratic fit
$a = 0.26;
$b = 11.60;
$c = -0.382;

$chisq = 0;
for ($i = 0; $i < $ndata; $i++) {
  $model_pos = $a + $b*$time_array[$i] + $c*$time_array[$i]*$time_array[$i];
  $diff = $model_pos - $pos_array[$i];
  $top = $diff*$diff;
  $bot = $posunc_array[$i]*$posunc_array[$i];
  $chisq += $top / $bot;
  if ($debug > 0) {
    printf " quad:    data %4d  top %8.2f bot %8.2f  this %8.2f  sum %8.2f \n",
           $i, $top, $bot, $top/$bot, $chisq;
  }
}
$red_chisq = $chisq / ($ndata - 3);
printf " chisq quad   %8.2f  %8.2f \n", $chisq, $red_chisq;







exit 0;


