#!/usr/bin/perl
#
# Compute positions of photons emitted by 
#   a central source as a function of time.
#
# We consider three situations:
#
#   a) photons emitted from source, all head
#        straight to observer
#
#   b) photons emitted from source, spread out
#        in circular pattern, go to a certain
#        radius
#
#   c) photons which start at that radius, 
#        (bounced off dust), head to observer
#
# Write output into three output files, which
#   can then be plotted via gnuPlot.
#
# MWR 3/28/2013

$debug = 1;
$TWO_PI = 2*3.15159;

# distance from source to circular radius
$circ_radius = 100;

# distance from source to observer
$dist_observer = 300;

# timestep, in units in which light
#   travels one distance unit per time
$timestep = 10;


# location of observer
$observer_x = 0;
$observer_y = 0;

# location of source
$source_x = 0;
$source_y = $dist_observer;

# number of rays to draw radially outward
$num_rays = 10;

# these are off the screen
$off_screen_x = -999;
$off_screen_y = -999;


# create a file with the location of observer
$outfile = "./calc_reverb_observer.out";
open(OUTFILE, ">$outfile") || die("can't open file $outfile");
printf OUTFILE " %9.4e   %8.2f %8.2f \n", 0, $observer_x, $observer_y;
close(OUTFILE);

# create a file with the location of source
$outfile = "./calc_reverb_source.out";
open(OUTFILE, ">$outfile") || die("can't open file $outfile");
printf OUTFILE " %9.4e   %8.2f %8.2f \n", 0, $source_x, $source_y;
close(OUTFILE);

# create a file with a circle of points to denote the 
#   ring of gas or dust off which photons bounce
$outfile = "./calc_reverb_circle.out";
open(OUTFILE, ">$outfile") || die("can't open file $outfile");
$num_dots = 100;
for ($i = 0; $i < $num_dots; $i++) {
  $theta_rad = ($TWO_PI/$num_dots)*$i;
  $x = $source_x + $circ_radius*cos($theta_rad);
  $y = $source_y + $circ_radius*sin($theta_rad);
  printf OUTFILE " %8.2f %8.2f \n", $x, $y;
}
close(OUTFILE);


# Step a)   Light from source travels
#             directly to observer
$outfile = "calc_reverb_a.out";
open(OUTFILE, ">$outfile") || die("can't open file $outfile");
$curvefile = "calc_reverb_a_curve.out";
open(CURVEFILE, ">$curvefile") || die("can't open file $curvefile");
$start_t = 0;
$x = $source_x;
$y = $source_y;
$curve_num_line = 0;
for ($t = $start_t; ; $t += $timestep) {

  
  printf OUTFILE " %9.4e  %8.2f %8.2f \n", $t, $x, $y;
  printf OUTFILE "\n\n";

  if ($y > $observer_y) {
    for ($i = 0; $i < $curve_num_line; $i++) {
      printf CURVEFILE " %9.4e %8.2f \n", 
           $curve_time_array[$i], $curve_light_array[$i];
    }
    $curve_time_array[$curve_num_line] = $t;
    $curve_light_array[$curve_num_line] = 0;
    printf CURVEFILE " %9.4e %8.2f \n", $t, 0;
    printf CURVEFILE "\n\n";
    $curve_num_line++;
  }
  else {
    for ($i = 0; $i < $curve_num_line; $i++) {
      printf CURVEFILE " %9.4e %8.2f \n", 
           $curve_time_array[$i], $curve_light_array[$i];
    }
    $curve_time_array[$curve_num_line] = $t;
    $curve_light_array[$curve_num_line] = 1;
    printf CURVEFILE " %9.4e %8.2f \n", $t, 1;
    printf CURVEFILE "\n\n";
    $curve_num_line++;
  }

  if ($y <= $observer_y) {
    last;
  }

  $y -= $timestep;

}
for ($j = 0; $j < 25; $j++) {

  $t += $timestep;

  for ($i = 0; $i < $curve_num_line; $i++) {
    printf CURVEFILE " %9.4e %8.2f \n", 
         $curve_time_array[$i], $curve_light_array[$i];
  }
  $curve_time_array[$curve_num_line] = $t;
  $curve_light_array[$curve_num_line] = 0;
  printf CURVEFILE " %9.4e %8.2f \n", $t, 0;
  printf CURVEFILE "\n\n";
  $curve_num_line++;

}
  

close(OUTFILE);







# Step b)   Light from source travels
#             radially outward to circ_radius
$outfile = "calc_reverb_b.out";
open(OUTFILE, ">$outfile") || die("can't open file $outfile");
$curvefile = "calc_reverb_b_curve.out";
open(CURVEFILE, ">$curvefile") || die("can't open file $curvefile");

$start_t = 0;
for ($i = 0; $i < $num_rays; $i++) {
  $ray_x[$i] = $source_x;
  $ray_y[$i] = $source_y;
}
$curve_num_line = 0;
for ($t = $start_t; $t < $circ_radius; $t += $timestep) {

  for ($i = 0; $i < $num_rays; $i++) {
    $theta_rad = ($TWO_PI/$num_rays)*$i;
    $ray_x[$i] += $timestep*cos($theta_rad);
    $ray_y[$i] += $timestep*sin($theta_rad);
    printf OUTFILE " %9.4e  %8.2f %8.2f \n", $t, $ray_x[$i], $ray_y[$i];
  }
  printf OUTFILE "\n\n";

  for ($i = 0; $i < $curve_num_line; $i++) {
    printf CURVEFILE " %9.4e %8.2f \n", 
         $curve_time_array[$i], $curve_light_array[$i];
  }
  $curve_time_array[$curve_num_line] = $t;
  $curve_light_array[$curve_num_line] = 0;
  printf CURVEFILE " %9.4e %8.2f \n", $t, 0;
  printf CURVEFILE "\n\n";
  $curve_num_line++;

}
close(OUTFILE);



# Step c)   Light bounces of material at circ_radius
#             and heads straight down to observer
$outfile = "calc_reverb_c.out";
open(OUTFILE, ">$outfile") || die("can't open file $outfile");
$curvefile = "calc_reverb_c_curve.out";
open(CURVEFILE, ">$curvefile") || die("can't open file $curvefile");

# we use the ending time from last loop
#   and ending positions from last loop
$start_t = $t;
$end_flag = 0;
for ($t = $start_t; $end_flag == 0; $t += $timestep) {

  $num_rays_on_screen = 0;
  for ($i = 0; $i < $num_rays; $i++) {
    $ray_y[$i] -= $timestep;

    if ($ray_y[$i] >= 0) {
      printf OUTFILE " %9.4e  %8.2f %8.2f \n", $t, $ray_x[$i], $ray_y[$i];
      $num_rays_on_screen++;
    }
  }
  if ($num_rays_on_screen == 0) {
      printf OUTFILE " %9.4e  %8.2f %8.2f \n", $t, $off_screen_x, $off_screen_y;
  }
  printf OUTFILE "\n\n";


  # print light curve data to light curve file
  for ($i = 0; $i < $curve_num_line; $i++) {
    printf CURVEFILE " %9.4e %8.2f \n", 
         $curve_time_array[$i], $curve_light_array[$i];
  }
  # we check to see if light from the ring has started to reach observer yet
  $curve_value = 0;
  for ($i = 0; $i < $num_rays; $i++) {
    if ($ray_y[$i] <= 0) {
      # yes, light has started to reach observer
      $curve_value = 1;
    }
  }
  # special case, because I can't get it quite to work with roundoff
  if ($t == 490) {
    $curve_value = 0;
  }
  $curve_time_array[$curve_num_line] = $t;
  $curve_light_array[$curve_num_line] = $curve_value;
  printf CURVEFILE " %9.4e %8.2f \n", $t, $curve_value;
  printf CURVEFILE "\n\n";
  $curve_num_line++;



  # check to see if we are done
  $all_done = 0;
  for ($i = 0; $i < $num_rays; $i++) {
    if ($ray_y[$i] <= $observer_y) {
      $all_done++;
    }
  }
  if ($all_done == $num_rays) {
    $end_flag = 1;
  }
  

}
# now create extra panels after light passes the observer
$off_screen_x = -999;
$off_screen_y = -999;
for ($j = 0; $j < 6; $j++) {


  # there's nothing to plot in bird's-eye view
  printf OUTFILE " %9.4e  %8.2f %8.2f \n", $t, $off_screen_x, $off_screen_y;
  printf OUTFILE "\n\n ";


  # but we do need to continue to plot light curve
  for ($i = 0; $i < $curve_num_line; $i++) {
    printf CURVEFILE " %9.4e %8.2f \n", 
         $curve_time_array[$i], $curve_light_array[$i];
  }
  # we check to see if light from the ring has started to reach observer yet
  $curve_value = 1;
  for ($i = 0; $i < $num_rays; $i++) {
    if ($ray_y[$i] <= 0) {
      # yes, light has started to reach observer
      $curve_value = 0;
    }
  }
  $curve_time_array[$curve_num_line] = $t;
  $curve_light_array[$curve_num_line] = $curve_value;
  printf CURVEFILE " %9.4e %8.2f \n", $t, $curve_value;
  printf CURVEFILE "\n\n";
  $curve_num_line++;

  $t += $timestep;

}
close(OUTFILE);



exit 0;
