# Compute the signals which would be seen by an astronomer
#   looking at a system with 
# 
#       a solar-twin star
#       a Jupiter-twin planet in circular orbit at 1 AU 
#
# We compute the radial velocity signal and the transit signal,
#   assuming one measurement per night.
#
# MWR 9/29/2016

$debug = 1;

$mypi = 3.14159;

# masses in kg
$mass_star = 1.99E30;
$mass_planet = 1.89E27;
# radii in m
$radius_star = 6.96E8;
$radius_planet = 7.15E7;
# orbital radius in m
$orb_radius = 1.5E11;
# orbital period in seconds
$orb_period = 3.15E7;

# this is the (approximate) radius of star's reflex motion
$d = ($mass_planet/$mass_star) * $orb_radius;
# speed of planet (m/s)
$speed_planet = 29800;
# speed of star
$speed_star = (2.0*$mypi*$d) / ($orb_period);

# depth of transit
$transit_depth = ($radius_planet*$radius_planet) / ($radius_star*$radius_star);
# duration of transit (sec)
$duration_transit = (2.0*$radius_star) / $speed_planet;


# planet starts at western elongation at time t = 0,
#   moves toward Earth at start.
$beta = ($radius_star + $radius_planet) / $orb_radius;
# the _approximate_ start of transit is when angle theta reaches this value
$start_transit_theta = ($mypi/2.0) - $beta;
# and this is the end of the transit
$end_transit_theta = ($mypi/2.0) + $beta;


# loop over time, go by days
$start_day = 0;
$end_day = 750;
for ($i_day = $start_day; $i_day < $end_day; $i_day++) {

  # time in seconds since the start
  $t = $i_day*86400;

  # angle of planet in its orbit (radians)
  $theta = 2*$mypi*($t / $orb_period);

  # radial velocity of star (m/s)
  $rv = (1.0)*$speed_star*cos($theta);

  # brightness of star -- we simply determine in/out of transit
     # first, we force the angle into range 0 - 2pi
     $theta_fix = $theta;
     while ($theta_fix > 2*$mypi) {
       $theta_fix -= 2.0*$mypi;
     }
  if (($theta_fix >= $start_transit_theta) && 
                ($theta_fix <= $end_transit_theta)) {
    $flux = 1.0 - $transit_depth;
  } else {
    $flux = 1.0;
  }
  


  # print one line per timestep: time(days)  time(sec)  RV   flux 
  $t_day = $t / 86400.0;
  printf "%-10.3f %12.3f    %7.2f  %12.5f  \n",
      $t_day, $t, $rv, $flux;


}




