#!/usr/bin/perl
#
# Create simulated datasets which can be used in lab exercises.
#    Idea is to allow students to convert intensities into magnitudes,
#    maybe add a zero point,
#    create light curves,
#    and so forth.
#
# MWR 2/7/2022

use POSIX;
$debug = 1;

$MYPI = 3.14159265;
srand();


#######################
# Create a datafile with INTENSITY measurements of three stars 
#    over a stretch of several hours on a single night
#    No variations as a function of time
#    Write output to a text file
#

$outfile = "./intensity_a.dat";

$JDBASE = 2459610.0;

$inten_a = 1000.0;
$inten_b = $inten_a*0.55;
$inten_c = $inten_a*0.16;

$jd_start = $JDBASE + 0.0;
$jd_end = $JDBASE + 0.40;
$d_jd = (2.0 / 1440.0);


open(OUTFILE, ">$outfile") || die("can't open file $outfile");
printf OUTFILE "# measurements of the intensity of three stars \n";
printf OUTFILE "#%12s  %10s  %10s  %10s \n", "JulianDay  ", "A", "B", "C";
printf OUTFILE "# \n";
for ($jd = $jd_start; $jd <= $jd_end; $jd += $d_jd) {

  # compute an intensity for each star, adding random Poissonian noise
  $obs_a = (1.0 + (gaussian_rand()*(1.0/sqrt($inten_a))))*$inten_a;
  $obs_b = (1.0 + (gaussian_rand()*(1.0/sqrt($inten_b))))*$inten_b;
  $obs_c = (1.0 + (gaussian_rand()*(1.0/sqrt($inten_c))))*$inten_c;

  printf OUTFILE "%12.5f  %10.2f  %10.2f  %10.2f  \n", 
       $jd, $obs_a, $obs_b, $obs_c;

}
close(OUTFILE);



#######################
# Create a datafile with INTENSITY measurements of three stars 
#    over a stretch of several hours on a single night
#    All stars follow a gradual increase, then decrease, in brightness
#    Write output to a text file
#

$outfile = "./intensity_b.dat";

$JDBASE = 2459598.2;

$inten_a = 3300.0;
$inten_b = $inten_a*0.77;
$inten_c = $inten_a*0.32;

# throughput/clouds/etc. vary on this period (days)
$cloud_p = 2.5;

$jd_start = $JDBASE + 0.0;
$jd_end = $JDBASE + 0.40;
$d_jd = (2.0 / 1440.0);


open(OUTFILE, ">$outfile") || die("can't open file $outfile");
printf OUTFILE "# measurements of the intensity of three stars \n";
printf OUTFILE "#%12s  %10s  %10s  %10s \n", "JulianDay  ", "A", "B", "C";
printf OUTFILE "# \n";
for ($jd = $jd_start; $jd <= $jd_end; $jd += $d_jd) {

  # compute a multiplicative factor by which all stars are scaled
  $factor = cos(2.0*$MYPI* ($jd-$JDBASE-2.8) /$cloud_p);

  # compute an intensity for each star, adding random Poissonian noise
  $obs_a = $factor*(1.0 + (gaussian_rand()*(1.0/sqrt($inten_a))))*$inten_a;
  $obs_b = $factor*(1.0 + (gaussian_rand()*(1.0/sqrt($inten_b))))*$inten_b;
  $obs_c = $factor*(1.0 + (gaussian_rand()*(1.0/sqrt($inten_c))))*$inten_c;

  printf OUTFILE "%12.5f  %10.2f  %10.2f  %10.2f  \n", 
       $jd, $obs_a, $obs_b, $obs_c;

}
close(OUTFILE);






#######################
# Create a datafile with INTENSITY measurements of three stars 
#    over a stretch of several hours on a single night
#    A cloud briefly dims all stars equally, then disappears
#    Write output to a text file
#

$outfile = "./intensity_c.dat";

$JDBASE = 2459598.2;

$inten_a = 3300.0;
$inten_b = $inten_a*0.77;
$inten_c = $inten_a*0.32;


$jd_start = $JDBASE + 0.0;
$jd_end = $JDBASE + 0.40;
$d_jd = (2.0 / 1440.0);


# cloud parameters
$cloud_start = ($jd_start + 0.3*($jd_end - $jd_start));
$cloud_mid = ($jd_start + 0.40*($jd_end - $jd_start));
$cloud_end = ($jd_start + 0.45*($jd_end - $jd_start));
# maximum fractional cloud decrease in brightness
$cloud_dip = 0.40;

open(OUTFILE, ">$outfile") || die("can't open file $outfile");
printf OUTFILE "# measurements of the intensity of three stars \n";
printf OUTFILE "#%12s  %10s  %10s  %10s \n", "JulianDay  ", "A", "B", "C";
printf OUTFILE "# \n";
for ($jd = $jd_start; $jd <= $jd_end; $jd += $d_jd) {

  # compute a multiplicative factor by which all stars are scaled
  $factor = 1.0;
  if (($jd > $cloud_start) && ($jd <= $cloud_mid)) {
    $x = ($jd - $cloud_start)/($cloud_mid - $cloud_start);
    $factor = 1.0 - ($cloud_dip*$x);
  }
  if (($jd > $cloud_mid) && ($jd <= $cloud_end)) {
    $x = ($cloud_mid - $jd)/($cloud_end - $cloud_mid);
    $factor = (1.0 - $cloud_dip) - ($cloud_dip*$x);
  }

  # compute an intensity for each star, adding random Poissonian noise
  $obs_a = $factor*(1.0 + (gaussian_rand()*(1.0/sqrt($inten_a))))*$inten_a;
  $obs_b = $factor*(1.0 + (gaussian_rand()*(1.0/sqrt($inten_b))))*$inten_b;
  $obs_c = $factor*(1.0 + (gaussian_rand()*(1.0/sqrt($inten_c))))*$inten_c;

  printf OUTFILE "%12.5f  %10.2f  %10.2f  %10.2f  \n", 
       $jd, $obs_a, $obs_b, $obs_c;

}
close(OUTFILE);




exit 0;






######################################################################
# PROCEDURE: gaussian_rand
#
# DESCRIPTION: Returns a random number drawn from a gaussian (normal)
#              distribution with mean 0.0 and standard deviation 1.0.
#
#              Stolen from 
#
#                 http://www.unix.com.ua/orelly/perl/cookbook/ch02_11.htm
#
#
sub gaussian_rand {
    my ($u1, $u2);  # uniformly distributed random numbers
    my $w;          # variance, then a weight
    my ($g1, $g2);  # gaussian-distributed numbers

    do {
        $u1 = 2 * rand() - 1;
        $u2 = 2 * rand() - 1;
        $w = $u1*$u1 + $u2*$u2;
    } while ( $w >= 1 );

    $w = sqrt( (-2 * log($w))  / $w );
    $g2 = $u1 * $w;
    $g1 = $u2 * $w;
    # could return both if wanted
    #   return wantarray ? ($g1, $g2) : $g1;
    # but I need just one
    return($g1);
}