#!/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

use POSIX;

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

# the base date for calculations is Jan 0, 2000.0
$jd_base = 2451543.5;
# this corresponds to May 26, 2020 = JD 2458995.0;
$start_date = 2458995.5;
# this is several years later
$end_date = $start_date + (365.25*3.0);
# time step (in days)
$d_date = 10.0;

# obliquity of Earth's orbit (degrees)
$obliq = 23.440972;
$cos_obliq = cos($obliq*$DEGTORAD);
$sin_obliq = sin($obliq*$DEGTORAD);
# 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)
$mu_ra = 2891.5;
# proper motion in Dec (milliarcsec per year)
$mu_dec = 411.9;
# parallax (milliarcsec)
$parallax = 280.7;

# cosine of declination of star
$cosdec = cos($DEGTORAD*$dec_start);

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

for ($date = $start_date - 1000; $date <= $end_date; $date += $d_date) {

  # days since $start_date
  $day_offset = $date - $start_date;
  # years since $start_date
  $year_offset = ($day_offset)/365.25;

  # compute the solar ecliptic longitude for this date
  #   days since Jan 0, 2000
  $d = $date - $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;

    printf " 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;

    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;

    printf " x1 %8.2f x3 %8.2f delta_dec %8.2f ", 
            $x1, $x3, $delta_dec;

    printf "\n";

}


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