#!/usr/bin/perl
#
# Create a copy of a stellar profile, and add it to the original
#
# MWR 10/26/2016

$debug = 0;

$num_dat = 0;
$infile = "star_prof.dat";
open(INFILE, "$infile") || die("can't open file $infile");
while (<INFILE>) {
  $line = $_;
  if ($line =~ /^#/) {
    next;
  }
  @words = split(/\s+/, " " . $line);
  $rad[$num_dat] = $words[1];
  $val[$num_dat] = $words[2];
  if ($debug > 0) {
    printf " read item %5d   rad %8.2f  val %8.2f \n",
         $num_dat, $rad[$num_dat], $val[$num_dat];
  }

  $num_dat++;
}
close(INFILE);
if ($debug > 0) {
  printf "# read %6d points from file $infile \n", $num_dat;
}

$xpos = 3.0;
$div = 100.0;
$floor = 3.0;
$gauss_amp = 100.0;
$gauss_fwhm = 4.3;

if (0 == 1) {
  # make a copy of the data, adding a constant to position,
  #    and dividing the data values by a second constant

		  for ($i = 0; $i < $num_dat; $i++) {
			 if ($debug > 0) {
				printf "  rad [ %5d ] is %8.2f   val [ %5d ] is %9.2f \n",
					 $i, $rad[$i], $i, $val[$i];
			 }
			
			 if ($i % 2 == 1) {
				$planet_x[$i] = $xpos + $rad[$i];
			 }
			 else {
				$planet_x[$i] = $xpos - $rad[$i];
			 }
			
			 $planet_val[$i] = ($val[$i] + 0.0) / $div;

			 if ($planet_val[$i] > $floor) {
				printf " %8.2f  %9.2f \n", $planet_x[$i], $val[$i] + $planet_val[$i];
			 }

        }
}
else {

  # add a small gaussian to the data

        $gauss_sigma = $gauss_fwhm / 2.35;

		  for ($i = 0; $i < $num_dat; $i++) {
			 if ($debug > 0) {
				printf "  rad [ %5d ] is %8.2f   val [ %5d ] is %9.2f \n",
					 $i, $rad[$i], $i, $val[$i];
			 }

          $xx = $rad[$i] - $xpos;
          $amp = $gauss_amp * 
                exp( 0.0 - ($xx*$xx)/(2.0*$gauss_sigma*$gauss_sigma));
			
		    $planet_x[$i] = $rad[$i];
			 $planet_val[$i] = $amp;

				printf " %8.2f  %9.2f  %9.2f \n", 
                $planet_x[$i], $val[$i] + $planet_val[$i], $val[$i];

        }
}




exit 0;

