#!/usr/bin/perl
#
# Given a set of N files, each of which contains a spectrum with format
#
#      lambda   inten
#
# we want to add the spectra together.  The tricky thing is that the
# wavelengths won't necessarily match.
#
# So, the algorithm is simple and stupid.
#
#   1. pick the first given spectrum as fiducial.  We'll produce output
#         with its wavelengths.  Read in all data from this file
#
#   2. for every other spectrum (loop)
#
#       3. read in all data from new spectrum
#
#       4. walk through list of fiducial spectrum wavelengths.  
#          for each wavelength
#
#             5. find wavelengths in new spectrum which bracket it
#                      (if no bracketing pair, skip it)
#
#             6. linearly interpolate to find value of new spectrum
#                      at location of fiducial spectrum
#
#             7. add interpolated value of new spectrum to total intensity
#                      in the summed output spectrum
#
#        end loop
#
#   end loop
#
#
# usage:   perl add_spec.pl spec1.dat spec2.dat spec3.dat ....
#
#
# MWR 3/21/2011

$debug = 0;


$nfiles = $#ARGV + 1;
if ($nfiles < 2) {
  printf STDERR "usage: perl add_spec.pl spec1.dat spec2.dat .... \n";
  exit(1);
}
if ($debug > 0) {
  for ($i = 0; $i < $nfiles; $i++) {
    printf " nfiles [%d] is ..%s.. \n", $i, $ARGV[$i];
  }
}

# read in all data from all datafiles
for ($i = 0; $i < $nfiles; $i++) {
  open(INFILE, "$ARGV[$i]") || die("can't open file %s \n");
  $num_lines[$i] = 0;
  while (<INFILE>) {
    $line = $_;
    if ($line =~ /^#/) {
      next;
    }
    if ($line =~ /^$/) {
      next;
    }
    
    $line = " " . $line;
    @words = split(/\s+/, $line);
    $lambda = $words[1];
    $inten = $words[2];
    if ($debug > 0) {
      printf " lambda %s   inten %s \n", $lambda, $inten;
    }

    $lambda_array[$i][$num_lines[$i]] = $lambda;
    $inten_array[$i][$num_lines[$i]] = $inten;

    $num_lines[$i]++;

  }
  if ($debug > 0) {
    printf " read %d lines from file %d = %s \n", 
        $num_lines[$i], $i, $ARGV[$i];
  }

}


# Now, we add up the spectra.  Use wavelengths of the
#   fiducial, first datafile.  Interpolate in others to find intensity
#   at that wavelength and place sum in sum_array[].
#
for ($i = 0; $i < $num_lines[0]; $i++) {
  $sum_lambda_array[$i] = $lambda_array[0][$i];
  $sum_inten_array[$i] = $inten_array[0][$i];

  if ($debug > 0) {
    printf " afore adding, sum lambda %8.1f   =  %9.4e \n",
            $sum_lambda_array[$i], $sum_inten_array[$i];
  }

  for ($i_file = 1; $i_file < $nfiles; $i_file++) {

     $interp_inten = interpolate_inten($i_file, $lambda_array[0][$i]);
     if ($interp_inten != -1) {
       $sum_inten_array[$i] += $interp_inten;
     }

  }
  if ($debug > 0) {
    printf " after adding, sum lambda %8.1f   =  %9.4e \n",
            $sum_lambda_array[$i], $sum_inten_array[$i];
  }

  # divide the sum by the number of spectra added together
  $sum_inten_array[$i] /= $nfiles;

  printf " %8.1f  %9.4e \n", $sum_lambda_array[$i], $sum_inten_array[$i];

}






exit 0;


######################################################################
# PROCEDURE: interpolate_inten
#
# DESCRIPTION: Given an index number, and a wavelength, 
#              we look for the wavelength entries in the index'th
#              datafile which surround the given wavelength.
#              We interpolate linearly to find the intensity 
#              at the given wavelength.  Return it.
#
#              If we can't find values on either side of the
#              given wavelength, return -1.0
#
# RETURNS:
#                  interpolated intensity at given wavelength, if can be found
#                  -1 if not
#
sub interpolate_inten {

 my($index, $fid_lambda);
 my($i, $prev_lambda, $next_lambda, $prev_inten, $next_inten);
 my($x, $frac, $this_inten);

   $index = $_[0];
   $fid_lambda = $_[1];
  
   if ($debug > 0) {
     printf " enter interpolate_inten for file %d lambda %8.1f \n",
                 $index, $fid_lambda;
   }

   for ($i = 0; $i < $num_lines[$index] - 1; $i++) {
     if (($lambda_array[$index][$i] <= $fid_lambda) &&
         ($lambda_array[$index][$i+1] > $fid_lambda)) {
       # okay, found the pair which span the desired wavelength
       $frac = ($fid_lambda - $lambda_array[$index][$i]) /
                   ($lambda_array[$index][$i+1] - $lambda_array[$index][$i]);
       $x = $inten_array[$index][$i] +
                $frac*($inten_array[$index][$i+1] - $inten_array[$index][$i]);
       if ($debug > 0) {
         printf "  %8.1f %8.1f %8.1f   %9.4e %9.4e %9.4e \n",
               $lambda_array[$index][$i],
               $fid_lambda,
               $lambda_array[$index][$i+1],
               $inten_array[$index][$i],
               $x,
               $inten_array[$index][$i+1];
       }
       return($x);
     }

   }


   # should never get here
   return(-1);
}
#
#               
