#!/usr/bin/perl
#
# Compare measurements of the position of GX And to a model
#   for its motion.  Find the best-fit values of a number
#   of parameters.
#
# This script uses a least-squares matrix inversion approach 
#   to determine the best-fit values.
#
# Usage:    perl matrix_gxand.pl  positions.dat
#
#   where
#       positions.dat     is an ASCII text file with positions
#                            of the star on a series of dates.
#                            Format is
#                                      
#              JD   ra_pos ra_uncert   dec_pos dec_uncert
#
#                         where the "ra_pos" and "dec_pos" are
#                         offsets from a reference position (in mas),
#                         and "ra_uncert" and "dec_uncert" are
#                         uncertainties in those offsets (in mas)
#
# MWR 6/18/2020
#
# This version uses the uncertainties in each measurement of position
#   to give weights in the least-squares calculation.
#
# MWR 6/25/2021
#
# Found and fixed an error in one term of the matrix calculation, 
#   for element matrix_A(3, 3)
#
# MWR 6/29/2021



use POSIX;
use Math::MatrixReal;

my $debug = 0;
my $mypi = 3.14159265;
my $DEGTORAD = ($mypi/180.0);

##########################################################
#    Time variables
#
# the base date for calculations is Jan 0, 2000.0
$jd_base = 2451543.5;
# this corresponds to May 26, 2020 = JD 2458995.0,
#    which is roughly the start date for the observations
$start_date = 2458995.5;
# this is several years later than the start 
#    of observations
$end_date = $start_date + (365.25*3.0);
# time step (in days)
$d_date = 10.0;
# length of year, in days
$year_length = 365.2422;


###########################################################
#    Parameters for the Earth's orbit
#
# obliquity of Earth's orbit (degrees)
$obliq = 23.440972;
$cos_obliq = cos($obliq*$DEGTORAD);
$sin_obliq = sin($obliq*$DEGTORAD);


############################################################
#    Properties of the star's motion
#
# RA of star at starting date (decimal degrees)
$ra_start = 4.618213;
$ra_rad = $ra_start*$DEGTORAD;
# Dec of star at starting date (decimal degrees)
$dec_start = 44.025266;
$dec_rad = $dec_start*$DEGTORAD;
# proper motion in RA (milliarcsec per year)
#     This is Gaia value
$mu_ra = 2891.5;
#$mu_ra = 2893.0;
# proper motion in Dec (milliarcsec per year)
#     This is Gaia value
$mu_dec = 411.9;
#$mu_dec = 411.0;
# parallax (milliarcsec)
#     This is Gaia value
$parallax = 280.7;
#$parallax = 258.2;
#
# cosine of declination of star
$cosdec = cos($DEGTORAD*$dec_start);



###############################################################
#        Arbitrary offset to position calculated by model
#
#   Add these amounts in RA and Dec to every position calculated
#   by our theoretical model.
#
# RA zeropoint offset (mas)
$ra_offset = -29.4;
#$ra_offset = 0.0;
# Dec zeropoint offset (mas)
$dec_offset = 74.4;
#$dec_offset = 0.0;


###############################################################
#       Pre-computed Modifications to the input measurements
#
# these are the offsets of the star from its starting position
#    in milliarcsec.  We have taken their values from the 
#    results of the grid-method approach.
#    Subtract these values from the measurements as we read them in.
$delta_ra = 0.0;
$delta_dec = 0.0;



##########################################################
#    Parse args
#
if ($#ARGV != 0) {
  printf "usage: fit_gxand.pl  positions.dat  \n";
  exit(1);
}
my $pos_data_file = $ARGV[0];
if ($debug > 1) {
  printf "# pos_data_file is %s \n", $pos_data_file;
}


##########################################################
#    Read data from the given datafile
#
$num_pos = 0;
open(POSITION, $pos_data_file) || die("can't open file $pos_data_file");
while (<POSITION>) {
  $line = $_;
  if ($line =~ /^#/) {
    next;
  }
  @words = split(/\s+/, $line);
  if ($#words != 4) {
    printf " invalid input line has %d words, not expected 5 \n", $#words;
    printf "%s", $line;
    exit(1);
  }
  $jd_array[$num_pos] = $words[0];
  $ra_pos_array[$num_pos] = $words[1];
  $ra_uncert_array[$num_pos] = $words[2];
  $dec_pos_array[$num_pos] = $words[3];
  $dec_uncert_array[$num_pos] = $words[4];

  # subtract our guessed zero-point offsets from the RA and Dec
  #   positions
  $ra_pos_array[$num_pos] -= $delta_ra;
  $dec_pos_array[$num_pos] -= $delta_dec;
  $num_pos++;
}
close(POSITION);
if ($debug > 1) {
  printf " read %d entries in file %s \n", $num_pos, $pos_data_file;
}



############################################################
#     Compute quantities in matrix A (on left-hand side)
#
$matrix_A = new Math::MatrixReal(3, 3);

for ($i = 0; $i < $num_pos; $i++) {
  $date = $jd_array[$i];
  $lambda_rad = compute_lambda($date);
  $ty = ($date - $start_date) / $year_length;
  if ($debug > 1) {
    printf "  date %12.5f  ty %6.3f   lambda_rad %9.4f \n",
           $date, $ty, $lambda_rad;
  }

  my $p = -sin($ra_rad)*cos($lambda_rad) + 
                   $cos_obliq*cos($ra_rad)*sin($lambda_rad);
  my $q = cos($dec_rad)*$sin_obliq*sin($lambda_rad)
                   - cos($ra_rad)*sin($dec_rad)*cos($lambda_rad)
                   - sin($ra_rad)*sin($dec_rad)*$cos_obliq*sin($lambda_rad);

  # compute a weight corresponding to the measurement in each direction
  my $weight_ra = 1.0 / $ra_uncert_array[$i];
  my $weight_dec = 1.0 / $dec_uncert_array[$i];
  # and since we'll always be using the squares of each weight, 
  #   let's compute those now
  my $wrasq = $weight_ra * $weight_ra;
  my $wdecsq = $weight_dec * $weight_dec;

  $matrix_A->assign(1, 1, $matrix_A->element(1, 1) + 
                           ($p*$p*$wrasq + $q*$q*$wdecsq) );
  $matrix_A->assign(1, 2, $matrix_A->element(1, 2) + ($p*$ty*$wrasq) );
  $matrix_A->assign(1, 3, $matrix_A->element(1, 3) + ($q*$ty*$wdecsq) );

  $matrix_A->assign(2, 1, $matrix_A->element(2, 1) + ($p*$ty*$wrasq) );
  $matrix_A->assign(2, 2, $matrix_A->element(2, 2) + ($ty*$ty*$wrasq) );
  $matrix_A->assign(2, 3, 0.0);

  $matrix_A->assign(3, 1, $matrix_A->element(3, 1) + ($q*$ty*$wdecsq) );
  $matrix_A->assign(3, 2, 0.0);
  $matrix_A->assign(3, 3, $matrix_A->element(3, 3) + ($ty*$ty*$wdecsq) );
}
if ($debug > 0) {
  printf "   here comes 3x3 matrix A \n";
  printf "   %9.4f %9.4f %9.4f \n", 
               $matrix_A->element(1, 1), $matrix_A->element(1, 2), 
                    $matrix_A->element(1, 3);
  printf "   %9.4f %9.4f %9.4f \n", 
               $matrix_A->element(2, 1), $matrix_A->element(2, 2), 
                    $matrix_A->element(2, 3);
  printf "   %9.4f %9.4f %9.4f \n", 
               $matrix_A->element(3, 1), $matrix_A->element(3, 2), 
                    $matrix_A->element(3, 3);
}





############################################################
#     Compute quantities in vector b (on right-hand side)
#
$vector_b = new Math::MatrixReal(1, 3);

for ($i = 0; $i < $num_pos; $i++) {
  $date = $jd_array[$i];
  $lambda_rad = compute_lambda($date);
  $ty = ($date - $start_date) / $year_length;
  if ($debug > 1) {
    printf "  date %12.5f  ty %6.3f   lambda_rad %9.4f \n",
           $date, $ty, $lambda_rad;
  }

  # compute a weight corresponding to the measurement in each direction
  my $weight_ra = 1.0 / $ra_uncert_array[$i];
  my $weight_dec = 1.0 / $dec_uncert_array[$i];
  # and since we'll always be using the squares of each weight, 
  #   let's compute those now
  my $wrasq = $weight_ra * $weight_ra;
  my $wdecsq = $weight_dec * $weight_dec;

  $p = -sin($ra_rad)*cos($lambda_rad) + 
                   $cos_obliq*cos($ra_rad)*sin($lambda_rad);
  $q = cos($dec_rad)*$sin_obliq*sin($lambda_rad)
                   - cos($ra_rad)*sin($dec_rad)*cos($lambda_rad)
                   - sin($ra_rad)*sin($dec_rad)*$cos_obliq*sin($lambda_rad);
  
  # I may need to add or subtract the "$ra_offset" and "$dec_offset"
  #      values to/from the positions before using them 
  #      in the following calculations ...
  my $ra_fix = $ra_pos_array[$i];
  my $dec_fix = $dec_pos_array[$i];

  $vector_b->assign(1, 1, $vector_b->element(1, 1) + 
                          ($p*$ra_fix*$wrasq + $q*$dec_fix*$wdecsq)  );
  $vector_b->assign(1, 2, $vector_b->element(1, 2) +  
                          ($ra_fix*$ty*$wrasq)  );
  $vector_b->assign(1, 3, $vector_b->element(1, 3) +  
                          ($dec_fix*$ty*$wdecsq) );

}
if ($debug > 0) {
  printf "   here comes 1x3 vector b \n";
  printf "   %9.4f %9.4f %9.4f \n", 
           $vector_b->element(1, 1), $vector_b->element(1, 2), 
                         $vector_b->element(1, 3);
}




############################################################
#     Invert maxtrix A
$inverse_A = $matrix_A->inverse();
#if ($inverse_A == undef) {
#  printf " can't invert matrix_A ... quitting \n";
#  exit(1);
#}
if ($debug > 0) {
  printf " succeeded in inverting matrix_A --- here comes inverse \n";
  printf "   %9.4f %9.4f %9.4f \n", 
               $inverse_A->element(1, 1), $inverse_A->element(1, 2), 
                    $inverse_A->element(1, 3);
  printf "   %9.4f %9.4f %9.4f \n", 
               $inverse_A->element(2, 1), $inverse_A->element(2, 2), 
                    $inverse_A->element(2, 3);
  printf "   %9.4f %9.4f %9.4f \n", 
               $inverse_A->element(3, 1), $inverse_A->element(3, 2), 
                    $inverse_A->element(3, 3);
}





############################################################
#     Compute values of vector x, which contains the desired
#         parameters we are fitting
#
$vector_x = $vector_b->multiply($inverse_A);
if ($debug > 0) {
  printf "   here comes 1x3 solution vector x \n";
  printf "   %9.4f %9.4f %9.4f \n", 
           $vector_x->element(1, 1), $vector_x->element(1, 2), 
                         $vector_x->element(1, 3);
}


##################################################################
#      Compute the chisq value of the model
$chisq = compute_chisq($vector_x->element(1, 1), $vector_x->element(1, 2),
                       $vector_x->element(1, 3), $ra_offset, $dec_offset);
$red_chisq = $chisq / ($num_pos - 3.0);
printf "#    chisq %9.2f    red_chisq %9.2f \n", $chisq, $red_chisq;





# print out a summary of all parameter values,
#    as well as a list of all data and the best-fit model
printf "#    solution   parallax  %8.1f  \n", 
         $vector_x->element(1, 1);
printf "#    solution   pm_ra     %8.1f  \n", 
         $vector_x->element(1, 2);
printf "#    solution   pm_dec    %8.1f  \n", 
         $vector_x->element(1, 3);

  compare_model_and_data($vector_x->element(1, 1), $vector_x->element(1, 2),
         $vector_x->element(1, 3), $ra_offset, $dec_offset);

# we can do this stuff if desired
if (0 == 1) {

  calc_cov_matrix($matrix_A);

  $this_parallax = $vector_x->element(1, 1);
  $this_mu_ra = $vector_x->element(1, 2);
  $this_mu_dec = $vector_x->element(1, 3);
  my $s_squared = compute_s_squared($this_parallax, $this_mu_ra, $this_mu_dec);
}


exit 0;


############################################################################
# PROCEDURE: degrees_in_range
#
# DESCRIPTION: Given an angle in degrees, add or subtract 360.0
#              as needed to bring it into the range 0 - 360 degrees.
#
# RETURNS:   the converted angle
#
sub degrees_in_range {
  my $angle = $_[0];

  while ($angle < 0.0) {
    $angle += 360.0;
  }
  while ($angle >= 360.0) {
    $angle -= 360.0;
  }

  return($angle);
}



############################################################
# PROCEDURE: compute_star_ra_dec   
#
# DESCRIPTION: Given a date (Julian Date),
#              compute the expected position of the star.
#              Calculate the position in milliarcsec
#              away from the starting position.
#
#              Input parameters are
#
#                   date            Julian Date 
#                   parallax        parallax, in mas
#                   mu_ra           proper motion RA, mas/year
#                   mu_dec          proper motion Dec, mas/year
#                   ra_offset       RA zeropoint offset, mas
#                   dec_offset      Dec zeropoint offset, mas
#
#              We use a bunch of hard-coded values 
#              in this routine to compute the solar
#              longitude and ecliptic longitude;
#              in addition, we refer to global variables
#              defined earlier.
#
# RETURNS:     (ra_mas, dec_mas)
#
#                     ra_mas        offset in RA from the starting
#                                      position (mas)
#
#                     dec_mas       offset in Dec frm the starting
#                                      position (mas)
#

sub compute_star_ra_dec {
  my ($jd, $day_offset, $year_offset);
  my ($parallax, $mu_ra, $mu_dec);
  my ($ra_offset, $dec_offset);
  my ($d, $L, $g, $lambda, $lambda_rad);
  my ($x1, $x2, $x3, $delta_ra, $delta_dec);

  $jd = $_[0];
  $parallax = $_[1];
  $mu_ra = $_[2];
  $mu_dec = $_[3];
  $ra_offset = $_[4];
  $dec_offset = $_[5];

  # days since $start_date
  $day_offset = $jd - $start_date;
  # years since $start_date
  $year_offset = ($day_offset)/$year_length;

  # compute the solar ecliptic longitude for this date
  #   days since Jan 0, 2000
  $d = $jd - $jd_base;
  #   mean longitude
  $L = 280.461 + 0.9856474*$d;
  $L = degrees_in_range($L);
  #   mean anomaly 
  $g = 357.528 + 0.9856003*$d;
  $g = degrees_in_range($g);
  #   ecliptic longitude
  $lambda = $L + 1.915*sin($g*$DEGTORAD) + 0.020*sin(2.0*$g*$DEGTORAD);
  $lambda = degrees_in_range($lambda);
  $lambda_rad = $lambda*$DEGTORAD;

  if ($debug > 1) {
    printf " Q: day %8d  L %8.2f  g %8.2f  lambda %8.2f ", $d, $L, $g, $lambda;
  }

  # now compute the offset of the star from its starting point 
  #    in RA (milliarcseconds)
  $x1 = $parallax*(-sin($ra_rad)*cos($lambda_rad) + 
                   $cos_obliq*cos($ra_rad)*sin($lambda_rad));
  $x2 = $mu_ra*$year_offset;
  $delta_ra = $x1 + $x2;
  $delta_ra += $ra_offset;

  if ($debug > 1) {
    printf " x1 %8.2f x2 %8.2f  delta_ra %8.2f ", $x1, $x2, $delta_ra;
  }

  # now compute the offset of the star from its starting point 
  #    in Dec (milliarcseconds)
  $x1 = $parallax*(cos($dec_rad)*$sin_obliq*sin($lambda_rad)
                   - cos($ra_rad)*sin($dec_rad)*cos($lambda_rad)
                   - sin($ra_rad)*sin($dec_rad)*$cos_obliq*sin($lambda_rad) );
  $x3 = $mu_dec*$year_offset;
  $delta_dec = $x1 + $x3;
  $delta_dec += $dec_offset;

  if ($debug > 1) {
    printf " x1 %8.2f x3 %8.2f delta_dec %8.2f ", 
            $x1, $x3, $delta_dec;
    printf "\n";
  }


  # return the position: offset in milli-arcsec (delta_ra, delta_dec)
  return($delta_ra, $delta_dec);
}





####################################################
# PROCEDURE: compute_chisq
#
# DESCRIPTION: Compute the value of chisq for a given
#              set of parameters.  
#
# USAGE:     compute_chisq  parallax mu_ra mu_dec ra_offset dec_offset
#
# RETURNS:   The value of chisq
#
# 
sub compute_chisq {

  my($parallax, $mu_ra, $mu_dec, $ra_offset, $dec_offset);
  my($chisq);
  my($t_ra, $t_dec);

  $parallax = $_[0];
  $mu_ra = $_[1];
  $mu_dec = $_[2];
  $ra_offset = $_[3];
  $dec_offset = $_[4];



  #########################################################
  #     Test main loop: print measurements and predictions
  #

  $chisq = 0.0;

  for ($i = 0; $i < $num_pos; $i++) {

    my $date = $jd_array[$i];

	 # compute theoretical position
	 ($t_ra, $t_dec) = compute_star_ra_dec($date, $parallax,
								 $mu_ra, $mu_dec, $ra_offset, $dec_offset);
	 if ($debug > 1) {
		printf "  on date %12.5f   t_ra %7.1f  t_dec %7.1f \n",
				 $date, $t_ra, $t_dec;
	 }

	 # print out the prediction, then the measurement
    if ($debug > 1) {
      printf "%12.5f  %7.1f %7.1f   %7.1f %6.1f  %7.1f %6.1f \n",
		  			$date, $t_ra, $t_dec, 
					$ra_pos_array[$i], $ra_uncert_array[$i],
					$dec_pos_array[$i], $dec_uncert_array[$i];
    }


    my $dra = ($t_ra - $ra_pos_array[$i]) / $ra_uncert_array[$i];
    my $ddec = ($t_dec - $dec_pos_array[$i]) / $dec_uncert_array[$i];
    my $diff_sq = $dra*$dra + $ddec*$ddec;
    $chisq += $diff_sq;

  }

  if ($debug > 1) { 
    printf "# chisq %8.2f \n", $chisq;
  }


  return($chisq);
}


###########################################################
# PROCEDURE: compare_model_and_data
#
# DESCRIPTION: For each measured position, compute the 
#              corresponding theoretical position, and print
#              out one line per date.
#
#              Arguments are:
#
#                  parallax mu_ra mu_dec ra_offset dec_offset
#
# RETURNS:     nothing
#
#
sub compare_model_and_data {

  my($parallax, $mu_ra, $mu_dec, $ra_offset, $dec_offset);
  my($date, $t_ra, $t_dec);

  $parallax = $_[0];
  $mu_ra = $_[1];
  $mu_dec = $_[2];
  $ra_offset = $_[3];
  $dec_offset = $_[4];


  for ($i = 0; $i < $num_pos; $i++) {

    $date = $jd_array[$i];

	 # compute theoretical position
	 ($t_ra, $t_dec) = compute_star_ra_dec($date, $parallax,
								 $mu_ra, $mu_dec, $ra_offset, $dec_offset);
	 if ($debug > 1) {
		printf "  on date %12.5f   t_ra %7.1f  t_dec %7.1f \n",
				 $date, $t_ra, $t_dec;
	 }

	 # print out the prediction, then the measurement
      printf "%12.5f  %7.1f %7.1f   %7.1f %6.1f  %7.1f %6.1f \n",
		  			$date, $t_ra, $t_dec, 
					$ra_pos_array[$i], $ra_uncert_array[$i],
					$dec_pos_array[$i], $dec_uncert_array[$i];

  }

}


############################################################
# PROCEDURE: compute_lambda
#
# DESCRIPTION: Given a date (Julian Date),
#              compute the ecliptic longitude of the Sun
#              on that date, in radians.
#
#              Input parameters are
#
#                   date            Julian Date 
#
#              We use a bunch of hard-coded values 
#              in this routine to compute the solar
#              longitude and ecliptic longitude;
#              in addition, we refer to global variables
#              defined earlier.
#
# RETURNS:     lambda
#
#                     lambda        solar longitude (radians)

sub compute_lambda {
  my ($jd, $day_offset, $year_offset);
  my ($d, $L, $g, $lambda, $lambda_rad);
  my ($x1, $x2, $x3, $delta_ra, $delta_dec);

  $jd = $_[0];

  # days since $start_date
  $day_offset = $jd - $start_date;
  # years since $start_date
  $year_offset = ($day_offset)/$year_length;

  # compute the solar ecliptic longitude for this date
  #   days since Jan 0, 2000
  $d = $jd - $jd_base;
  #   mean longitude
  $L = 280.461 + 0.9856474*$d;
  $L = degrees_in_range($L);
  #   mean anomaly 
  $g = 357.528 + 0.9856003*$d;
  $g = degrees_in_range($g);
  #   ecliptic longitude
  $lambda = $L + 1.915*sin($g*$DEGTORAD) + 0.020*sin(2.0*$g*$DEGTORAD);
  $lambda = degrees_in_range($lambda);
  $lambda_rad = $lambda*$DEGTORAD;

  if ($debug > 1) {
    printf " Q: day %8d  L %8.2f  g %8.2f  lambda %8.2f ", $d, $L, $g, $lambda;
  }

  return($lambda_rad);
}




############################################################
# PROCEDURE: calc_cov_matrix
#
# DESCRIPTION: Given a matrix A, compute the covariance
#              matrix C.
#
# RETURNS:     the new matrix C
#

sub calc_cov_matrix {
  my ($matrix_A, $matrix_At, $matrix_H, $matrix_C);

  $matrix_A = $_[0];

  if ($debug > 0) {
    printf " calc_cov_matrix: here comes input matrix A \n";
    printf "   %12.6f %12.6f %12.6f \n", 
               $matrix_A->element(1, 1), $matrix_A->element(1, 2), 
                    $matrix_A->element(1, 3);
    printf "   %12.6f %12.6f %12.6f \n", 
               $matrix_A->element(2, 1), $matrix_A->element(2, 2), 
                    $matrix_A->element(2, 3);
    printf "   %12.6f %12.6f %12.6f \n", 
               $matrix_A->element(3, 1), $matrix_A->element(3, 2), 
                    $matrix_A->element(3, 3);
  }

  # first, we compute the transpose of A
  #    create a matrix of the same size as A ....
  $matrix_At = $matrix_A->shadow();
  #    and then fill it with the transpose of A
  $matrix_At->transpose($matrix_A);
  if ($debug > 0) {
    printf " calc_cov_matrix: here comes transposed matrix At \n";
    printf "   %12.6f %12.6f %12.6f \n", 
               $matrix_At->element(1, 1), $matrix_At->element(1, 2), 
                    $matrix_At->element(1, 3);
    printf "   %12.6f %12.6f %12.6f \n", 
               $matrix_At->element(2, 1), $matrix_At->element(2, 2), 
                    $matrix_At->element(2, 3);
    printf "   %12.6f %12.6f %12.6f \n", 
               $matrix_At->element(3, 1), $matrix_At->element(3, 2), 
                    $matrix_At->element(3, 3);
  }

  # the next step is to multiply the two transposed matrices 
  #    to create the intermediate matrix H
  $matrix_H = $matrix_At->multiply($matrix_A);
  if ($debug > 0) {
    printf " calc_cov_matrix: here comes temporary matrix H \n";
    printf "   %12.6f %12.6f %12.6f \n", 
               $matrix_H->element(1, 1), $matrix_H->element(1, 2), 
                    $matrix_H->element(1, 3);
    printf "   %12.6f %12.6f %12.6f \n", 
               $matrix_H->element(2, 1), $matrix_H->element(2, 2), 
                    $matrix_H->element(2, 3);
    printf "   %12.6f %12.6f %12.6f \n", 
               $matrix_H->element(3, 1), $matrix_H->element(3, 2), 
                    $matrix_H->element(3, 3);
  }

  # finally, by inverting the temporary matrix H, we create
  #   the desired covariance matrix C
  $matrix_C = $matrix_H->inverse;
  if ($debug > 0) {
    printf " calc_cov_matrix: here comes covariance matrix C \n";
    printf "   %12.6f %12.6f %12.6f \n", 
               $matrix_C->element(1, 1), $matrix_C->element(1, 2), 
                    $matrix_C->element(1, 3);
    printf "   %12.6f %12.6f %12.6f \n", 
               $matrix_C->element(2, 1), $matrix_C->element(2, 2), 
                    $matrix_C->element(2, 3);
    printf "   %12.6f %12.6f %12.6f \n", 
               $matrix_C->element(3, 1), $matrix_C->element(3, 2), 
                    $matrix_C->element(3, 3);
  }

  return($matrix_C);
}



############################################################
# PROCEDURE: compute_s_squared
#
# DESCRIPTION: Given some set of parameters -- 
# 
#                   parallax, mu_ra, mu_dec
#
#              cycle through all the dates of measurement.
#              For each date, calculate the theoretical position
#              of the target (RA_t, Dec_t), and compare to the
#              measured (RA_m, Dec_m).  Compute the sum of the squares
#              of the differences:
#
#                 s_squared  =  Sum [  (RA_t - RA_m)^2 + (Dec_t - Dec_m)^2 ]
#                              ----------------------------------------------
#                                              N - 3
#
#              where N is the number of measurements.  
#
# RETURNS:     s_squared
#

sub compute_s_squared {
  my ($parallax, $mu_ra, $mu_dec);
  my ($diff, $sum);
  my ($i, $date);
  my ($t_ra, $t_dec);
  my ($dra, $ddec, $diff_sq, $sum);

  $parallax = $_[0];
  $mu_ra = $_[1];
  $mu_dec = $_[2];

  if ($debug > 0) {
    printf " compute_s_squared: parallax  %9.2f \n", $parallax;
    printf " compute_s_squared: mu_ra     %9.2f \n", $mu_ra;
    printf " compute_s_squared: mu_dec    %9.2f \n", $mu_dec;
  }

  if ($num_pos < 4) {
    printf " compute_s_squared: invalid, num_pos = %d, must be >= 4 \n",
        $num_pos;
    exit(1);
  }


  $sum = 0;

  for ($i = 0; $i < $num_pos; $i++) {

    $date = $jd_array[$i];

	 # compute theoretical position
	 ($t_ra, $t_dec) = compute_star_ra_dec($date, $parallax,
								 $mu_ra, $mu_dec, $ra_offset, $dec_offset);
	 if ($debug > 1) {
		printf "  on date %12.5f   t_ra %7.1f  t_dec %7.1f \n",
				 $date, $t_ra, $t_dec;
	 }

    $dra = $t_ra - $ra_pos_array[$i];
    $ddec = $t_dec - $dec_pos_array[$i];

    $dra = ($t_ra - $ra_pos_array[$i]) / $ra_uncert_array[$i];
    $ddec = ($t_dec - $dec_pos_array[$i]) / $dec_uncert_array[$i];


    $diff_sq = ($dra*$dra + $ddec*$ddec);
    $sum += $diff_sq;

    if ($debug > 0) {
      printf "%12.5f  %7.1f %7.1f  %7.1f %7.1f  %6.1f %6.1f   %8.2f %8.2f   \n",
		  			$date, $t_ra, $t_dec, 
					$ra_pos_array[$i], $dec_pos_array[$i], 
               $dra, $ddec, $diff_sq, $sum;
    }

  }

  my $s_squared = $sum / ($num_pos - 3.0);
  if ($debug > 0) {
    printf " s_squared is %9.2f \n", $s_squared;
  }


  return($s_squared);
}

