#!/usr/bin/perl
#
# Calculate a temperature for a star, given information about it
#    (color, or multiple colors, or spectral type, or all of the above).
#
# The secret is to use data from Eric Mamajek's table of data on stars,
#    which can be found at
#
#    http://www.pas.rochester.edu/~emamajek/EEM_dwarf_UBVIJHK_colors_Teff.txt
#
# Also, calculate the absolute magnitude of the star (using apparent V mag
#    and parallax).
#
# Print one line per star, consisting of the entire original line, 
#    then temperature, then absolute magnitude.
#
# MWR 5/22/2015
# MWR 6/4/2015

use POSIX;

$debug = 0;

# letters in the spectral sequence
$spec_sequence = "OBAFGKMLT";

# the datafile with known stellar properties
$star_data_file = "./EEM_dwarf_UBVIJHK_colors_Teff.txt";


#### step 1: read data from the table of known stellar properties
$nstar = 0;
open(STAR_DATA_FILE, "$star_data_file") || die("can't open file $star_data_file");
while (<STAR_DATA_FILE>) {
  $line = $_;
  if ($line =~ /^#/) {
    next;
  }
  if ($debug > 0) {
    printf "line is %s", $line;
  }

  @words = split(/\s+/, " " . $line);
  $sp_array[$nstar] = $words[1];
  $teff_array[$nstar] = $words[2];
  $bv_array[$nstar] = $words[7];
  $ub_array[$nstar] = $words[9];
  $vi_array[$nstar] = $words[11];

  # use the spectral type to compute a numerical value: spec_numerical
  $spec_numerical = get_spec_numerical($sp_array[$nstar]);
  $spnum_array[$nstar] = $spec_numerical;
 
  if ($debug > 0) {
    printf " star %4d  sp %5s  %6.1f  Teff %6.0f  B-V %6.3f  U-B %6.3f  V-I %6.3f \n",
         $nstar,
         $sp_array[$nstar],
         $spec_numerical,
         $teff_array[$nstar],
         $bv_array[$nstar],
         $ub_array[$nstar],
         $vi_array[$nstar];
  }

  $nstar++;
}
close(STAR_DATA_FILE);



# now we read the _real_ file, full of stars with (B-V) and sptype
#    information ... but not temperatures.  We will try to compute
#    the temperature of each star, then print out one line per
#    star with the original info, plus the temperature
#

$cat_file = "./gliese_318.txt";
open(CAT_FILE, "$cat_file") || die("can't open $cat_file");
while (<CAT_FILE>) {

  $line = $_;
  if ($line =~ /^#/) {
    next;
  }

  # read the spectral type
  $sp_str = substr($line, 68, 4);
     #
     # get rid of any leading 'd'
     $sp_str =~ s/^d//;
  @words = split(/\s+/, $sp_str);
  $sp = $words[0];
  if ($debug > 0) { 
    printf "  sp_str ..%s..  sp ..%s.. \n", $sp_str, $sp;
  }
  # make sure the spectral type is valid; if not, discard the star
  $first_letter = substr($sp, 0, 1); 
  $sp_i = index($spec_sequence, $first_letter);
  if (($first_letter eq "") || ($sp_i == -1)) {
    if ($debug > 0) {
      printf " can't find letter ..%s.., skipping \n", $first_letter;
    }
    next;
  }

  # read the (B-V) color
  $bv_color = -999;
  $bv_str = substr($line, 89, 5);
  @words = split(/\s+/, " " . $bv_str);
  if ($#words > 0) {
    $bv_color = $words[1];
  }
  if ($debug > 0) {
    printf "  sp_str ..%s..   bv_str ..%s..  bv_color %6.3f \n", 
           $sp_str, $bv_str, $bv_color;
  }
  if ($bv_color == -999) {
    if ($debug > 0) {
      printf " can't find B-V color from ..%s.. skipping \n", $bv_str;
    }
    next;
  }

  # read the parallax value.  If we can't find a good one, skip it
  $plx = -1;
  $plx_str = substr($line, 94, 5);
  if ($plx_str <= 0) {
    if ($debug > 0) {
      printf " bad plx_mas value ..%s.. skipping  \n", $plx_str;
    }
    next;
  } 
  $plx = $plx_str*0.001;

  # read the apparent V-mag.
  $vmag = -999;
  $vmag = substr($line, 82, 5);  

  # use parallax and apparent V-mag to compute absolute V-mag
  $dist = 1.0 / $plx;
  $abs_mag = $vmag - (5.0*log10($dist) - 5.0);
  if ($debug > 0) {
    printf "   plx_mas %6.1f  dist %6.1f  vmag %6.3f  absmag %6.3f \n",
         $plx_str, $dist, $vmag, $abs_mag;
  }
   

    
  if (1 == 1) {    
      $temp_from_bv = temp_from_bv($bv_color);
      if ($debug > 0) {
        printf "  (B-V) %6.3f   temp %9.1f \n", $bv_color, $temp_from_bv;
      }
    
      $spec_numerical = get_spec_numerical($sp);
      if ($debug > 1) {
        printf "  sp %6s  spec_numerical %5.2f \n", $sp, $spec_numerical;
      }
      $temp_from_spnum = temp_from_spnum($spec_numerical);
      if ($debug > 0) {
        printf "  spnum %6.3f   temp %9.1f \n", $spec_numerical, $temp_from_spnum;
      }
    
      # now, compute the average temperature of the B-V and spnum estimates
      $num = 0;
      $sum = 0;
      if ($temp_from_bv > 0) {
        $sum += $temp_from_bv;
        $num++;
      }
      if ($temp_from_spnum > 0) {
        $sum += $temp_from_spnum;
        $num++;
      }
    
      if ($num == 0) {
        printf " Q  could not compute temp from colors \n";
      }
      else {
        $temp = $sum / $num;
      }
    
      if ($debug > 0) {
         printf "   color_t %6.0f  sp_t %6.0f  avg_t  %6.0f \n", 
                $temp_from_bv, $temp_from_spnum, $temp;
      }

      # finally, print the original line, with our new temperature
      #     appended at the end
      $new_line = $line;
      chomp($new_line);
      printf " %s     %6.0f  %6.2f  \n", $new_line, $temp, $abs_mag;

   }

    

}
    




exit 0;



####################################################################
# PROCEDURE: get_spec_numerical
#
# DESCRIPTION: Turn a spectral type "A3V" into a numerical value
#              which runs in numerical order from coolest to hottest.
#
# RETURNS:   the numerical value
#
sub get_spec_numerical {

  my($sp_string);
  my($sp_letter, $spec_num, $n);
  my($sp_numerical);

  $sp_string = $_[0];

  $spec_letter = substr($sp_string, 0, 1);
  $spec_num = substr($sp_string, 1, length($sp_string)-1);
  $n = index($spec_sequence, $spec_letter);
  if ($n == -1) {
    printf " error, spec_letter %s not found ?! \n";
    exit(1);
  }
  $n = 10*(length($spec_sequence) - $n);
  $spec_numerical = $n + (10 - $spec_num);
  if ($debug > 1) {
    printf "  spec_letter ..%s..  spec_num ..%s.. n %3d  sp_numerical %6.2f \n", 
           $spec_letter, $spec_num, $n, $spec_numerical;
  }


  return($spec_numerical);
}



####################################################################
# PROCEDURE: temp_from_bv
#
# DESCRIPTION: Given a (B-V) color, interpolate among the 
#              Mamajek stars to estimate a temperature.
#
#              Search for the B-V colors which are closest to the 
#              given B-V color;  then interpolate between their temperatures 
#              to find the  given star's temperature
#
# RETURNS:     temperature         if all goes well
#              -1                  if an error occurs
#
sub temp_from_bv {

  my($bv);
  my($temp_bv, $found);
  my($i);
  my($frac); 

  $bv = $_[0];

   $found = 0;
   $temp_bv = -1;
   for ($i = 0; $i < ($nstar - 1); $i++) {
     if (($bv_array[$i] <= $bv) && ($bv_array[$i+1] > $bv)) {
       if ($debug > 0) {
         printf "   %3d %6.3f   %6.3f    %3d %6.3f  \n",
               $i, $bv_array[$i], $bv, $i+1, $bv_array[$i+1];
       }
       $frac = ($bv - $bv_array[$i]) / ($bv_array[$i+1] - $bv_array[$i]);
       $temp_bv = $teff_array[$i] + $frac*($teff_array[$i+1] - $teff_array[$i]);
       if ($debug > 0) {
         printf "   %3d %6.0f   %6.0f    %3d %6.0f  \n",
               $i, $teff_array[$i], $temp_bv, $i+1, $teff_array[$i+1];
       }
       $found = 1;
       last;
     }
  }
  if ($found == 0) {
    # complain that it wasn't found
    printf " temp_from_bv: couldn't find match for B-V = %6.3f \n", $bv;
    return(-1);
  }

  return($temp_bv);
}


####################################################################
# PROCEDURE: temp_from_spnum
#
# DESCRIPTION: Given a spec_numerical value, interpolate among the 
#              Mamajek stars to estimate a temperature.
#
#              Search for the sp_numerical which are closest to the 
#              given sp_numerical; then interpolate between their temperatures 
#              to find the  given star's temperature
#
# RETURNS:     temperature         if all goes well
#              -1                  if an error occurs
#
sub temp_from_spnum {

  my($spnum);
  my($temp_spnum, $found);
  my($i);
  my($frac); 

  $spnum = $_[0];

   $found = 0;
   $temp_spnum = -1;
   for ($i = 0; $i < ($nstar - 1); $i++) {
     if (($spnum_array[$i] >= $spnum) && ($spnum_array[$i+1] < $spnum)) {
       if ($debug > 0) {
         printf "   %3d %6.3f   %6.3f    %3d %6.3f  \n",
               $i, $spnum_array[$i], $spnum, $i+1, $spnum_array[$i+1];
       }
       $frac = ($spnum - $spnum_array[$i]) / ($spnum_array[$i+1] - $spnum_array[$i]);
       $temp_spnum = $teff_array[$i] + $frac*($teff_array[$i+1] - $teff_array[$i]);
       if ($debug > 0) {
         printf "   %3d %6.0f   %6.0f    %3d %6.0f  \n",
               $i, $teff_array[$i], $temp_spnum, $i+1, $teff_array[$i+1];
       }
       $found = 1;
       last;
     }
  }
  if ($found == 0) {
    # complain that it wasn't found
    printf " temp_from_spnum: couldn't find match for spnum = %6.3f \n", $spnum;
    return(-1);
  }

  return($temp_spnum);
}
