set role gis;

drop schema if exists hex cascade;
create schema hex;

set search_path to hex,public,gis;

create or replace function hexgrid(
	xdim integer,
	ydim integer,
	angle numeric,
	coverage geometry,
	srid	integer,
	origin	geometry(POINT),
	tilesize numeric
) returns table (x integer, y integer, cantor numeric, hex geometry) as $$
declare
	rcover geometry;
	stride double precision;
	x integer;
	y integer;
	xmin double precision;
	ymin double precision;
	xoff double precision;
begin
	-- the origin will be the center of the 1,1 hex.
	-- if null, then the origin will be set to fit in the
	-- covering box.

	-- xdim, and ydim can be null, in which case they are
	-- calculated from the tilesize and the covering box
	
	-- calculate ydim

	-- precedence: ydim, xdim, tilesize

	-- counter rotate geometry, fill, rotate grid
	rcover = ST_Rotate(st_transform(coverage,srid), -angle*3.1415926/180.0);
	
	xmin := ST_XMin(rcover);
	ymin := ST_YMin(rcover);
	
	stride := tilesize * 2.0 / sqrt(3.0);

	--raise notice 'stride = %', stride;
	xdim := ((ST_XMax(rcover) -  ST_XMin(rcover))/tilesize*2/sqrt(3.0))::integer;
	ydim := ((ST_YMax(rcover) -  ST_YMin(rcover))/tilesize)::integer;

	--raise notice 'coverage = % % % %', xmin, ymin, xdim, ydim;
	
        -- TODO reorient these so right hand rule is inside
        hex := ST_Scale(ST_GeomFromText('POLYGON((
.577350269189625764509148780502 0.0,
.288675134594812882254574390251 0.5,
 -.288675134594812882254574390251 0.5,
 -.577350269189625764509148780502 0.0,
 -.288675134594812882254574390251 -0.5,
 .288675134594812882254574390251 -0.5,
 .577350269189625764509148780502 0.0))',ST_SRID(rcover)),tilesize,tilesize,1.0);

        xoff := .577350269189625764509148780502 * 1.5;

        for x in 0 .. xdim loop
                for y in 0 .. ydim loop
			-- TODO actually calculate cantor
			return query select x, ydim-y, hex.cantor(ARRAY[x,ydim-y]), 
				ST_Translate(hex,
				tilesize * (x * xoff) + xmin,
				--tilesize * (y + 0.5 * (x%2) - (x+1)/2) + ymin
				tilesize * (y + 0.5 * (x%2) ) + ymin
				)
				--st_rotate(ST_Translate(hex,
				--tilesize * (x * xoff) + xmin,
				--tilesize * (y + 0.5 * (x%2) - (x+1)/2) + ymin
				--), -3.1415926/6)
			;
                end loop;
        end loop;
END;
$$
language 'plpgsql'
;

-- create some hexagon ops

create function ijk(xy integer[]) returns integer[] language 'plpgsql' as $$
declare
	ijk integer[];
begin
	if array_length(xy,1) = 3 then
		return xy;
	end if;

	ijk[1] := xy[1];
	ijk[2] := -xy[2];

	if xy[1] < 0 then
		ijk[2] := ijk[2] + (-xy[1]+1) / 2;
	else
		ijk[2] := ijk[2] - xy[1]/2;
	end if;

	ijk[3] := -ijk[1] - ijk[2];

	return ijk;
end
$$ immutable strict;

create function ijk(x integer, y integer) returns integer[] language 'sql' as $$
	select hex.ijk(ARRAY[x,y]);
$$ immutable strict;

-- invert the cantor paring
create function xy(hex numeric) returns integer[] language 'plpgsql' as $$
declare
	xy integer[];
	w integer;
	t integer;
	py integer;
begin
	if hex >= 0 then
		hex := hex-1;
		w := floor((sqrt(8.0 * hex + 1.0) - 1.0)/2.0)::integer;
		t := (w * w + w)/2;
		py := hex - t;

		xy = ARRAY[w-py,py];
	-- handle the negative cases
	else 
		hex := -hex;
		hex := hex-1;
		w := floor((sqrt(8.0 * hex + 1.0) - 1.0)/2.0)::integer;
		t := (w * w + w)/2;
		py := hex - t;

		xy = ARRAY[w-py,py];

		if xy[1] % 2 = 1 then
			xy[1] = -(xy[1]+1)/2;
		else
			xy[1] = xy[1]/2;
		end if;
		if xy[2] % 2 = 1 then
			xy[2] = -(xy[2]+1)/2;
		else
			xy[2] = xy[2]/2;
		end if;
	end if;

	return xy;
end;
$$ immutable strict;

create function ijk(hex numeric) returns integer[] language 'sql' as $$
	select hex.ijk(hex.xy(hex));
$$ immutable strict;

create function xy(ijk integer[]) returns integer[] language 'plpgsql' as $$
declare
	xy integer[];
begin
	-- make sure this is really an ijk
	-- just no-op if it's already an
	-- xy, so this is sort of a coercion function
	if array_length(ijk,1) = 2 then
		return ijk;
	end if;
	xy[1] := ijk[1];
	xy[2] := -ijk[2];

	if ijk[1] < 0 then
		xy[2] := xy[2] + (-ijk[1] + 1) / 2;
	else
		xy[2] := xy[2] - ijk[1]/2;
	end if;
	return xy;
end;
$$;

CREATE OR REPLACE FUNCTION vec_add(arr1 integer[], arr2 integer[])
RETURNS integer[] AS
$$
SELECT array_agg(result)
FROM (SELECT tuple.val1 + tuple.val2 AS result
      FROM (SELECT UNNEST($1) AS val1
                   ,UNNEST($2) AS val2
                   ,generate_subscripts($1, 1) AS ix) tuple
      ORDER BY ix) inn;
$$ LANGUAGE SQL IMMUTABLE STRICT;

create function adjhex(dir integer, hex integer[]) returns integer[] language 'plpgsql' as $$
declare
	adj integer[];
begin
	if array_length(hex,1) = 2 then
		hex = hex.ijk(hex);
	end if;

	dir := dir % 6;
	if dir = 2 then
		hex := hex.vec_add(hex, ARRAY[-1, 1, 0]);
	elsif dir = 1 then
		hex := hex.vec_add(hex, ARRAY[ 0, 1,-1]);
	elsif dir = 0 then
		hex := hex.vec_add(hex, ARRAY[ 1, 0,-1]);
	elsif dir = 5 then
		hex := hex.vec_add(hex, ARRAY[ 1,-1, 0]);
	elsif dir = 4 then
		hex := hex.vec_add(hex, ARRAY[ 0,-1, 1]);
	elsif dir = 3 then
		hex := hex.vec_add(hex, ARRAY[-1, 0, 1]);
	else
		-- error
	end if;
	return hex;
end;
$$ immutable strict;

create function adjhex(dir integer, hex numeric) returns integer[] language 'sql' as $$
	select hex.adjhex(dir, hex.ijk(hex));
$$ immutable strict;

-- TODO make this work for any, or inline it
create function natcantor(hex integer[]) returns numeric language 'sql' as $$
	select cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
$$ immutable strict;

create function cantor(hex integer[]) returns numeric language 'plpgsql' as $$
begin
	-- convert to xy
	if array_length(hex,1) = 3 then
		hex = hex.xy(hex);
	end if;

	if hex[1] < 0 or hex[2] < 0 then
		hex[1] = hex[1] * 2;
		hex[2] = hex[2] * 2;
		if hex[1] < 0 then
			hex[1] = -hex[1] - 1;
		end if;
		if hex[2] < 0 then
			hex[2] = -hex[2] - 1;
		end if;
		return -cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
	end if;

	return cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
end;
$$ immutable strict;

create function adjacent(hex numeric) returns setof numeric language 'sql' as $$
select hex.cantor(hex.adjhex(dir, hex)) from generate_series(0,5) as dir;
$$ immutable strict;

create function adjacent(dir integer, hex numeric) returns numeric language 'sql' as $$
	select hex.cantor(hex.adjhex(dir, hex.ijk(hex)));
$$ immutable strict;

create function within(hex numeric, range integer = 1, include boolean = false) returns setof numeric language 'plpgsql' as $$
declare
	c integer[]; -- cantor coordinates
	base integer[];
begin
	if range <= 0 then
		return;
	end if;
	if include then
		return next hex;
	end if;

	-- we could in theory eliminate the inner loop and
	-- do a return query
	base = hex.ijk(hex);
	for r in 1 .. range loop
		c = base;
		c[1] := c[1] + r;
		c[3] := -c[1] - c[2]; -- keep the invariant
		hex := hex.cantor(c);
		
		for h in 0 .. r * 6 - 1 loop
			return next hex;
			hex = hex.adjacent(h/r+2,hex);
		end loop;
	end loop;
	return;
end;
$$ immutable strict;

create function within_array(hex numeric, range integer = 1, include boolean = false) returns numeric[] language 'sql' as $$
select array_agg(c) from hex.within(hex, range, include) c;
$$ immutable strict;

create function atrange(hex numeric, range integer = 1) returns setof numeric language 'plpgsql' as $$
declare
	c integer[]; -- cantor coordinates
begin
	if range = 0 then
		return next hex;
	end if;
	if range <= 0 then
		return;
	end if;

	c = ijk(hex);
	c[1] := c[1] + range;
	c[3] := -c[1] - c[2]; -- keep the invariant
	hex := hex.cantor(c);
	
	for q in 0 .. range * 6 - 1 loop
		return next hex;
		hex = hex.adjacent(q/range+2,hex);
	end loop;
	return;
end;
$$ immutable strict;

select
xy(ARRAY[1,2]),
ijk(xy(ARRAY[1,2])),
xy(ijk(xy(ARRAY[1,2])))
;

select xy(ARRAY[2,2]), dir, xy(adjhex(dir, ijk(ARRAY[2,2]))) from generate_series(0,5) as dir;
select
xy(adjhex(dir, ijk(ARRAY[1,2]))),
cantor(xy(adjhex(dir, ijk(ARRAY[1,2])))),
xy(cantor(xy(adjhex(dir, ijk(ARRAY[1,2])))))
from generate_series(0,5) as dir;

select cantor(ARRAY[4,3]);
--select * from within(cantor(ARRAY[4,4]), 2);
select h, xy(h) from within(cantor(ARRAY[4,3]), 2) as h;