#!/usr/bin/perl
#
# Compute a rough approximation to the motion (in RA and Dec)
#   of the star Ross 248 over a period of several years.
#
# MWR 5/28/2020
#
# Fixed two mistakes in calculation delta_dec;
#   don't use Nakanishi et al. (2015); instead, rely upon
#   Smart and Taff books.
#
# MWR 6/8/2020
#
# Modified to compute the motion of GX And.
#
# MWR 6/14/2020
#
# Compare measurements of the position of GX And to a model
#   for its motion.  Find the best-fit values of a number
#   of parameters.
#
# Usage:    perl fit_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/15/2020



use POSIX;

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;


###########################################################
#    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 = 180.0;
# proper motion in Dec (milliarcsec per year)
#     This is Gaia value
#$mu_dec = 411.9;
$mu_dec = 220.0;
# parallax (milliarcsec)
#     This is Gaia value
#$parallax = 280.7;
$parallax = 90.0;
#
# cosine of declination of star
$cosdec = cos($DEGTORAD*$dec_start);
# RA zeropoint offset (mas)
$ra_offset = -10.0;
# Dec zeropoint offset (mas)
$dec_offset = 30.0;


###############################################################
#       Free Parameters for the Fit
#
# these will the offsets of the star from its starting position
#    in milliarcsec
$delta_ra = 0.0;
$delta_dec = 0.0;


###########################################################
#       Initial ranges for fitting parameters
#
# We'll use a grid search to find the best-fit values
#   for each parameter.  These are the initial sizes
#   of the range on each side of the starting value
#   that we'll check.  The ranges will decrease
#   after each iteration.
#
# For example, if  
#                        parallax       = 300
#                        parallax_range = 100
#
# then we'll check values of 200 < parallax < 400.
#
$parallax_range = 50.0;    
$mu_ra_range = 80.0;
$mu_dec_range = 80.0;
$ra_offset_range = 60.0;
$dec_offset_range = 60.0;
#
# after each iteratio, shrink each range by this factor
$SHRINK_FACTOR = 0.85;


##########################################################
#    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];
  $num_pos++;
}
close(POSITION);
if ($debug > 1) {
  printf " read %d entries in file %s \n", $num_pos, $pos_data_file;
}



################################################################
#     Set up for iterating to find best values

$num_param = 5;
$param_val[0] = $parallax;
$param_range[0] = $parallax_range;
$param_val[1] = $mu_ra;
$param_range[1] = $mu_ra_range;
$param_val[2] = $mu_dec;
$param_range[2] = $mu_dec_range;
$param_val[3] = $ra_offset;
$param_range[3] = $ra_offset_range;
$param_val[4] = $dec_offset;
$param_range[4] = $dec_offset_range;

# number of iterations to perform
$num_iter = 10;
# break up range into this many pieces during each iteration
$num_piece = 21;



for ($i_iter = 0; $i_iter < $num_iter; $i_iter++) {

  if ($debug > 1) {
    printf " top of loop for iter %d \n", $i_iter;
    for ($i_param = 0; $i_param < $num_param; $i_param++) {
      printf "    param %2d  val %8.1f  range %7.1f \n", 
             $i_param, $param_val[$i_param], $param_range[$i_param];
    }
  }

  for ($i_param = 0; $i_param < $num_param; $i_param++) {
    if ($debug > 1) {
      printf "    checking param %d \n", $i_param;
    }

    # start by assigning the current values to all temp param variables
    for ($i = 0; $i < $num_param; $i++) {
      $test_param[$i] = $param_val[$i];
    }
    # now we will loop over the range for param $i_param
    $min = $test_param[$i_param] - $param_range[$i_param];
    $max = $test_param[$i_param] + $param_range[$i_param];
    $dx = ($max - $min) / $num_piece;
    if ($debug > 1) {
      printf "     param %2d  loop from %7.1f to %7.1f  by %7.1f \n",
               $i_param, $min, $max, $dx;
    }

    $test_param[$i_param] = $min;
    for ($i_piece = 0; $i_piece < $num_piece; $i_piece++) {

        $chisq = compute_chisq($test_param[0], $test_param[1],
                          $test_param[2], $test_param[3], $test_param[4]);
        if ($debug > -1) { 
          printf "# chisq %8.2f  for param %2d = %8.1f  \n", 
                   $chisq, $i_param, $test_param[$i_param];
        }

        if ($i_piece == 0) {
          $min_chisq = $chisq;
          $best_val = $test_param[$i_param];
          if ($debug > 1) {
            printf "#      initial min_chisq  %10.1f     for val %6.1f \n", 
                           $min_chisq, $best_val;
          }
        }
        else {
          if ($chisq < $min_chisq) {
            $min_chisq = $chisq;
            $best_val = $test_param[$i_param];
            if ($debug > 1) {
              printf "#      new     min_chisq  %10.1f   for val %6.1f \n", 
                           $min_chisq, $best_val;
            }
          }
        }

        $test_param[$i_param] += $dx;

    }

    # save the best value for use in next iteration
    $new_val[$i_param] = $best_val;

  }  # end of loop over all parameters

  # now, we replace the value of each parameter with the new best value
  for ($i_param = 0; $i_param < $num_param; $i_param++) {
    if ($debug > 1) {
      printf "#    replacing param %2d  %7.1f  with %7.1f \n",
             $i_param, $param_val[$i_param], $new_val[$i_param];
    }
    $param_val[$i_param] = $new_val[$i_param];
  }

  # finally, we shrink the range of each parameter
  for ($i_param = 0; $i_param < $num_param; $i_param++) {
    $new_range = $param_range[$i_param]*$SHRINK_FACTOR;
    if ($debug > 1) {
      printf "#    replacing param %2d  range  %7.1f  with %7.1f \n",
             $i_param, $param_range[$i_param], $new_range;
    }
    $param_range[$i_param] = $new_range;

  }

} # end of loop over all iterations


# print out a summary of all parameter values,
#    as well as a list of all data and the best-fit model
printf "# end of fitting after %d iterations \n", $num_iter;
$chisq = compute_chisq($param_val[0], $param_val[1],
                          $param_val[2], $param_val[3], $param_val[4]);
printf "# N = %6d data, chisq is %7.1f \n", $num_pos, $chisq;
for ($i_param = 0; $i_param < $num_param; $i_param++) {
  printf "#    param %2d  val %8.1f  range %7.1f \n", 
         $i_param, $param_val[$i_param], $param_range[$i_param];
  # special case -- add an extra line for the params
  #          parallax   pm_ra   pm_dec
  if ($i_param == 0) { 
    printf "#     solution   parallax   %9.1f \n", $param_val[$i_param];
  }
  if ($i_param == 1) { 
    printf "#     solution   pm_ra      %9.1f \n", $param_val[$i_param];
  }
  if ($i_param == 2) { 
    printf "#     solution   pm_dec     %9.1f \n", $param_val[$i_param];
  }
}
compare_model_and_data($param_val[0], $param_val[1],
         $param_val[2], $param_val[3], $param_val[4]);


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)/365.2422;

  # 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];

  }

}
