#!/usr/bin/perl -w

# pfit: linear least square fit to quadratic function
# self contained, needs only perl. Should run on most any
# unix or linux or Mac OS X machine.
# John Lapeyre Tue Jan 28 18:45:19 MST 2003
# lapeyre@physics.arizona.edu

# Usage: pfit < datafile
# read data from file, two columns separated by white space.
# Ignore comments beginning with '#'. 
# suitable for data files that are not too large, maybe 1000 pairs or
# so ?

# fit to A x^2  + B x + C
# Works by solving matrix equation M.c = v. 
# c = (A,B,C)
# v = (v1,v2,v3), where vn = \sum_i y_i x_i^(n-1)
# m = ( (m4,m3,m2), (m3,m2,m1), (m2,m1,m0))
# where mn = \sum_i x_i^n

@x=@y=(); # data arrays
$i=0;
while(<>) {   # read lines from standard input
  next  if /\#/;   # skip comments
  (@nums) = split;
  if (@nums == 1) {   # one column
      $x = $i;        # not efficient to check for nums every time
      $y = $nums[0];
  }
  elsif (@nums == 2) {  # two columns
      $x = $nums[0];
      $y = $nums[1];
  }
  else { die "Need one or two white space separated columns."}
  push @x,$x;
  push @y,$y;
#  print "$x $y\n";
  $i++;
}


$N= $#x + 1;  # number of data pairs

$v1=$v2=$v3=$m1=$m2=$m3=$m4=0;

$m0 = $N;

for($i=0;$i<$N;$i++) {
    $x=$x[$i];
    $y=$y[$i];
    $v1 += $y*$x**2;
    $v2 += $y*$x;
    $v3 += $y;
    $m1 += $x;
    $m2 += $x*$x;
    $m3 += $x**3;
    $m4 += $x**4;
}

# used mathematica to invert the matrix and translated to perl

$A = ($m1**2*$v1 - $m0*$m2*$v1 + $m0*$m3*$v2 + $m2**2*$v3 - $m1*($m2*$v2 + $m3*$v3))/
  ($m2**3 + $m0*$m3**2 + $m1**2*$m4 - $m2*(2*$m1*$m3 + $m0*$m4));
$B = (-($m1*$m2*$v1) + $m0*$m3*$v1 + $m2**2*$v2 - $m0*$m4*$v2 - $m2*$m3*$v3 + $m1*$m4*$v3)/
  ($m2**3 + $m0*$m3**2 + $m1**2*$m4 - $m2*(2*$m1*$m3 + $m0*$m4));
$C =  ($m2**2*$v1 - $m1*$m3*$v1 + $m1*$m4*$v2 + $m3**2*$v3 - $m2*($m3*$v2 + $m4*$v3))/
    ($m2**3 + $m0*$m3**2 + $m1**2*$m4 - $m2*(2*$m1*$m3 + $m0*$m4));

# A x^2  + B x + C
printf("%e\t%e\t%e\n",$A,$B,$C);