#!/usr/bin/perl
#
# Compute the variations in measured radial velocity due to
#     a) the spin of the Earth
#     b) the orbital motion of the Earth
#
# for a "typical" star as seen from a "typical" observatory.

use POSIX;

$debug = 1;

$mypi = 3.1415926;
$DEGTORAD = ($mypi / 180.0);


# equatorial spin velocity of earth (m/s)
$spin_v = 465.0;
# latitude of observatory (degrees)
$obs_lat = 40.0;
$obs_lat_rad = $obs_lat*$DEGTORAD;
# hour angle of target at midnight (degrees) 
$ha_at_midnight = 0;
$ha_at_midnight_rad = $ha_at_midnight*$DEGTORAD;
# time at midnight (days)
$jd_at_midnight = 0.0;
# declination of star (degrees)
$star_dec = 30.0;
$star_dec_rad = $star_dec*$DEGTORAD;
# length of sidereal day, in days
$sidereal_length = 23.9344696/24.0;

# max orbital speed of Earth (m/s)
$orb_v = 29800;
# orbital period of Earth (days)
$orb_period = 365.25;
# ecliptic latitude of the target (degrees)
$ecl_lat = 30.0;
$ecl_lat_rad = $ecl_lat*$DEGTORAD;
# ecliptic longitude of Earth at starting time (degrees)
$ecl_long_earth_start = 20.0;
$ecl_long_earth_start_rad = $ecl_long_earth_start*$DEGTORAD;
# ecliptic longitude of star (degrees)
$ecl_long_star = 20.0;
$ecl_long_star_rad = $ecl_long_star*$DEGTORAD;


$start_jd = $jd_at_midnight;
$end_jd = $start_jd + 600;
$d_jd = 0.05;


# note we'll use time in days
for ($jd = $start_jd; $jd <= $end_jd; $jd += $d_jd) {

  # compute spin component of radial velocity here
  $ha_cycles = ($jd - $start_jd) / $sidereal_length;
  $ha_rad = $ha_cycles*2.0*$mypi;
  $rv_spin = $spin_v*cos($star_dec_rad)*sin($ha_rad);


  # compute orbital component of radial velocity here
  $orb_cycles = ($jd - $start_jd) / ($orb_period);
  $orb_rad = $orb_cycles*2.0*$mypi;
  $ecl_long_earth_rad = $ecl_long_earth_start_rad + $orb_rad;
  $az_angle_rad = $ecl_long_earth_rad - $ecl_long_star_rad;
  $rv_orb = $orb_v*cos($ecl_lat_rad)*sin($az_angle_rad);
 

 
  # compute total radial velocity
  $rv_tot = $rv_spin + $rv_orb;


  # print results
  printf " %12.5f   %10.4f  %10.4f   %10.4f   \n", 
        $jd, $rv_spin, $rv_orb, $rv_tot;




}



exit 0;



