#!/usr/bin/perl
# Compute the apparent angular sizes of Sun and Moon as seen from Earth
#    over the next 5 billion years or so
#    Include the following effects:
#
#      a) Moon's recession from Earth due to tidal effects
#          assumed constant at 3.8 cm per year
#
#      b) Sun's increase in radius due to stellar evolution
#
#    Assume that the Moon's orbit retains its current eccentricity
#    during that period, simply changing semi-major axis.
#
#    Ignore slight increase in Earth's orbital radius -- is too small.
#
# MWR 8/19/2017

use POSIX;
use Class::Struct;

$debug = 0;

$RADTOARCMIN = 60*(180.0/3.14159);


struct ( Solar => { age => '$', 
                    teff => '$', 
                    logg => '$', 
                    logL => '$', 
                    absV => '$', 
                    rad => '$', 
                    rad_ratio => '$', 
                    xxx => '$'  } );
my @solar_array = ();
my $num_solar = 0;


# file containing information solar radius vs. time
$solar_radius_file = "./calc_radius.out";
# the distance from Sun to Earth, assumed constant (m)
$sun_dist = 1.496E11;
# the current age of Sun (years)
$current_solar_age = 4.64E9;
# eccentricity of Earth's orbit (assumed fixed)
$earth_eccen = 0.0167;
$earth_eccen = 0.058;
# eccentricity of Moon's orbit (assumed fixed)
$moon_eccen = 0.0549;
$moon_eccen = 0.0700;
# current value of Moon's orbital semi-major axis (meters)
$moon_a_now = 3.8475E8;
# rate at which Moon's orbital radius grows (meters per year)
$moon_a_rate = 0.038;
# radius of Moon (meters)
$moon_rad = 1738.0E3;





# read information from the file with solar radius vs. time,
#    save information in arrays for interpolation
if (read_solar_data($solar_radius_file) != 0) {
  printf STDERR "error: read_solar_data fails for file $infile \n";
  exit(1);
}



$start_time = 0;
$end_time = 1.0e9;
$time_step = 0.5e7;
for ($time = $start_time; $time <= $end_time; $time += $time_step) {

  # compute the physical radius of the Sun at this time
  $sun_rad = interp_solar_radius($time);
  if ($debug > 0) {
    printf "     %9.4e  Sun_rad   %12.5e \n", $time, $sun_rad;
  }

  # compute the angular size of the Sun
  $sun_min_ang_size_rad = atan2($sun_rad, (1.0 + $earth_eccen)*$sun_dist);
  $sun_min_ang_size_arcmin = $sun_min_ang_size_rad*$RADTOARCMIN;
  $sun_max_ang_size_rad = atan2($sun_rad, (1.0 - $earth_eccen)*$sun_dist);
  $sun_max_ang_size_arcmin = $sun_max_ang_size_rad*$RADTOARCMIN;

  # compute the orbital radius of the Moon at this time
  $moon_dist = $moon_a_now + $moon_a_rate*$time;
  
     $moon_max_dist = $moon_dist*(1.0 + $moon_eccen);
     $moon_min_dist = $moon_dist*(1.0 - $moon_eccen);
     if ($debug > 0) {
       printf "  moon extremes  %12.6e %12.6e \n", $moon_min_dist, $moon_max_dist;
     }


  # compute the angular sizes of the Moon
     # first, for apogee
     $moon_min_ang_size_rad = atan2($moon_rad, (1 + $moon_eccen)*$moon_dist);
     $moon_min_ang_size_arcmin = $moon_min_ang_size_rad*$RADTOARCMIN;
     # next, for perigee
     $moon_max_ang_size_rad = atan2($moon_rad, (1 - $moon_eccen)*$moon_dist);
     $moon_max_ang_size_arcmin = $moon_max_ang_size_rad*$RADTOARCMIN;

  # now, print one line with summary of information
  printf " %9.4e  sunrad %12.5e  %6.3f %6.3f   moon_dist %12.5e  %6.3f %6.3f \n",
      $time, $sun_rad, $sun_min_ang_size_arcmin, $sun_max_ang_size_arcmin,
      $moon_dist, $moon_min_ang_size_arcmin, $moon_max_ang_size_arcmin;



}


exit 0;


######################################################################
# Read information about the change in solar radius with time
#   into an array of structures, which we'll use for interpolation
#   to arbitrary times.
#
# 
sub read_solar_data {
  my $infile = $_[0];

  my ($line, @words);
  my ($age, $teff, $logg, $logL, $absV, $rad, $rad_ratio);

  open(INFILE, "$infile") || die("read_solar_data: can't open $infile");
  while (<INFILE>) {
    $line = $_;
    if ($line =~ /^#/) {
      next;
    }

    @words = split(/\s+/, " " . $line);
    $age = $words[1];
    $teff = $words[2];
    $logg = $words[3];
    $logL = $words[4];
    $absV = $words[5];
    $rad = $words[6];
    $rad_ratio = $words[7];
    if ($debug > 0) {
      printf " read age %9.4e  teff %7.0f  rad_ratio %6.3f \n",
          $age, $teff, $rad_ratio;
    }

    $solar_array[$num_solar] = Solar->new();
    $solar_array[$num_solar]->age($age);
    $solar_array[$num_solar]->teff($teff);
    $solar_array[$num_solar]->logg($logg);
    $solar_array[$num_solar]->logL($logL);
    $solar_array[$num_solar]->absV($absV);
    $solar_array[$num_solar]->rad($rad);
    $solar_array[$num_solar]->rad_ratio($rad_ratio);

    $num_solar++;
  }

  if ($debug > 0) {
    printf " read %5d items from %s \n", $num_solar, $infile;
  }


  return(0);
}



######################################################################
# PROCEDURE: interp_solar_radius
#
# DESCRIPTION: Given a time (in years from the present time),
#              interpolate into the @solar_array to figure out
#              the radius of the Sun at that time, in meters.
#
# RETURNS:     radius of sun               if all goes well
#              -1                          if an error occurs
#
# 
sub interp_solar_radius {

  my ($age) = $_[0];
  my ($prev_age, $next_age);
  my ($prev_rad, $next_rad, $this_rad);
  my ($frac);

  if ($age < 0) {
    printf "interp_solar_radius: invalid age %s \n", $age;
    return(-1);
  }

  # add the current age of the Sun to the given number of years,
  #   so that it represents the age of the Sun from its birth
  $age += $current_solar_age;

  if ($debug > 0) {
    printf " inter_solar_radius: looking for age %9.4e \n", $age;
  }
  
  # find the two elements of the solar_array which bracket 
  #   the desired time
  for (my $i = 0; $i < $num_solar - 1; $i++) {
    $prev_age = $solar_array[$i]->age;
    $next_age = $solar_array[$i + 1]->age;
    if ($debug > 0) {
      printf "  i %6d  prev_age %9.4e  next_age %9.4e \n",
           $i, $prev_age, $next_age;
    }

    if (($age >= $prev_age) && ($age < $next_age)) {

      $frac = ($age - $prev_age) / ($next_age - $prev_age);
      $prev_rad = $solar_array[$i]->rad;
      $next_rad = $solar_array[$i + 1]->rad;
      $this_rad = $prev_rad + $frac*($next_rad - $prev_rad);

      return($this_rad);
    }
  }
  
  # if we get here, we could not find the right spot to interpolate
  printf STDERR "interp_solar_radius: fell out of loop for age %9.5e \n",
      $age;
  return(-1);
}

