
    #   Read background radiation file from standard input
    #   and output histogram binned by local sidereal time in hours.

    require 'readrad.pl';

    $Datasel = 0;                     # Nonzero to randomly select 1/4 data set

    $Prune = 1;                       # Nonzero to prune outliers
    $PruneMu = 5665;                  # Mean (× 100) of composite distribution
    $PruneSigma = 931;                # Std. dev. (× 100) of composite distribution
    $PruneZ = 1;                      # Standard deviations to prune
    $PruneInverse = 0;                # Keep only outliers if nonzero

    $Longitude = -7;                  # Longitude (- = East, + = West)

    $day = 24 * 60 * 60;
    $Binsize = 60 * 60;               # Bin size in seconds
    $SecPerDay = 24 * 60 * 60;

    require 'csv.pl';

    srand(time());

    while (($counts, $time) = &readrad) {
        if ($time == -1) {
            last;
        }

        #   Random data subset selection
        if ($Datasel) {
            if (int(rand(32768)) & 3 != 2) {
                next;
            }
        }

        #   Outlier pruning
        if ($Prune) {
            if ($PruneInverse ^
                (($counts < ($PruneMu - ($PruneSigma * $PruneZ))) ||
                 ($counts > ($PruneMu + ($PruneSigma * $PruneZ))))) {
#               print(STDERR "Prune $counts\n");
                next;
            }
        }

        #   Calculate local sidereal time from UTC

        ($sec, $min, $hour, $mday, $mon, $year) = gmtime($time);
#printf("%d/%d/%d %2d:%02d", $year + 1900, $mon + 1, $mday, $hour, $min);

        #   Calculate Greenwich Mean Sidereal time and bin by hour

        $julianDate = &utctoj($year + 1900, $mon, $mday, $hour, $min, $sec);
        $gmSt = &gmst($julianDate);
        $stime = int($gmSt * 60 * 60);

        #   Adjust Greenwich mean sidereal time to local meridian
        $stime -= ($Longitude / 15) * (60 * 60);
#printf("  JD = $julianDate GMST = $gmSt stime = %d\n", $stime);

        $btime = int($stime / $Binsize) * $Binsize;
        if ($btime < 0) {
            $btime += $SecPerDay;
        } elsif ($btime >= $SecPerDay) {
            $btime -= $SecPerDay;
        }
#print("Binned time = $btime\n");

        $hour = int(($btime % $day) / (60 * 60));
        $histo[$hour] += $counts;
        $points[$hour]++;
#print("$time $hour\n");
    }

    for ($i = 0; $i < 24; $i++) {
        if (!(defined $histo[$i])) {
            $rad = 0;
            $histo[$i] = 0;
        } else {
            $rad = $histo[$i] / ($points[$i] * 354.0);
        }
        print("$i,$histo[$i],$points[$i],$rad\n");
    }

    #   UTCTOJ  --  Convert GMT date and time to astronomical
    #               Julian time (i.e. Julian date plus day fraction,
    #               expressed as a double).
    #
    #               IMPORTANT:  The month argument ($mon) is a number
    #               between 0 (January) and 11 (December), *not* a
    #               1-based index as usually written in dates.  This
    #               is compatible with the month index returned by
    #               gmtime().

    sub utctoj {
        local($year, $mon, $mday, $hour, $min, $sec) = @_;
        local($a, $b, $m, $y);

        #  Algorithm as given in Meeus, Astronomical Algorithms, Chapter 7, page 61

        die "utctoj: argument out of range" unless ($mon  >= 0 && $mon  < 12);
        die "utctoj: argument out of range" unless ($mday >  0 && $mday < 32);
        die "utctoj: argument out of range" unless ($hour >= 0 && $hour < 24);
        die "utctoj: argument out of range" unless ($min  >= 0 && $min  < 60);
        die "utctoj: argument out of range" unless ($sec  >= 0 && $sec  < 60);

        $m = $mon + 1;
        $y = $year;

        if ($m <= 2) {
            $y--;
            $m += 12;
        }

        #  Determine whether date is in Julian or Gregorian calendar based on
        #  canonical date of calendar reform.

        if (($year < 1582) || (($year == 1582) &&
            (($mon < 9) || ($mon == 9 && $mday < 5)))) {
            $b = 0;
        } else {
            $a = int($y / 100);
            $b = 2 - $a + int($a / 4);
        }

        ((int(365.25 * ($y + 4716))) + (int(30.6001 * ($m + 1))) +
                    $mday + $b - 1524.5) +
                (($sec + 60 * ($min + 60 * $hour)) / 86400.0);
    }

    #  GMST  --  Calculate Greenwich Mean Sidereal Time for a given
    #            instant expressed as a Julian date and fraction.

    sub gmst {
        local ($jd) = @_;
        local ($t, $theta0);

        #  Time, in Julian centuries of 36525 ephemeris days,
        #  measured from the epoch 1900 January 0.5 ET.

        $t = ((int($jd + 0.5) - 0.5) - 2415020.0) / 36525.0;

        $theta0 = 6.6460656 + 2400.051262 * $t + 0.00002581 * $t * $t;

        $t = ($jd + 0.5) - (int($jd + 0.5));

        $theta0 += ($t * 24.0) * 1.002737908;

        ($theta0 - 24.0 * (int($theta0 / 24.0)));
    }
