#!/usr/bin/perl
#
# Read data from a table of stellar isochrones.  Use the data
#    to compute a grid of absolute magnitudes of combined light
#    for all possible mass combinations in a binary star pair.
#
# Print the resulting to as ASCII table to output,
#    and create a color-map graphic with the results.
#
# MWR 4/12/2022

use POSIX;

$debug = 0;
$infile = "./isochrones_a.csv";

# only choose stars with this age
$want_age = 1.0E9;
# only include entries with an index less than this value
#    stars with larger indices have evolved off the main sequence
$max_index = 458;

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

  @words = split(/,/, $line);
  $index = $words[1];
  $log_age = $words[2];
  $initial_mass = $words[3];
  $current_mass = $words[4];
  $log_Teff = $words[5];
  $log_g = $words[6];  
  $log_L = $words[7];
  $fe_h_initial = $words[8];
  $fe_h_current = $words[9];
  $bess_U = $words[10];
  $bess_B = $words[11];
  $bess_V = $words[12];
  $bess_R = $words[13];
  $bess_I = $words[14];
  

  $age = pow(10.0, $log_age);
  $Teff = pow(10.0, $log_Teff);

  if ($age != $want_age) {
    next;
  }
  if ($index > $max_index) {
    next;
  }

  # add this star to our array of possible binary components
  $index_array[$num] = $index;
  $age_array[$num] = $age;
  $mass_array[$num] = $current_mass;
  $teff_array[$num] = $Teff;
  $absmagB_array[$num] = $bess_B;
  $absmagV_array[$num] = $bess_V;

  if ($debug > 0) {
    printf " %5d  %10.5e  %9.4f  %6.0f %6.3f %6.3f   \n", 
         $index, $age, $current_mass, $Teff, $bess_B, $bess_V;
  }

  $num++;

}
close(INFILE);
if ($debug > 0) {
  printf " read a total of %6d main sequence stars \n", $num;
}


# now, walk through the list of stars.  
#   For each star, walk through the list again, choosing each star
#   in turn as a companion.
#
#     For each choice of pairs, compute the absolute magnitude 
#     of the light in B and V.  Print out one line per pair
#     with info on that pair.
#
for ($i_star = 0; $i_star < $num; $i_star++) {

  $inten_B_i = pow(10.0, (-0.4*$absmagB_array[$i_star]));
  $inten_V_i = pow(10.0, (-0.4*$absmagV_array[$i_star]));

  for ($j_star = 0; $j_star < $num; $j_star++) {

    $inten_B_j = pow(10.0, (-0.4*$absmagB_array[$j_star]));
    $inten_V_j = pow(10.0, (-0.4*$absmagV_array[$j_star]));

    $tot_inten_B = $inten_B_i + $inten_B_j;
    $tot_absmagB = -2.5*log10($tot_inten_B);
    $tot_inten_V = $inten_V_i + $inten_V_j;
    $tot_absmagV = -2.5*log10($tot_inten_V);
    $tot_BminusV = $tot_absmagB - $tot_absmagV;


    printf "  %3d %8.4f %6.0f %6.3f %6.3f   %3d %8.4f %6.0f %6.3f %6.3f   %6.3f %6.3f  %5.3f \n",
          $index_array[$i_star], $mass_array[$i_star], $teff_array[$i_star],
          $absmagB_array[$i_star], $absmagV_array[$i_star],
          $index_array[$j_star], $mass_array[$j_star], $teff_array[$j_star],
          $absmagB_array[$j_star], $absmagV_array[$j_star],
          $tot_absmagB, $tot_absmagV, $tot_BminusV;

  }

  # add a blank line for gnuplot
  printf "\n";

}

exit 0;

