#!/bin/perl

{
$eol="\n";
$npts=0;
$junk=<STDIN>;

print "\start slope intercept= $junk \n";

print "\nx,y pairs - raw data\n";
while(<STDIN>)
{
   chomp;
  ($x1,$y1)=split/[,\t\s]+/;
   $x[$npts]=$x1;
   $y[$npts]=$y1;
   $npts=$npts+1;
   print "$x1 , $y1\n";
}

  print "npts= $npts\n\n";

  ($a,$b,$siga,$sigb,$variance,$r2)=&ols();

  print "slope = $a \n";
  print "intercept = $b \n";
  print "sigslope =  $siga \n";
  print "sigint =  $sigb \n";
  print "r-squared = $r2 \n";
  print "variance = $variance \n";


  print $eol;
  ($a,$siga,$variance,$r2right,$r2wrong)=&rto();

  print "slope = $a \n";
  print "sigslope = $siga \n";
  print "r-squared = $r2right\n";
  print "r2wrong = $r2wrong\n";
  print "variance = $variance \n";
  print $eol;

}
####################################################################

sub rto
{
#
# least squares fit to straight line through origin  y = a*x
#
# initialize summing variables
#
my $xsum=0.0;
my $ysum=0.0;
my $xysum=0.0;
my $xxsum=0.0;
my $yysum=0.0;
#
my $i=0;
my $x1=0.0; my $y1=0.0;
#
# calculate the sums needed
#
for($i=0;$i<$npts;$i++)
{
  $x1=$x[$i];
  $y1=$y[$i];
  $xsum=$xsum+$x1;
  $ysum=$ysum+$y1;
  $xxsum=$xxsum+$x1*$x1;
  $xysum=$xysum+$x1*$y1;
  $yysum=$yysum+$y1*$y1;
}
#
# calculate slope
# and error statistics
#
my $a=$xysum/$xxsum;
my $ybar=$ysum/$npts;
my $bottomsum=$yysum-2.0*$ybar*$ysum+$npts*$ybar*$ybar; # sum((y-yavg)**2)
my $topsum=$yysum-2.0*$a*$xysum+$a*$a*$xxsum;      # sum(residual**2)
my $variance=($yysum-$a*$a*$xxsum)/($npts-1);
my $sigma=$variance**0.5;
my $vara=$variance/$xxsum;
my $siga=$vara**0.5;
my $r2wrong=1.0-($topsum/$bottomsum);
my $r2right=$a*$a*$xxsum/$yysum;
my @retval=($a,$siga,$variance,$r2right,$r2wrong);
}

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

sub ols
{
#
# least squares fit to straight line  y = a*x + b
#
# initialize summing variables
#
my $sum=0.0;
my $xsum=0.0;
my $ysum=0.0;
my $xysum=0.0;
my $xxsum=0.0;
my $yysum=0.0;
my $a=0.0;
my $b=0.0;
my $i=0;
my $x1=0.0; my $y1=0.0;
#
# calculate the sums needed
#
for($i=0;$i<$npts;$i++)
{
  $x1=$x[$i];
  $y1=$y[$i];
  $sum=$sum+1.0;
  $xsum=$xsum+$x1;
  $ysum=$ysum+$y1;
  $xxsum=$xxsum+$x1*$x1;
  $xysum=$xysum+$x1*$y1;
  $yysum=$yysum+$y1*$y1;
}
#
# calculate slope and intercept
# and error statistics
#
my $det=$npts*$xxsum-$xsum*$xsum;
$b=($xxsum*$ysum-$xsum*$xysum)/$det;
$a=($npts*$xysum-$xsum*$ysum)/$det;
my $variance=($yysum+$b*$b*$sum+$a*$a*$xxsum-2.0*($b*$ysum+$a*$xysum-$b*$a*$xsum))/($npts-2);
my $sigma=$variance**0.5;
my $vara=$variance/$det*$npts;
my $varb=$variance/$det*$xxsum;
my $siga=$vara**0.5;
my $sigb=$varb**0.5;
my $r2=($a*($xysum-$xsum*$ysum/$npts))/($yysum-$ysum*$ysum/$npts);
my @retval=($a,$b,$siga,$sigb,$variance,$r2);
}