#!/usr/bin/perl
#
# Compute a simple model of the motion of a star,
#    with combination of proper motion and orbital motion around CoM


$TWOPI = 2.0*3.14159;
# amplitude of orbital motion of star (millarcsec)
$star_orb_radius = 10.0;
# amplitude of orbital motion of planet (milliarcsec)
$planet_orb_radius = 100.0;
# period of orbital motion (years)
$period = 20.0;
# proper motion in X-dir (milliarcsec per year)
$mu_x = 5.0;
# proper motion in Y-dir (milliarcsec per year)
$mu_y = 0.0;
# starting position of center of mass in X-dir
$start_x = 0.0;
# starting position of center of mass in Y-dir
$start_y = 15.0;
# starting time (years)
$start_time = 0.0;
# ending time (years)
$end_time = 100.0;
# timestep (years)
$dt = 0.2;

for ($time = $start_time; $time <= $end_time; $time += $dt) {

  $sox = $star_orb_radius*cos( ($TWOPI*$time)/$period );
  $soy = $star_orb_radius*sin( ($TWOPI*$time)/$period );

  $pox = $planet_orb_radius*(0 - cos( ($TWOPI*$time)/$period ));
  $poy = $planet_orb_radius*(0 - sin( ($TWOPI*$time)/$period ));

  $star_x = $start_x + ($mu_x*$time) + $sox;
  $star_y = $start_y + ($mu_y*$time) + $soy;

  $planet_x = $start_x + ($mu_x*$time) + $pox;
  $planet_y = $start_y + ($mu_y*$time) + $poy;

  printf " t %8.3f   sx %10.4f  sy %10.4f  px %10.4f  py %10.4f \n", 
      $time, $star_x, $star_y, $planet_x, $planet_y;

}


