#!/usr/bin/perl

my $VERSION = '$Revision: 1.5 $' ;

#-- ----------------------------------------------------------------------------
##-- From http://egb13.net/
##-- ----------------------------------------------------------------------------
##-- Copyright (C) 2011 by EGB13.net
##-- 
##-- Permission is hereby granted, free of charge, to any person obtaining a copy
##-- of this software and associated documentation files (the "Software"), to deal
##-- in the Software without restriction, including without limitation the rights
##-- to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
##-- copies of the Software, and to permit persons to whom the Software is
##-- furnished to do so, subject to the following conditions:
##-- 
##-- The above copyright notice and this permission notice shall be included in
##-- all copies or substantial portions of the Software.
##-- 
##-- THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
##-- IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
##-- FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
##-- AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
##-- LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
##-- OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
##-- THE SOFTWARE.
##-- ----------------------------------------------------------------------------

$scale = 3.0 ;

use GD;
use POSIX qw(mktime) ;
use Math::Trig ;


#-- Date widget initializations ---------------------------------------
my @moname = ( 'J', 'F', 'M', 'A', 'M', 'J', 'J', 'A', 'S', 'O', 'N', 'D') ;
my @mopos ;
my $morad = 45 ;
my $dwx = 55 ;
my $dwy = 355 ;
foreach my $mo (0 .. 11) {
	my $deg = ( ( (12-$mo)/12.0) * 360 ) ;
	my $x = ($dwx-3) + sin(deg2rad($deg)) * $morad ;
	my $y = ($dwy-6) - cos(deg2rad($deg+180)) * $morad ;
	$mopos[$mo] = [$x, $y] ;
}
my $degperday = 360.0 / 365.0 ;
my $dps = 3 ; #degrees per slice
#----------------------------------------------------------------------


# my @mo ;
# my @tn ;
# &ReadState ("mo4.txt", \@mo) ;
# &ReadState ("tn2.txt", \@tn) ;

my @states ;
&ReadState ("states.txt", \@states) ;
my @cities ;
&ReadCities ("cities.txt", \@cities) ;



##East North DATE-O.T.-(UTC) DEP MAG
#833.81 4030.31 1974/06/29-09:27:09.5000 5.0000 1.6
#813.61 4017.33 1974/07/04-13:58:06.3000 5.0000 1.8

my $dstart = mktime(0, 0, 0, 29, 5, 74, 0,0,-1) ;
my $daysec = 24 * 60 * 60 ;
my $dend = $dstart + $daysec ;
my @day ;
my @qlist ;
my $dix = 30 ;

my ($min_x, $min_y, $max_x, $max_y) = (9999999, 9999999, -9999999, -9999999) ;

#-
#- Group quake locations & magnitudes by day
#-
open (Q, "qcat") || die ("qcat: $!\n") ;
while (<Q>) {
	next if (m/^\s*\D/) ;
	chomp ;
	my ($east, $north, $tstring, $depth, $mag) = split (/\s+/) ;
	my ($year, $mon, $mday, $hour, $min, $sec) = ($tstring =~ m=(\d\d\d\d)/(\d\d)/(\d\d)-(\d\d):(\d\d):(\d\d)=) ;
	my $time_t = mktime($sec, $min, $hour, $mday, $mon, ($year-1900), 0,0,-1) ;
	my  @gmt = gmtime($time_t);
	while ($time_t >= $dend) {
		$dstart = $dend ;
		$dend = $dstart + $daysec ;
		$dix++ ;
	}
	
	push @{$day[$dix]}, [$east, $north, $mag, $gmt[5]+1900, $gmt[7] ] ;	# gmt year & yday
	push @qlist, [$east, $north, $mag] ; # ungrouped for summary slide only
	$min_x = $east if ($east < $min_x) ;
	$min_y = $north if ($north < $min_y) ;
	$max_x = $east if ($east > $max_x) ;
	$max_y = $north if ($north > $max_y) ;
}

#- example access to the data structure we just built
if (defined ($opt_d) ) {
	foreach my $ix (0 .. $#day) {
		printf "%d: %d\n"
		,	$ix
		,	$#{$day[$ix]}
		;
		foreach my $jx (0 .. $#{$day[$ix]}) {
			printf "\t%0.1f %0.1f   %0.1f\n"
			,	${$day[$ix]}[$jx]->[0]
			,	${$day[$ix]}[$jx]->[1]
			,	${$day[$ix]}[$jx]->[2]
			;
		}
	}
}


#- scaling for map
my ($mwide, $mhigh) = (640, 480) ;
my $xscale = ( $mwide / ( $max_x - $min_x ) ) ;
my $yscale = ( $mhigh / ( $max_y - $min_y ) ) ;
$xscale = $yscale if ($yscale < $xscale) ;
$yscale = $xscale if ($xscale < $yscale) ;
my $xoffs  = 0 - $min_x  + ( ( ($mwide / $xscale) - ($max_x - $min_x) ) / 2) ;
my $yoffs  = 0 - $min_y ;


printf "Bounding box: %0.3f %0.3f  ::  %0.3f %0.3f\n"
	, $min_x
	, $min_y
	, $max_x
	, $max_y
	;


#- for each day print the next month's worth

my $year ;
my $yday ;
my $lnd ;  #last known day
my $lndx ; #last known day index

foreach my $dx (0 .. $#day) {
	if (${$day[$dx]}[0]->[3] ne "") {
		$year = ${$day[$dx]}[0]->[3] ;
		$yday = ${$day[$dx]}[0]->[4] ;
		$lndx = $dx ;
		$lnd = $yday ;
		last ;
	}
}

#-------------------------------------------------------------------------------
#- title slide
#-
my $map = new GD::Image ($mwide, $mhigh, 1);
$map->trueColor(1) ;
$map->alphaBlending(1) ;
my $white = $map->colorAllocateAlpha (255,255,255, 0);
my $black = $map->colorAllocateAlpha (0,0,128, 0) ;
my $grey = $map->colorAllocateAlpha (96,96,96, 0) ;
my $bgcolor = $map->colorAllocate(0xd8,0xd8,0xd8) ;

$map->filledRectangle (0,0, $mwide,$mhigh, $bgcolor) ;
$map->filledRectangle (
		  int ( ( ($min_x + $xoffs) * $xscale) + 0.5)
		, int ( ( ($min_y + $yoffs) * $xscale) + 0.5)
		, int ( ( ($max_x + $xoffs) * $xscale) + 0.5)
		, int ( ( ($max_y + $yoffs) * $xscale) + 0.5)
		, $white
		) ;
&DrawState (\@states, $map, $black) ;
&DrawCities (\@cities, $map, $grey) ;
&DrawMagScale ($map, $black) ;
&DrawTitle ($map, $black) ;
&DrawLegendScale ($map, $black) ;

my $fname = "frames/f-0a.png" ;
open (PNG, ">".$fname) || die ("test.png: $!\n") ;
print PNG $map->png ;
close PNG ;
undef $map ;

#--

my $map = new GD::Image ($mwide, $mhigh, 1);
$map->trueColor(1) ;
$map->alphaBlending(1) ;
my $white = $map->colorAllocateAlpha (255,255,255, 0);
my $black = $map->colorAllocateAlpha (0,0,128, 0) ;
my $grey = $map->colorAllocateAlpha (96,96,96, 0) ;
my $bgcolor = $map->colorAllocate(0xd8,0xd8,0xd8) ;

$map->filledRectangle (0,0, $mwide,$mhigh, $bgcolor) ;
$map->filledRectangle (
		  int ( ( ($min_x + $xoffs) * $xscale) + 0.5)
		, int ( ( ($min_y + $yoffs) * $xscale) + 0.5)
		, int ( ( ($max_x + $xoffs) * $xscale) + 0.5)
		, int ( ( ($max_y + $yoffs) * $xscale) + 0.5)
		, $white
		) ;
&DrawState (\@states, $map, $black) ;
&DrawCities (\@cities, $map, $grey) ;
&DrawMagScale ($map, $black) ;
&DrawTitle ($map, $black) ;
&DrawLegendScale ($map, $black) ;

my @qq = sort { $a->[2] <=> $b->[2]; } @qlist ;
foreach my $qref (@qq) {
	my $x = int ( ( ($qref->[0] + $xoffs) * $xscale) + 0.5) ;
	my $y = $mhigh - int ( ( ($qref->[1] + $yoffs) * $xscale) + 0.5) ;
	my $r = $scale * $qref->[2] ;
	my $fade = (5.0-$qref->[2]) / 5.0 ;
	$fade = 1 if ($fade > 1) ;

	$r = 0.5 if ($r < 0.5) ;
	my $color = $map->colorAllocateAlpha(
		  128+(127*$fade)
		, 255*$fade
		, 255*$fade
		, 0
		) ;
	$map->filledArc ($x,$y, $r,$r, 0,360, $color);
}

my $fname = "frames/f-0b.png" ;
open (PNG, ">".$fname) || die ("test.png: $!\n") ;
print PNG $map->png ;
close PNG ;
undef $map ;

exit 0 ;

#-------------------------------------------------------------------------------
#- Day-by-day frames
#-

foreach my $dx (0 .. $#day) {

	printf "[ DAY %d of %d] %0.1f%% \n"
	,	$dx
	,	$#day
	,	($dx / $#day) * 100.0
	if ( ($dx % 10) == 0) ;

	my $map = new GD::Image ($mwide, $mhigh, 1);
	$map->trueColor(1) ;
	$map->alphaBlending(1) ;
	my $white = $map->colorAllocateAlpha (255,255,255, 0);
	my $black = $map->colorAllocateAlpha (0,0,128, 0) ;
	my $grey = $map->colorAllocateAlpha (96,96,96, 0) ;
	my $bgcolor = $map->colorAllocate(0xd8,0xd8,0xd8) ;
	my $fade ;

	$map->filledRectangle (0,0, $mwide,$mhigh, $bgcolor) ;
	$map->filledRectangle (
			  int ( ( ($min_x + $xoffs) * $xscale) + 0.5)
			, int ( ( ($min_y + $yoffs) * $xscale) + 0.5)
			, int ( ( ($max_x + $xoffs) * $xscale) + 0.5)
			, int ( ( ($max_y + $yoffs) * $xscale) + 0.5)
			, $white
			) ;

	foreach my $tomorrow (0 .. 29) {
		my $ix = $dx + $tomorrow ;
		foreach my $jx (0 .. $#{$day[$ix]}) {
			#$fade = int ( ( (30.0 - $tomorrow) / 30.0) * 127) ;
			$fade = ( (30.0 - $tomorrow) / 30.0) ;
			if (defined ($opt_d) ) {
				printf "%d\t%0.1f %0.1f   %0.1f (%0.2f)\n"
				,	$tomorrow
				,	${$day[$ix]}[$jx]->[0]
				,	${$day[$ix]}[$jx]->[1]
				,	${$day[$ix]}[$jx]->[2]
				,	$fade
				;
			}

			my $x = int ( ( (${$day[$ix]}[$jx]->[0] + $xoffs) * $xscale) + 0.5) ;
			my $y = $mhigh - int ( ( (${$day[$ix]}[$jx]->[1] + $yoffs) * $xscale) + 0.5) ;
			##my $r = int ( (2.0 * ${$day[$ix]}[$jx]->[2]) + 0.5) ;
			my $r = $scale * ${$day[$ix]}[$jx]->[2] ;
			$r = 0.5 if ($r < 0.5) ;
			if (defined ($opt_d) ) {
				printf "	+%d+%d %d (%d)\n"
					,	$x
					,	$y
					,	$r
					,	$fade
					;
			}
			my $color = $map->colorAllocateAlpha(
				  128+(127*$fade)
				, 255*$fade
				, 255*$fade
				, int ($fade * 127)
				) ;
			$map->filledArc ($x,$y, $r,$r, 0,360, $color);
		}

		if ($tomorrow == 29) {
			if (${$day[$dx+$tomorrow]}[$jx]->[3] ne "") {
				$lndx = $dx+$tomorrow ;
				$year = ${$day[$lndx]}[$jx]->[3] ;
				$yday = ${$day[$lndx]}[$jx]->[4] ;
				$lnd = $yday ;
			} else {
				$yday = $lnd + $tomorrow + $dx - $lndx ;
				$year++ if ($yday == 366) ;
			}
			&DateWidget ($map, $year, $yday, $bgcolor, $white, $black) ;
		}

	}
	if ($dx >= 30) {
		&DrawState (\@states, $map, $black) ;
		&DrawCities (\@cities, $map, $grey) ;
		&DrawMagScale ($map, $black) ;
		&DrawTitle ($map, $black) ;
		&DrawLegendScale ($map, $black) ;

		my $fname = sprintf "frames/f-%d.png", $dx-29 ;
		open (PNG, ">".$fname) || die ("test.png: $!\n") ;
		print PNG $map->png ;
		close PNG ;
	}
	undef $map ;

	$dx++ ;
}
print "\n" ;



exit 0 ;

################################################################################
################################################################################


sub ReadState {	my ($stfile, $stref) =@_ ;
	#1	841580 4079370 0
	#2	842050 4079000 0

	open (ST, $stfile) || die ("$stfile: $!\n") ;
	my $mx = 0 ;
	while (<ST>) {
		if (m/^\s*$/) {
			$mx++ ;
			next ;
		}
		my ($x, $y) = m/\s*(\d+\.?\d*)\s(\d+\.?\d*)\b/ ;
		push @{$$stref[$mx]}, [$x, $y] ;
	}
	close ST ;

}
#----------------------------------------------------------------------

sub DrawState {	my ($stref, $mapref, $coloref) = @_ ;

	foreach my $seg (0 .. $#{$stref}) {
		foreach my $pt (1 .. $#{$$stref[$seg]}) {

			my $x1 = int ( ( ($$stref[$seg]->[$pt-1]->[0] + $xoffs) * $xscale) + 0.5) ;
			my $y1 = $mhigh - int ( ( ($$stref[$seg]->[$pt-1]->[1] + $yoffs) * $xscale) + 0.5) ;
			my $x2 = int ( ( ($$stref[$seg]->[$pt]->[0] + $xoffs) * $xscale) + 0.5) ;
			my $y2 = $mhigh - int ( ( ($$stref[$seg]->[$pt]->[1] + $yoffs) * $xscale) + 0.5) ;

#			printf "%d,%d %d,%d  "
#				, $x1
#				, $y1
#				, $x2
#				, $y2
#				;


			$mapref->line	( $x1
					, $y1
					, $x2
					, $y2
					, $coloref
					) ;
		}
#		print "\n" ;
	}
}
#----------------------------------------------------------------------

sub ReadCities {	my ($ctfile, $ctref) =@_ ;

	#724918.837693	3993166.218603	Paragould

	open (CT, $ctfile) || die ("$ctfile: $!\n") ;
	while (<CT>) {
		my ($x, $y, $n) = m/\s*(\d+\.?\d*)\s(\d+\.?\d*)\s+(\S+.*\S)\s*$/ ;
		push @{$ctref}, [$x, $y, $n] ;
	}
	close CT ;

}
#----------------------------------------------------------------------

sub DrawCities { my ($ctref, $mapref, $coloref) = @_ ;

	foreach my $city (0 .. $#{$ctref}) {
		my $x = int ( ( (${$ctref}[$city]->[0] + $xoffs) * $xscale) + 0.5) ;
		my $y = $mhigh - int ( ( (${$ctref}[$city]->[1] + $yoffs) * $xscale) + 0.5) ;
		$mapref->arc ($x,$y, 5,5, 0,360, $coloref);
		$mapref->string(gdMediumBoldFont, $x+5, $y-4, ${$ctref}[$city]->[2], $coloref);
	}
}
		

#----------------------------------------------------------------------


sub DateWidget {	my ($mapref, $year, $day, $bgc, $dial, $line) = @_ ;

	printf "<%d/%d>\n", $year, $day if (defined ($opt_d) ) ;

	$day -= 30 ;
	if ($day < 0) {
		$day += 365 ;
		$year -= 1;
	}

	# create a new image

	# allocate some colors
	# fill 
	$mapref->filledRectangle(0,0,99,99,$bgc);
	$mapref->filledArc($dwx,$dwy, 75,75, 0,360, $dial) ;

	$ddeg = 90 + ($day * $degperday) ;
	my $dbak = 0 ;
	while ($dbak < 30) {
		my $pct = $dbak / 30 ;
		my $color = $mapref->colorAllocate (
				  128 + int(127*$pct)
				, 255 * $pct
				, 255 * $pct
				) ;
		my $start = int ($ddeg - ($dbak * $degperday) ) ;
		my $end = int ($ddeg - ( ($dbak + $dps) * $degperday) ) ;
		printf "%d:%d > \$mapref->filledArc($dwx,$dwy, 75,75, %d,%d,\$color);\n" , $day, $dbak, $end , $start if (defined ($opt_d) ) ;
		$mapref->filledArc($dwx,$dwy, 75,75, $end,$start,$color);
		$dbak += 1 ;
	}

	#- color center
	$mapref->filledArc($dwx,$dwy, 40,40, 0,360, $bgc) ;
	$mapref->arc($dwx,$dwy, 40,40, 0,360, $black) ;
	$mapref->arc($dwx,$dwy, 75,75, 0,360, $black) ;
	#- year
	$mapref->string(gdMediumBoldFont,($dwx-13),($dwy-6),$year,$black);
	#- months
	foreach my $mo (0 .. 11) {
		$mapref->string(gdMediumBoldFont, $mopos[$mo]->[0], $mopos[$mo]->[1], $moname[$mo], $black) ;
	}

}
#----------------------------------------------------------------------

sub DrawMagScale { my ($img, $blk) = @_ ;

	my $x = 550 ;
	my $dy = 25 ;
	my $y0 = 300 ;
	my $color = $img->colorAllocate (128,0,0) ;
	foreach my $mag (1 .. 5) {

			my $r = $scale * $mag ;
			$img->filledArc	( $x
					, $y0 + $dy * ($mag - 1)
					, $r, $r
					, 0, 360
					, $color
					);
			$img->string	( gdMediumBoldFont
					, $x + 12
					, $y0 + $dy * ($mag - 1) - 6
					, "Mag. ".$mag
					, $black
					) ;
	}
}

#----------------------------------------------------------------------

sub DrawTitle {	my ($img, $clr) = @_ ;

	my $tX = 7 ;
	my $tY = 20 ;
	my $dX = 8 ;
	my $dY = 15 ;

	$img->string	( gdGiantFont
			, $tX + $dX
			, $tY
			, "New Madrid"
			, $clr
			) ;
	$img->string	( gdGiantFont
			, $tX
			, $tY + $dY
			, "Seismic Zone"
			, $clr
			) ;

	$img->string	( gdSmallFont
			, $tX
			, $tY + $dY * 3
			, "Data from:"
			, $clr
			) ;

	$img->string	( gdSmallFont
			, $tX + 3
			, $tY + $dY * 4
			, "www.ceri."
			, $clr
			) ;
	$img->string	( gdSmallFont
			, $tX+6
			, $tY + $dY * 5
			, "memphis.edu"
			, $clr
			) ;

}
#----------------------------------------------------------------------

sub DrawLegendScale {	my ($img, $clr) = @_ ;


	#- 277 km (172 mi) wide
	#- 332 km (207 mi) high
	#- 89.406681 37.990597

	my $tX = 530 ;
	my $tY = 20 ;
	my $dX = 8 ;
	my $dY = 15 ;

	$img->string	( gdMediumBoldFont
			, $tX
			, $tY
			, "Map area:"
			, $clr
			) ;

	$img->string	( gdMediumBoldFont
			, $tX
			, $tY + $dY*1
			, " 227 x 332 km"
			, $clr
			) ;
	$img->string	( gdMediumBoldFont
			, $tX
			, $tY + $dY*2
			, " 172 x 207 mi"
			, $clr
			) ;

	$img->string	( gdMediumBoldFont
			, $tX
			, $tY + $dY*4
			, "Center:"
			, $clr
			) ;
	$img->string	( gdMediumBoldFont
			, $tX
			, $tY + $dY*5
			, " 36.495621 N"
			, $clr
			) ;
	$img->string	( gdMediumBoldFont
			, $tX
			, $tY + $dY*6
			, " 89.477055 W"
			, $clr
			) ;

	$img->string	( gdMediumBoldFont
			, $tX
			, $tY + $dY*9
			, "Projection:"
			, $clr
			) ;
	$img->string	( gdMediumBoldFont
			, $tX
			, $tY + $dY*10
			, " UTM Zone 15N"
			, $clr
			) ;

}