#!/usr/bin/perl
#
# Read information from a datafile containing sizes and distances
#    of planets and moons, then compute the apparent angular
#    size of Sun and (best) moon as seen from each planet.
#
# MWR 8/15/2017

use POSIX;
use Class::Struct;

struct ( Body_orb => { name => '$', 
                   number => '$', 
                   parent => '$', 
                   distance => '$', 
                   period => '$', 
                   incl => '$', 
                   eccen => '$', 
                   phys_radius => '$', 
                   is_planet => '$', 
                   nomatch => '$',
                   sun_ang_size => '$',
                   bestmoon_ang_size => '$',
                   bestmoon_name => '$',
                   bestmoon_ratio => '$',
                   xxx => '$' });
my @body_orb_array = ();

struct ( Body_phys => { name => '$', 
                   radius => '$', 
                   mass => '$', 
                   dens => '$', 
                   parent => '$', 
                   sun_ang_size => '$',
                   xxx => '$' });
my @body_phys_array = ();


# convert from radians to arcminutes
$RADTOARCMIN = (180.0/3.14159)*60.0;
# Sun's radius in meters
$solar_radius_m = 6.95E8;


$debug = 0;

# file containing orbital information
$orbital_datafile = "./orbital_table.dat";
# file containing physical information
$physical_datafile = "./physical_table.dat";

# first, read in all orbital info
$num_body_orb = 0;
open(DATAFILE, "$orbital_datafile") || die("can't read file $orbital_datafile");
while (<DATAFILE>) {

  $line = $_;
  # skip comments and blank lines
  if ($line =~ /^#/) {
    next;
  }
  if ($line =~ /^$/) {
    next;
  }

  @words = split(/\s+/, $line);
  if ($#words < 7) {
    printf " bad line has %d words, skipping .. \n", $#words;
    printf "%s", $line;
    next;
  }

  if ($debug > 0) {
    printf " next good line is: %s", $line;
  }

  $name = $words[0];
  $number = $words[1];
  $parent = $words[2];
  $distance = $words[3];
  $period = $words[4];
  $incl = $words[5];
  $eccen = $words[6];

  if ($debug > 0) {
    printf "  name ..%s..   number ..%s..  parent ..%s.. \n",
          $name, $number, $parent;
  }

  $body_orb_array[$num_body_orb] = Body_orb->new(name => $name);
  $body_orb_array[$num_body_orb]->number($number);
  $body_orb_array[$num_body_orb]->parent($parent);
  $body_orb_array[$num_body_orb]->distance($distance);
  $body_orb_array[$num_body_orb]->period($period);
  $body_orb_array[$num_body_orb]->incl($incl);
  $body_orb_array[$num_body_orb]->eccen($eccen);
  $body_orb_array[$num_body_orb]->nomatch(0);
  $body_orb_array[$num_body_orb]->is_planet(0);
  $body_orb_array[$num_body_orb]->sun_ang_size(-1);
  $body_orb_array[$num_body_orb]->bestmoon_ang_size(-1);

  $num_body_orb++;
}
close(DATAFILE);

if ($debug > 0) {
  printf " read %5d bodies \n", $num_body_orb;
  for (my $i = 0; $i < $num_body_orb; $i++) {
    printf " body %5d is %20s %5s %20s \n", $i,
           $body_orb_array[$i]->name,
           $body_orb_array[$i]->number,
           $body_orb_array[$i]->parent;
  }
}
 


# next, read in all physical info
$num_body_phys = 0;
open(DATAFILE, "$physical_datafile") || die("can't read file $physical_datafile");
while (<DATAFILE>) {

  $line = $_;
  # skip comments and blank lines
  if ($line =~ /^#/) {
    next;
  }
  if ($line =~ /^$/) {
    next;
  }

  @words = split(/\s+/, $line);
  if ($#words < 3) {
    printf " bad line has %d words, skipping .. \n", $#words;
    printf "%s", $line;
    next;
  }

  if ($debug > 0) {
    printf " next good line is: %s", $line;
  }

  $name = $words[0];
  $radius = $words[1];
  $mass = $words[2];
  $dens = $words[3];

  if ($debug > 0) {
    printf "  name ..%s..   radius ..%s..  mass ..%s.. \n",
          $name, $radius, $mass;
  }

  $body_phys_array[$num_body_phys] = Body_phys->new(name => $name);
  $body_phys_array[$num_body_phys]->radius($radius);
  $body_phys_array[$num_body_phys]->mass($mass);
  $body_phys_array[$num_body_phys]->dens($dens);
  $body_phys_array[$num_body_phys]->parent("");

  $num_body_phys++;
}
close(DATAFILE);

if ($debug > 0) {
  printf " read %5d bodies \n", $num_body_phys;
  for (my $i = 0; $i < $num_body_phys; $i++) {
    printf " body %5d is %20s %5s %20s \n", $i,
           $body_phys_array[$i]->name,
           $body_phys_array[$i]->radius,
           $body_phys_array[$i]->mass;
  }
}
 

# Go through the two lists of information and look for matching objects.
#   Transfer useful information from the "body_phys" to the "body_orb"
#   structures
for ($i_orb = 0; $i_orb < $num_body_orb; $i_orb++) {

  $found = 0;
  for ($i_phys = 0; $i_phys < $num_body_phys; $i_phys++) {
    if ($body_orb_array[$i_orb]->name eq $body_phys_array[$i_phys]->name) {

      # copy the physical radius of the body to the body_orb_array[]
      #    entry
      $body_orb_array[$i_orb]->phys_radius($body_phys_array[$i_phys]->radius);

      $found = 1;
    }

  }

  if ($found == 0) {
    printf " warning: could not match i_phys %5d = %s -- disqualifying \n", $i_orb,
            $body_orb_array[$i_orb]->name;
    $body_orb_array[$i_orb]->nomatch(1);
  }

}


# Loop over all the bodies in the orbital array, to decide if each one
#    is a planet (orbits Sun) or a moon (orbits other body).
#    Set the "is_planet" field appropriately
for ($i_orb = 0; $i_orb < $num_body_orb; $i_orb++) {
  if ($body_orb_array[$i_orb]->parent eq "Sun") {
    $body_orb_array[$i_orb]->is_planet(1);
    if ($debug > 0) {
      printf " body %5d = %20s is a planet \n", 
           $i_orb, $body_orb_array[$i_orb]->name;
    }
  }
}
  




# now, compute the Sun's angular size as seen from 
#   each planet (units: minutes of arc),
#   and save in the "sun_ang_size" field
for ($i = 0; $i < $num_body_orb; $i++) {

  if ($body_orb_array[$i]->is_planet == 1) {

    # compute distance to Sun in meters
    $orb_dist_m = $body_orb_array[$i]->distance*1.0e6;

    # compute angular radius of Sun (radians)
    $sun_ang_radius = atan2($solar_radius_m, $orb_dist_m);
    # convert from radians to arcmin
    $sun_ang_radius *= $RADTOARCMIN;

    $body_orb_array[$i]->sun_ang_size($sun_ang_radius);
    if ($debug > 0) {
      printf " body %5d = %20s has sun ang radius %9.3lf arcmin \n",
            $i, $body_orb_array[$i]->name, $body_orb_array[$i]->sun_ang_size;
    }

  }

}



# For each planet, loop over all the moons associated with it.
#
#   For each moon, compute the angular radial size (arcmin) 
#   of the moon, as seen from the planet's surface/cloudtops.
#  
#   Compute the ratio of angular size of Sun to angular size of moon,
#   and find the moon which yields the ratio closest to 1.0.
#   Save the size and name of that "best" moon in the planet's
#   fields.
#
for ($i = 0; $i < $num_body_orb; $i++) {

  if ($body_orb_array[$i]->is_planet == 1) {
  
    $this_planet = $body_orb_array[$i]->name;
    $this_sun_size = $body_orb_array[$i]->sun_ang_size;
    if ($debug > 0) {
      printf "  looking for moons of planet %s Sun_size %7.3f \n", 
             $this_planet, $this_sun_size;
    }

    $best_moon_index = -1;
    $best_moon_ratio = 0.0;
    for ($i_moon = 0; $i_moon < $num_body_orb; $i_moon++) {

      if ($body_orb_array[$i_moon]->parent eq $this_planet) {

        $moon_dist_m = $body_orb_array[$i_moon]->distance*1.0e6 
                        -  $body_orb_array[$i]->phys_radius*1.0e3;

        $moon_radius_m = $body_orb_array[$i_moon]->phys_radius*1.0e3;

        $moon_ang_radius = atan2($moon_radius_m, $moon_dist_m);
        $moon_ang_radius *= $RADTOARCMIN;

        if ($debug > 0) {
          printf "     moon %20s has ang radius %12.9f arcmin \n", 
               $body_orb_array[$i_moon]->name, $moon_ang_radius;
        }

        if ($moon_ang_radius <= 0) {
          printf " skipping moon %20s because radius is zero \n", 
               $body_orb_array[$i_moon]->name;
          next;
        }

        $this_ratio = $this_sun_size / $moon_ang_radius;
        if ($this_ratio > 1.0) {
          $this_ratio = 1.0 / $this_ratio;
        }

        if ($this_ratio > 0.8) {
          if ($debug > 0) {
            printf " planet %12s  moon %5d = %20s has good ratio %7.3f \n", 
               $this_planet, $i_moon, $body_orb_array[$i_moon]->name,
               $this_ratio;
          }
        }


        if ($this_ratio > $best_moon_ratio) {
          $best_moon_ratio = $this_ratio;
          $best_moon_index = $i_moon;

          $body_orb_array[$i]->bestmoon_name($body_orb_array[$i_moon]->name);
          $body_orb_array[$i]->bestmoon_ratio($this_sun_size / $moon_ang_radius);
          $body_orb_array[$i]->bestmoon_ang_size($moon_ang_radius);

          if ($debug > 0) {
            printf "  moon %5d = %20s has best ratio %7.3f \n", 
               $i_moon, $body_orb_array[$i_moon]->name,
               $this_ratio;
          }
        }

      }

    }

  }

}


###########
# Finally, we print out the results: best moon for each planet
#
for ($i = 0; $i < $num_body_orb; $i++) {

  if ($body_orb_array[$i]->is_planet == 1) {
  
    $this_planet = $body_orb_array[$i]->name;
    $this_sun_size = $body_orb_array[$i]->sun_ang_size;
    $bestmoon_name = $body_orb_array[$i]->bestmoon_name;
    $bestmoon_ratio = $body_orb_array[$i]->bestmoon_ratio;
    $bestmoon_ang_size = $body_orb_array[$i]->bestmoon_ang_size;

    printf " %15s  Sun %9.5f  %15s %9.5f  ratio %7.3f \n",
         $this_planet, $this_sun_size, $bestmoon_name, $bestmoon_ang_size,
         $bestmoon_ratio;


  }
}


##############
# Make another set of checks, looking for eclipses which could occur
#    when one moon of a planet passes in front of the Sun, as seen
#    from a second moon of the same planet.
#
for ($i = 0; $i < $num_body_orb; $i++) {

  if ($body_orb_array[$i]->is_planet == 1) {
  
    $this_planet = $body_orb_array[$i]->name;
    $this_sun_size = $body_orb_array[$i]->sun_ang_size;
    if ($debug > 0) {
      printf "  looking for moon-moon events of planet %s Sun_size %7.3f \n", 
             $this_planet, $this_sun_size;
    }

    $best_moon_index = -1;
    $best_moon_ratio = 0.0;
    for ($i_moon = 0; $i_moon < $num_body_orb; $i_moon++) {

      if ($body_orb_array[$i_moon]->parent eq $this_planet) {

        $at_moon_name = $body_orb_array[$i_moon]->name;

        if ($debug > 0) {
          printf "   consider looking AT moon %s \n", $at_moon_name;
        }

        $at_moon_radius_m = $body_orb_array[$i_moon]->phys_radius*1.0e3;
        $moon_distance_m = $body_orb_array[$i_moon]->distance*1.0e6;
        $at_moon_incl = $body_orb_array[$i_moon]->incl;

        # figure out the distance from which this moon would subtend the 
        #   same angular radius as the Sun-as-seen-from-this-planet
        $angle_radians = $this_sun_size/$RADTOARCMIN;
        $want_distance_m = $at_moon_radius_m / tan($angle_radians);
        if ($debug > 0) { 
          printf "       want to be at distance %9.4e meters \n", 
                   $want_distance_m;
        }

        # now loop over all the other moons of this planet, checking to see
        #    if they ever reach this wanted distance from the moon
        #    under consideration
        for ($j_moon = 0; $j_moon < $num_body_orb; $j_moon++) {
 
          # skip the at-moon itself
          if ($j_moon == $i_moon) {
            next;
          }

          if ($body_orb_array[$j_moon]->parent eq $this_planet) {

            $from_moon_name = $body_orb_array[$j_moon]->name;

            if ($debug > 0) {
              printf "      consider looking FROM moon %s \n", $at_moon_name;
            }

            $from_moon_radius_m = $body_orb_array[$j_moon]->phys_radius*1.0e3;
            $from_moon_distance_m = $body_orb_array[$j_moon]->distance*1.0e6;
            $from_moon_incl = $body_orb_array[$j_moon]->incl;

            $min_dist = abs($from_moon_distance_m - $moon_distance_m);
            $max_dist = $from_moon_distance_m + $moon_distance_m;

            if (($want_distance_m > $min_dist) && 
                    ($want_distance_m < $max_dist)) {
              printf " planet %12s Sun %7.3f from %12s to %12s %6.0f dist %9.4e incls %2.0f %2.0f \n",
                     $this_planet, $this_sun_size, $from_moon_name,
                     $at_moon_name, $at_moon_radius_m/1.0e3, $want_distance_m,
                     $from_moon_incl, $at_moon_incl;
            } 

          }  # is the from-moon circling correct planet?

        }  # end of loop over other moon

      } # this the at-moon circling correct planet?

    } # end of loop over at-moon

  } # is this body a planet?

} # end of loop over all bodies
 



exit 0;
