#!/usr/bin/perl
#----------------------------------------------------------------------------------------------
# by Joseph A. DiVerdi
# Copyright 2006 by XTR Systems, LLC
#----------------------------------------------------------------------------------------------
# Includes and other external modules
	use warnings;
	use strict;
	use Math::Trig;
	use Math::Complex;
	use Time::Local;
	use POSIX qw(ceil floor);
#----------------------------------------------------------------------------------------------
# Main execution module
	$| = 1;
	
	# locations of important items
	my $directions_file_name = "directions to use";
	my $usable_files_file_name = "files to use";
	my $data_files_directory = "data files";
	my $processed_data_file_name = "processed_data_scaled";
	my $gridded_data_file_name = "gridded_data_scaled";
	
	my $dec_spacing = 1.0;
	my $ra_spacing = 2.0;
	
	# variable definitions
	my @file_names;
	my $azimuth_list;
	my $elevation_list;
	
	# create the arrays for gridding the data
	my $list;
	for (my $dec = 0; $ dec < 90; $dec += $dec_spacing) {
		for (my $ra = 0; $ ra < 360; $ra += $ra_spacing) {
			@{$list->{$ra}{$dec}} = ();
		}
	}
		
	# set up the names of the data files to process
	open INPUT, '<', $usable_files_file_name or
		die "Can't read '$usable_files_file_name': $!\n";
	while (my $line = <INPUT>) {
		next if $line =~ /^(#|\s*$)/;
		chomp $line;
		push @file_names, $line;
	}
	close INPUT;
	
	# set up the pointing directions
	open INPUT, '<', $directions_file_name or
		die "Can't read '$directions_file_name': $!\n";
	while (my $line = <INPUT>) {
		next if $line =~ /^(#|\s*$)/;
		chomp $line;
		my ($azimuth, $elevation, $day, $month, $year, $time, $remainder) = split /,?\s+/, $line;
		my $epoch_time = convert_time_and_date($day, $month, $year, $time);
		$azimuth_list->{$epoch_time} = $azimuth;
		$elevation_list->{$epoch_time} = $elevation;
	}
	close INPUT;
	
	print "Processing data files";
	open PROCESSED, '>', $processed_data_file_name or
		die "Can't create '$processed_data_file_name': $!\n";
	print PROCESSED "epoch_time\tdeclination\tright ascension\tintensity\n";
		
	# iterate over the raw data files
	print "Processing data files";
	foreach (@file_names) { print "|"; }
	foreach (@file_names) { print "\x08"; }
	foreach my $input_file_name (sort @file_names) {
		print ".";
		$input_file_name = join("/", $data_files_directory, $input_file_name);
		unless (open INPUT, '<', $input_file_name) {
			print "Can't process '$input_file_name': $!\n";
			next;
		}
		
		# read in the entire contents of a particular data file
		my $data_record;
		my $intensity_min = 10;
		my $intensity_max = 0;
		while (my $line = <INPUT>) {
			next if $line =~ /^(#|\s*$)/;
			chomp $line;
			
			my ($epoch_time, $el, $az, @values) = split /\s+/, $line;
			
			# compute elevation and azimuth for this epoch time
			($el, $az) = compute_horizontal($epoch_time, $elevation_list, $azimuth_list);
			my ($dec, $ra) = compute_equatorial($epoch_time, $el, $az);
			
 			my $integral = 0;
 			foreach my $value (@values) { $integral += $value; }
 			$integral /= scalar @values;
 			
 			# store the values in a suitable structure
 			$data_record->{$epoch_time}{'dec'} = $dec;
 			$data_record->{$epoch_time}{'ra'} = $ra;
 			$data_record->{$epoch_time}{'integral'} = $integral;
 			$data_record->{$epoch_time}{'values'} = [@values];
 			
			$intensity_min = $integral if $integral < $intensity_min;
			$intensity_max = $integral if $integral > $intensity_max;
		}
 			
		foreach my $epoch_time (sort keys %$data_record) {
			my $dec = $data_record->{$epoch_time}{'dec'};
			my $ra = $data_record->{$epoch_time}{'ra'};
#			my $integral = ($data_record->{$epoch_time}{'integral'} - $intensity_min) / ($intensity_max - $intensity_min);
			my $integral = ($data_record->{$epoch_time}{'integral'} - $intensity_min);
# 			my $integral = ($data_record->{$epoch_time}{'integral'} - $intensity_max + 1.0);
			my @values = @{$data_record->{$epoch_time}{'values'}};
			
			printf PROCESSED "%10d\t%+8.3f\t%+8.3f\t%7.5f", $epoch_time, $dec, $ra, $integral;
# 			printf PROCESSED "\t%s", join("\t", @values);
			printf PROCESSED "\n";
			
			# quantize the pointing direction on the specified grid
			# put the data value on the list of values for that pointing direction
			my $dec_quantized = int($dec / $dec_spacing + 0.5) * $dec_spacing;
			my $ra_quantized = int($ra / $ra_spacing + 0.5) * $ra_spacing;
			push @{$list->{$ra_quantized}{$dec_quantized}}, $integral;
		}
		
		close INPUT;
	}
	close PROCESSED;
	print "\n";
		
	# form the mean and (sample) std dev and write the gridded data file
	print "Gridding data";
	open GRIDDED, '>', $gridded_data_file_name or
		die "Can't create '$gridded_data_file_name': $!\n";
	print GRIDDED "declination\tright ascension\tcount\tintensity\n";
	
	foreach (my $dec_quantized = 0; $dec_quantized < 90; $dec_quantized += $dec_spacing) { print "|"; }
	foreach (my $dec_quantized = 0; $dec_quantized < 90; $dec_quantized += $dec_spacing) { print "\x08"; }
	for (my $dec_quantized = 0; $dec_quantized < 90; $dec_quantized += $dec_spacing) {
		print ".";
		for (my $ra_quantized = 0; $ra_quantized < 360; $ra_quantized += $ra_spacing) {
			my @list = sort @{$list->{$ra_quantized}{$dec_quantized}};
			
			my $median;
			if (scalar @list == 0) {
				$median = 0;
			} elsif ((scalar @list) % 2 == 0) {
				$median = ($list[floor($#list / 2)] + $list[ceil($#list / 2)]) / 2;
			} else {
				$median = $list[$#list / 2];
			}
			
			printf GRIDDED "%+8.3f\t%+8.3f\t%05d\t%7.5f", $dec_quantized, $ra_quantized, (scalar @list), $median;
			if (scalar @list > 0 and $list[0] != 0 and $list[$#list] != 0) {
	 			foreach my $value (@list) {
	 				printf GRIDDED "\t%7.5f", $value;
	 			}
			}
			printf GRIDDED "\n";
		}
	}
	
	close GRIDDED;
	print "\n";
	
	exit;
#----------------------------------------------------------------------------------------------
sub convert_time_and_date {
	my $day = shift;
	my $month = shift;
	my $year = shift;
	my $hour_minute = shift;

	if ($month =~ /jan/i) { $month = 0; }
	elsif ($month =~ /feb/i) { $month = 1; }
	elsif ($month =~ /mar/i) { $month = 2; }
	elsif ($month =~ /apr/i) { $month = 3; }
	elsif ($month =~ /may/i) { $month = 4; }
	elsif ($month =~ /jun/i) { $month = 5; }
	elsif ($month =~ /jul/i) { $month = 6; }
	elsif ($month =~ /aug/i) { $month = 7; }
	elsif ($month =~ /sep/i) { $month = 8; }
	elsif ($month =~ /oct/i) { $month = 9; }
	elsif ($month =~ /nov/i) { $month = 10; }
	elsif ($month =~ /dec/i) { $month = 11; }
	else { die "Bad month: '$month'\n"; }
	
	my ($hour, $minute) = $hour_minute =~ /(\d{2})(\d{2})/;
	
	return timelocal(0, $minute, $hour, $day, $month, $year);
}
#----------------------------------------------------------------------------------------------
sub compute_horizontal {
	my $epoch_time = shift;
	my $elevation_list = shift;
	my $azimuth_list = shift;
	
	my $time;
	foreach my $test_time (reverse sort keys %{$elevation_list}) {
		 if ($test_time <= $epoch_time) { $time = $test_time; last; }
	}
	
	return ($elevation_list->{$time}, $azimuth_list->{$time});
}
#----------------------------------------------------------------------------------------------
sub compute_equatorial {
	my $time = shift;
	my $el = shift;
	my $az = shift;
	
	my $lat = +40.14791; # coordinates for Table Mountain Upper Dish
	my $lon = -105.2324; # coordinates for Table Mountain Upper Dish
	
	# The following formulas used for computing sidereal time in this program have been published by 
	# Keith Burnett (kburnett@btinternet.com) on his web page: http://www.btinternet.com/~kburnett/kepler/sidereal.htm
	# And related material at: http://www.xylem.f2s.com/kepler/mean.html
	
	# All the formulas connecting UT with sidereal time depend on the time elapsed from some fundamental epoch. 
	# Older references have used a date in 1900 as the epoch, Meeus and most recent formulas use J2000.0 as the epoch. 
	# The formula below will find the days since J2000.0 for any date after 1900 and before 2099.
	# The date and time must be in UT.
	
	my @times = gmtime($time);
	my $year = $times[5] + 1900;	# four digit year
	my $month = $times[4] + 1;		# 1 <= month <= 12
	my $day = $times[3];					# 1 <= day <= 31
	
	my $epoch_days = 367 * $year - int(7 * ($year + int(($month + 9) / 12)) / 4) + int(275 * $month / 9) + $times[3] - 730531.5;
	$epoch_days += ($times[2] + $times[1] / 60 + $times[0] / 3600) / 24;
	
	# The mean sidereal time at zero longitude is often called
	# Greenwich Mean Sidereal Time or GMST. The formula below is based on
	# Meeus formula 11.4 with terms in the square and the cube of the time left out. 
	# epoch_days is the number of UT days since J2000.0, including fractional part of a day and GMST is in degrees.
	# The value may require rotation to get it in the range 0 <= GMST < 360
	my $GMST = rotate_into_range(280.46061837 + 360.98564736629 * $epoch_days);
	
	# To get the sidereal time at the local longitude, known as Local Mean
	# Sidereal Time or LMST, add the local longitude in degrees, taking East as Positive.
	# The value may require rotation to get it in the range 0 <= LMST < 360
	my $LMST = rotate_into_range($GMST + $lon);
	
	# calculate the declination
	my $dec = asin_deg(sin_deg($el) * sin_deg($lat) + cos_deg($el) * cos_deg($lat) * cos_deg($az));
	
	# compute the hour angle
	# take the real part of the potentially complex value (which can occur near singularities) and convert it to degrees
	my $hour_angle = rad2deg(Re(acos((sin_deg($el) - sin_deg($lat) * sin_deg($dec)) / (cos_deg($lat) * cos_deg($dec)))));
	
	# reflect the hour angle depending on which quadrant the azimuth is in then rotate into the principal range
	$hour_angle = 360 - $hour_angle if 0 <= $az and $az < 180;
	$hour_angle = rotate_into_range($hour_angle);
	
	# The RA is the difference between the LMST and the hour angle rotated into the principal range
	my $ra = rotate_into_range($LMST - $hour_angle);
	
	return ($dec, $ra);
}
#----------------------------------------------------------------------------------------------
sub rotate_into_range { 	# rotate the argument into the range 0 <= argument < 360
	my $angle = shift;
	while ($angle >= 360) { $angle -= 360; }
	while ($angle <    0) { $angle += 360; }
	return $angle;
}
#----------------------------------------------------------------------------------------------
sub sin_deg  { return  sin(deg2rad(shift)); }
sub cos_deg  { return  cos(deg2rad(shift)); }
sub asin_deg { return rad2deg(asin(shift)); }
sub acos_deg { return rad2deg(acos(shift)); }
#----------------------------------------------------------------------------------------------