#!/usr/bin/perl
#
# Given a set of parameters
#
#          parallax            parallax in mas
#          mu_ra               proper motion in RA, mas
#          mu_dec              proper motion in Dec, mas
#          ra_offset           zero-point offset in RA, mas
#          dec_offset          zero-point offset in Dec, mas
#          pos_uncert          uncertainty in each measurement, mas
#
# generate a set of (RA, Dec) positions over a range
# of JD values.  
#
# Print to stdout a list of positions in the following format:
#
#        JD           RA      sigma_RA     Dec      sigma_Dec
#   2459092.79887    924.3      40.9      400.9      59.8 
#
# where all the RA and Dec values are in mas.
#
# MWR 6/29/2021
#
# Modified to show the motion of a star near the equatorial equator.
#
# MWR 8/11/2021


use POSIX;
use Math::MatrixReal;

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

# set the random number generator here.
#    pick one of the two options for the value of seed:
#         a fixed value, will always generate identical data
#      or
#         value based on current time, will be different each run
$seed = 0;
$seed = time();
srand($seed);


##########################################################
#    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
#  how many years to generate?
$num_year = 3.0;
$end_date = $start_date + (365.25*$num_year);
# time step (in days)
$d_date = 10.0;
#     Should we skip measurements during solar conjunction?
#           1 = yes, 0 = no
$skip_conjunction = 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
#
#    This position is near the Ecliptic Equator,
#       so we should see edge-on parallax motion.
$ra_start = 116.0;
$dec_start = 21.0;
#
#
#    This position is chosen so that the ecliptic longitude is 
#       the same as that of GX And
$ra_start = 127.0;
$dec_start = 58.0;
#
#
$ra_start = 127.0;
$dec_start = 58.0;
# RA of star at starting date (decimal degrees)
$ra_rad = $ra_start*$DEGTORAD;
# Dec of star at starting date (decimal degrees)
$dec_rad = $dec_start*$DEGTORAD;
# cosine of declination of star
$cosdec = cos($dec_rad);



##########################################################
#    Parse args
#
if ($#ARGV != 5) {
  printf "usage: generate_ra_dec_faceon.pl  parallax mu_ra mu_dec ra_offset dec_offset  pos_uncert \n";
  exit(1);
}
my $parallax = $ARGV[0];
my $mu_ra = $ARGV[1];
my $mu_dec = $ARGV[2];
my $ra_offset = $ARGV[3];
my $dec_offset = $ARGV[4];
my $pos_uncert = $ARGV[5];
if ($debug > 1) {
  printf "# parallax   is %6.1f \n", $parallax;
  printf "# mu_ra      is %6.1f \n", $mu_ra;
  printf "# mu_dec     is %6.1f \n", $mu_dec;
  printf "# ra_offset  is %6.1f \n", $ra_offset;
  printf "# dec_offset is %6.1f \n", $dec_offset;
  printf "# pos_uncert is %6.1f \n", $pos_uncert;
}




# print out a summary of all parameter values,
#    as well as a list of all data and the best-fit model
printf "#    input      parallax   %8.1f  \n", $parallax;
printf "#    input      pm_ra      %8.1f  \n", $mu_ra;
printf "#    input      pm_dec     %8.1f  \n", $mu_dec;
printf "#    input      ra_zp      %8.1f  \n", $ra_offset;
printf "#    input      dec_zp     %8.1f  \n", $dec_offset;
printf "#    input      pos_uncert %8.1f  \n", $pos_uncert;


#####################################################################
#        Compute positions, and print out as we go
#
for ($jd = $start_date; $jd <= $end_date; $jd += $d_date) {

    # check to see if we should skip this date (because star isn't
    #     visible in the night sky)
    if ($skip_conjunction == 1) {
      my $year_since_start = ($jd - $start_date) / $year_length;
      my $frac_year = $year_since_start - int($year_since_start);
      if (($frac_year > 0.7) && ($frac_year < 0.9)) {
        next;
      }
      else {
      }
    }
      

  ($ra, $dec) = compute_star_ra_dec($jd, $parallax, $mu_ra, $mu_dec,
                         $ra_offset, $dec_offset);
  
  # generate random uncertainties in position
  $ra_uncert = $pos_uncert*gaussian_rand();
  $dec_uncert = $pos_uncert*gaussian_rand();

  $measured_ra = $ra + $ra_uncert;
  $measured_dec = $dec + $dec_uncert;

  printf "%12.5f  %8.2f %7.2f  %8.2f %7.2f  \n", 
            $jd, $measured_ra, $pos_uncert, $measured_dec, $pos_uncert;

}


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);
}






######################################################################
# PROCEDURE: gaussian_rand
#
# DESCRIPTION: Returns a random number drawn from a gaussian (normal)
#              distribution with mean 0.0 and standard deviation 1.0.
#
#              Stolen from 
#
#                 http://www.unix.com.ua/orelly/perl/cookbook/ch02_11.htm
#
#
sub gaussian_rand {
    my ($u1, $u2);  # uniformly distributed random numbers
    my $w;          # variance, then a weight
    my ($g1, $g2);  # gaussian-distributed numbers

    do {
        $u1 = 2 * rand() - 1;
        $u2 = 2 * rand() - 1;
        $w = $u1*$u1 + $u2*$u2;
    } while ( $w >= 1 );

    $w = sqrt( (-2 * log($w))  / $w );
    $g2 = $u1 * $w;
    $g1 = $u2 * $w;
    # could return both if wanted
    #   return wantarray ? ($g1, $g2) : $g1;
    # but I need just one
    return($g1);
}






